text:
# -*- coding: utf-8 -*-
"""
Reimplementation of matlab algorithms for fishlength detection in python.
The next step is to move them into kwiver.

TODO:
    fix hard coded paths for doctests
"""
from __future__ import division, print_function
from collections import namedtuple
import cv2
import itertools as it
import numpy as np
import scipy.optimize
from .imutils import (imscale, ensure_grayscale, from_homog, to_homog)
from os.path import splitext
import ubelt as ub
import logging
import warnings
from six.moves import zip

logger = logging.getLogger(__name__)

OrientedBBox = namedtuple('OrientedBBox', ('center', 'extent', 'angle'))

__OPENCV_VERSION_2__ = cv2.__version__.startswith('2')
if __OPENCV_VERSION_2__:
    warnings.warn('Using an old OpenCV version. '
                  'This may cause unexpected results. '
                  'Please update to 3.x')
    import skimage  # NOQA
    import skimage.measure  # NOQA

def dict_subset(dict_, keys):
    return {k: dict_[k] for k in keys}

def dict_update_subset(dict_, other):
    """
    updates items in `dict_` based on `other`. `other` is not allowed to
    specify any keys that do not already exist in `dict_`.
    """
    for k, v in other.items():
        if k not in dict_:
            raise KeyError(k)
        dict_[k] = v

class ParamInfo(object):
    def __init__(self, name, default, doc=None):
        self.name = name
        self.default = default
        self.doc = doc

class BoundingBox(ub.NiceRepr):
    def __init__(bbox, coords):
        bbox.coords = coords

    def __nice__(self):
        return 'center={}, wh={}'.format(self.center, (self.width, self.height))

    @classmethod
    def from_coords(self, xmin, ymin, xmax, ymax):
        coords = np.array([xmin, ymin, xmax, ymax])
        return BoundingBox(coords)

    @property
    def xmin(bbox):
        return bbox.coords[0]

    @property
    def ymin(bbox):
        return bbox.coords[1]

    @property
    def xmax(bbox):
        return bbox.coords[2]

    @property
    def ymax(bbox):
        return bbox.coords[3]

    @property
    def width(bbox):
        return bbox.xmax - bbox.xmin

    @property
    def height(bbox):
        return bbox.ymax - bbox.ymin

    @property
    def center(self):
        if self.coords is None:
            return None
        (xmin, ymin, xmax, ymax) = self.coords
        cx = (xmax + xmin) / 2
        cy = (ymax + ymin) / 2
        return cx, cy

    def scale(self, factor):
        """
        inplace upscaling of bounding boxes and points
        (masks are not upscaled)
        """
        self.coords = np.array(self.coords) * factor

class DetectedObject(ub.NiceRepr):
    """
    Example:
        >>> from viame.processes.camtrawl.algos import *
        >>> cc_mask = np.zeros((11, 11), dtype=np.uint8)
        >>> cc_mask[3:5, 2:7] = 1
        >>> self = DetectedObject.from_connected_component(cc_mask)
        >>> print(self)
        <DetectedObject(center=(4.0, 3.5), wh=(4, 1))>
    """

    def __init__(self, bbox, mask):
        # bbox is kept in the image coordinate frame
        self.bbox = bbox
        # mask is kept in its own coordinate frame
        self.mask = mask
        # keep track of the scale factor from mask to bbox
        self.bbox_factor = 1.0

    def __nice__(self):
        return self.bbox.__nice__()

    def num_pixels(self):
        """
        Example:
            >>> from viame.processes.camtrawl.algos import *
            >>> cc_mask = np.zeros((11, 11), dtype=np.uint8)
            >>> cc_mask[3:5, 2:7] = 1
            >>> self = DetectedObject.from_connected_component(cc_mask)
            >>> print(str(self.num_pixels()))
            10
        """
        # number of pixels in the mask
        # scale to the number that would be in the original image
        n_pixels = int((self.mask > 0).sum() * (self.bbox_factor ** 2))
        return n_pixels

    def hull(self):
        """
        Example:
            >>> from viame.processes.camtrawl.algos import *
            >>> cc_mask = np.zeros((11, 11), dtype=np.uint8)
            >>> cc_mask[3:5, 2:7] = 1
            >>> self = DetectedObject.from_connected_component(cc_mask)
            >>> print(self.hull().tolist())
            [[[6, 4]], [[2, 4]], [[2, 3]], [[6, 3]]]
        """
        if self.mask.ndim > 2:
            self.mask = np.squeeze(self.mask, axis=2)
        cc_y, cc_x = np.where(self.mask)
        points = np.vstack([cc_x, cc_y]).T
        # Find a minimum oriented bounding box around the points
        hull = cv2.convexHull(points)
        # move points from mask coordinates to image coordinates
        if self.bbox_factor != 1.0:
            hull = hull * self.bbox_factor
        hull = hull + [[self.bbox.xmin, self.bbox.ymin]]
        hull = np.round(hull).astype(np.int)
        return hull

    def oriented_bbox(self):
        """
        Example:
            >>> from viame.processes.camtrawl.algos import *
            >>> cc_mask = np.zeros((11, 11), dtype=np.uint8)
            >>> cc_mask[3:5, 2:7] = 1
            >>> self = DetectedObject.from_connected_component(cc_mask)
            >>> print(str(self.oriented_bbox()))
            OrientedBBox(center=(4.0, 3.5), extent=(1.0, 4.0), angle=-90.0)
        """
        hull = self.hull()
        oriented_bbox = OrientedBBox(*cv2.minAreaRect(hull))
        return oriented_bbox

    def box_points(self):
        """
        CommandLine:
            python -m viame.processes.camtrawl.algos DetectedObject.box_points

        Example:
            >>> from viame.processes.camtrawl.algos import *
            >>> _, B = 0, 1
            >>> cc_mask = np.array([
            >>>     [_, _, _, _, _, _, _, _],
            >>>     [_, _, _, B, _, _, _, _],
            >>>     [_, _, B, B, B, _, _, _],
            >>>     [_, B, B, B, B, B, _, _],
            >>>     [_, _, B, B, B, B, B, _],
            >>>     [_, _, _, B, B, B, _, _],
            >>>     [_, _, _, _, B, _, _, _],
            >>>     [_, _, _, _, _, _, _, _],
            >>> ])
            >>> self = DetectedObject.from_connected_component(cc_mask)
            >>> points = self.box_points()
            >>> print(ub.repr2(points.tolist(), precision=2, nl=0))
            [[4.00, 6.00], [1.00, 3.00], [3.00, 1.00], [6.00, 4.00]]
        """
        if __OPENCV_VERSION_2__:
            return np.array(cv2.cv.BoxPoints(self.oriented_bbox()),
                            dtype=np.float32)
        else:
            return cv2.boxPoints(self.oriented_bbox())

    def scale(self, factor):
        """ inplace """
        if factor != 1.0:
            self.bbox_factor *= factor
            self.bbox.scale(factor)

    @classmethod
    def from_connected_component(DetectedObject, cc_mask):
        # note, `np.where` returns coords in (r, c)
        ys, xs = np.where(cc_mask)
        xmin, xmax = xs.min(), xs.max()
        ymin, ymax = ys.min(), ys.max()
        bbox = BoundingBox.from_coords(xmin, ymin, xmax, ymax)
        yslice = slice(bbox.ymin, bbox.ymax + 1)
        xslice = slice(bbox.xmin, bbox.xmax + 1)
        mask = cc_mask[yslice, xslice]
        self = DetectedObject(bbox, mask)
        return self

class GMMForegroundObjectDetector(object):
    """
    Uses background subtraction and 4-way connected compoments algorithm to
    detect potential fish objects. Objects are filtered by size, aspect ratio,
    and closeness to the image border to remove bad detections.

    References:
        https://stackoverflow.com/questions/37300698/gaussian-mixture-model
        http://docs.opencv.org/trunk/db/d5c/tutorial_py_bg_subtraction.html
    """

    @staticmethod
    def default_params():
        default_params = {
            'preproc': [

                ParamInfo(name='factor', default=1.0,
                          doc='image downsample factor'),

                ParamInfo(name='smooth_ksize', default=(10, 10),
                          doc=('postprocessing filter size for noise removal '
                               '(wrt orig image size)'))

            ],

            'gmm': [

                ParamInfo(name='n_startup_frames', default=3,
                          doc=('number of frames before the background model '
                               'will wait before returning any detections')),

                ParamInfo(name='n_training_frames', default=300,
                          doc='number of frames to use for training'),

                ParamInfo(name='gmm_thresh', default=30,
                          doc='GMM variance threshold'),

            ],

            'filter': DetectionShapeFilter.default_params()['filter'],
        }
        return default_params

    def __init__(detector, **kwargs):
        # setup default config
        detector.config = {}
        default_params = detector.default_params()
        for pinfos in default_params.values():
            detector.config.update({pi.name: pi.default for pi in pinfos})
        # modify based on user args
        dict_update_subset(detector.config, kwargs)

        # Setup GMM background subtraction algorithm
        logger.debug('Using GMM from cv2.__version__ = {}'.format(cv2.__version__))
        if cv2.__version__.startswith('2'):
            # not sure about these params
            detector.background_model = cv2.BackgroundSubtractorMOG2(
                history=detector.config['n_training_frames'],
                varThreshold=detector.config['gmm_thresh'],
                bShadowDetection=False
            )
        else:
            detector.background_model = cv2.createBackgroundSubtractorMOG2(
                history=detector.config['n_training_frames'],
                varThreshold=detector.config['gmm_thresh'],
                detectShadows=False)

        # Setup detection filter algorithm
        filter_config = dict_subset(
            detector.config, [pi.name for pi in default_params['filter']])
        detector.filter = DetectionShapeFilter(**filter_config)

        detector.n_iters = 0
        # masks from previous iter are kept in memory for visualization
        detector._masks = {}

    def detect(detector, img):
        """
        Main algorithm step.
        Detects the objects in the image and update the background model.

        Args:
            img (ndarray): image to perform detection on

        Returns:
            detections : list of DetectedObjects

        CommandLine:
            python -m viame.processes.camtrawl.algos GMMForegroundObjectDetector.detect

        Doctest:
            >>> import matplotlib as mpl
            >>> mpl.use('agg')
            >>> from viame.processes.camtrawl.algos import *
            >>> from viame.processes.camtrawl.demo import *
            >>> detector, img = demodata_detections(target_step='detect', target_frame_num=7)
            >>> detections = detector.detect(img)
            >>> print('detections = {!r}'.format(detections))
            >>> masks = detector._masks
            >>> # xdoc: REQUIRES(--show)
            >>> draw_img = DrawHelper.draw_detections(img, detections, masks)
            >>> fpath = ub.ensure_app_cache_dir('camtrawl') + '/GMMForegroundObjectDetector.detect.png'
            >>> cv2.imwrite(fpath, draw_img)
            >>> ub.startfile(fpath)
            >>> #from matplotlib import pyplot as plt
            >>> #plt.imshow(cv2.cvtColor(draw_img, cv2.COLOR_BGR2RGB))
            >>> #plt.gca().grid(False)
            >>> #plt.show()
        """

        detector._masks = {}

        # Downsample and convert to grayscale
        img_, upfactor = detector.preprocess_image(img)

        # Run detection / update background model
        logger.debug('detect background')
        mask = detector.background_model.apply(img_)
        detector._masks['orig'] = mask.copy()

        if detector.n_iters < detector.config['n_startup_frames']:
            # Skip the first few frames while the model is learning
            detections = []
        else:
            # Remove noise
            if detector.config['smooth_ksize'] is not None:
                mask = detector.postprocess_mask(mask)
                detector._masks['post'] = mask.copy()

            # Find detections using CC algorithm
            detectgen = detector.detections_in_mask(mask)

            # Filter detections by shape and size
            img_dsize = tuple(img.shape[0:2][::-1])
            detectgen = detector.filter.filter_detections(detectgen, img_dsize)

            # return detections in a list
            detections = list(detectgen)

        detector.n_iters += 1
        logger.debug('made {} detections'.format(len(detections)))
        return detections

    def preprocess_image(detector, img):
        """
        Preprocess image before subtracting backround
        """
        logger.debug('preprocess image before detect')
        # Convert to grayscale
        img_ = ensure_grayscale(img)
        # Downsample image before running detection
        factor = detector.config['factor']
        if factor != 1.0:
            downfactor_ = 1 / factor
            img_, downfactor = imscale(img, downfactor_)
            upfactor = 1 / downfactor[0]
        else:
            upfactor = 1.0
            img_ = img
        return img_, upfactor

    def postprocess_mask(detector, mask):
        """ remove noise from detection intensity masks """
        logger.debug('postprocess mask')
        ksize = np.array(detector.config['smooth_ksize'])
        ksize = tuple(np.round(ksize / detector.config['factor']).astype(np.int))
        kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, ksize)
        # opening is erosion followed by dilation
        mask = cv2.morphologyEx(mask, cv2.MORPH_OPEN, kernel, dst=mask)
        # Do a second dilation
        mask = cv2.dilate(src=mask, kernel=kernel, dst=mask)
        return mask

    def detections_in_mask(detector, mask):
        """
        Find pixel locs of each cc and determine if its a valid detection

        Args:
            mask (ndarray): mask where non-zero pixels indicate all candidate
                objects

        Example:
            >>> from viame.processes.camtrawl.algos import *
            >>> detector = GMMForegroundObjectDetector()
            >>> detector.config['min_num_pixels'] = 2
            >>> x, y = np.indices((10, 10))
            >>> s = 2
            >>> mask = (((x // s) % s == 0) & ((y // s) % s == 0)).astype(np.uint8)
            >>> mask[0, -2:] = 0
            >>> mask[1, -1] = 0
            >>> mask[2, 2] = 1
            >>> mask[-1, 2] = 1
            >>> mask[-5:-1, -3] = 1
            >>> detections = list(detector.detections_in_mask(mask))
            >>> assert len(detections) == 7
        """
        # 4-way connected compoment algorithm
        if __OPENCV_VERSION_2__:
            # opencv2 doesnt have a builtin CC algo, need to use skimage
            cc_mask, n_ccs_sk = skimage.measure.label(mask, neighbors=8,
                                                      background=0,
                                                      return_num=True)
            # Be consistent with opencv, which always includes the background
            # label in the num (even if no background exists).
            n_ccs = n_ccs_sk + 1
        else:
            n_ccs, cc_mask = cv2.connectedComponents(mask, connectivity=8)

        factor = detector.config['factor']

        # Process only labels with enough points
        min_num_pixels = detector.config['min_num_pixels']
        if min_num_pixels is None:
            valid_labels = np.arange(1, n_ccs + 1)
        else:
            # speed optimization: quickly determine num pixels for each cc
            # using a histogram instead of checking in the filter func
            hist, bins = np.histogram(cc_mask[cc_mask > 0].ravel(),
                                      bins=np.arange(1, n_ccs + 1))
            min_num_pixels_ = min_num_pixels / (factor ** 2)
            # only consider large enough regions
            valid_labels = bins[0:-1][hist >= min_num_pixels_]

        # Filter ccs to generate only "good" detections
        # We may be able to speed this up using ndimage.find_objects
        for cc_label in valid_labels:
            cc = (cc_mask == cc_label)
            detection = DetectedObject.from_connected_component(cc)
            # Upscale back to input img coords (to agree with camera calib)
            detection.scale(factor)
            yield detection

class DetectionShapeFilter(object):
    """
    Filters masked detections based on their size and shape
    """

    @staticmethod
    def default_params():
        default_params = {
            'filter': [

                ParamInfo(name='min_num_pixels', default=800,
                          doc=('keep detections only with the number of pixels '
                               'wrt original image size')),

                ParamInfo(name='edge_trim', default=(12, 12),
                          doc=('limits accepable targets to be within the region '
                               '[padx, pady, img_w - padx, img_h - pady]. '
                               'These are wrt the original image size')),

                ParamInfo(name='aspect_thresh', default=(3.5, 7.5,),
                          doc='range of valid aspect ratios for detections')

            ]
        }
        return default_params

    def __init__(self, **kwargs):
        # setup default config
        self.config = {}
        default_params = self.default_params()
        for pinfos in default_params.values():
            self.config.update({pi.name: pi.default for pi in pinfos})
        # modify based on user args
        dict_update_subset(self.config, kwargs)

    def filter_detections(self, detections, img_dsize=None):
        logger.debug('filter detections')
        for detection in detections:
            if self.is_valid(detection, img_dsize):
                yield detection

    def is_valid(self, detection, img_dsize=None):
        """
        Checks if the detection passes filtering constraints

        Args:
            detection (DetectedObject): mask where non-zero pixels indicate a single
                candidate object
            img_dsize (tuple): w/h of the original image

        Returns:
            bool: True if the detection is valid else False
        """
        if self.config['min_num_pixels'] is not None:
            # Remove small regions
            if detection.num_pixels() < self.config['min_num_pixels']:
                return False

        if self.config['edge_trim'] is not None:
            # Define thresholds to filter edges
            img_width, img_height = img_dsize
            xmin_lim, ymin_lim = self.config['edge_trim']
            xmax_lim = img_width - (xmin_lim)
            ymax_lim = img_height - (ymin_lim)

            # Filter objects detected on the edge of the image region
            (xmin, ymin, xmax, ymax) = detection.bbox.coords
            if any([xmin < xmin_lim, xmax > xmax_lim,
                    ymin < ymin_lim, ymax > ymax_lim]):
                return None

        # Find a minimum oriented bounding box around the points
        w, h = detection.oriented_bbox().extent
        if w == 0 or h == 0:
            # assert w == 0 or h == 0, 'detection has no width or height'
            return False

        # Filter objects without fishy aspect ratios
        ar = max([(w / h), (h / w)])
        min_aspect, max_aspect = self.config['aspect_thresh']
        if any([ar < min_aspect, ar > max_aspect]):
            return False
        return True

class FishStereoMeasurments(object):
    """
    Algo for matching detections in left and right camera and determening the
    fish length in milimeters.
    """

    @staticmethod
    def default_params():
        default_params = {
            'thresholds': [
                ParamInfo('max_err', (6, 14), doc=(
                    'Threshold for errors between before & after projected '
                    'points to make matches between left and right')),

                ParamInfo('small_len', 150, doc=(
                    'length (in milimeters) to switch between high and low '
                    'max error thresholds ')),
            ]
        }
        return default_params

    def __init__(self, **kwargs):
        # setup default config
        self.config = {}
        default_params = self.default_params()
        for pinfos in default_params.values():
            self.config.update({pi.name: pi.default for pi in pinfos})
        # modify based on user args
        dict_update_subset(self.config, kwargs)

    def triangulate(self, cal, det1, det2):
        """
        Assuming, det1 matches det2, we determine 3d-coordinates of each
        detection and measure the reprojection error.

        References:
            http://answers.opencv.org/question/117141
            https://gist.github.com/royshil/7087bc2560c581d443bc
            https://stackoverflow.com/a/29820184/887074

        Doctest:
            >>> # Rows are detections in img1, cols are detections in img2
            >>> from viame.processes.camtrawl.algos import *
            >>> from viame.processes.camtrawl.demo import *
            >>> detections1, detections2, cal = demodata_detections(target_step='triangulate', target_frame_num=6)
            >>> det1, det2 = detections1[0], detections2[0]
            >>> self = FishStereoMeasurments()
            >>> assignment, assign_data, cand_errors = self.triangulate(cal, det1, det2)
        """
        logger.debug('triangulate')
        _debug = 0
        if _debug:
            # Use 4 corners and center to ensure matrix math is good
            # (hard to debug when ndims == npts, so make npts >> ndims)
            pts1 = np.vstack([det1.box_points(), det1.oriented_bbox().center])
            pts2 = np.vstack([det2.box_points(), det2.oriented_bbox().center])
        else:
            # Use only the corners of the bbox
            pts1 = det1.box_points()[[0, 2]]
            pts2 = det2.box_points()[[0, 2]]

        # Move into opencv point format (num x 1 x dim)
        pts1_cv = pts1[:, None, :]
        pts2_cv = pts2[:, None, :]

        # Grab camera parameters
        K1, K2 = cal.intrinsic_matrices()
        kc1, kc2 = cal.distortions()
        rvec1, tvec1, rvec2, tvec2 = cal.extrinsic_vecs()

        # Make extrincic matrices
        R1 = cv2.Rodrigues(rvec1)[0]
        R2 = cv2.Rodrigues(rvec2)[0]
        T1 = tvec1[:, None]
        T2 = tvec2[:, None]
        RT1 = np.hstack([R1, T1])
        RT2 = np.hstack([R2, T2])

        # Undistort points
        # This puts points in "normalized camera coordinates" making them
        # independent of the intrinsic parameters. Moving to world coordinates
        # can now be done using only the RT transform.
        unpts1_cv = cv2.undistortPoints(pts1_cv, K1, kc1)
        unpts2_cv = cv2.undistortPoints(pts2_cv, K2, kc2)

        # note: trinagulatePoints docs say that it wants a 3x4 projection
        # matrix (ie K.dot(RT)), but we only need to use the RT extrinsic
        # matrix because the undistorted points already account for the K
        # intrinsic matrix.
        world_pts_homog = cv2.triangulatePoints(RT1, RT2, unpts1_cv, unpts2_cv)
        world_pts = from_homog(world_pts_homog)

        # Compute distance between 3D bounding box points
        if _debug:
            corner1, corner2 = world_pts.T[[0, 2]]
        else:
            corner1, corner2 = world_pts.T

        # Length is in milimeters
        fishlen = np.linalg.norm(corner1 - corner2)

        # Reproject points
        world_pts_cv = world_pts.T[:, None, :]
        proj_pts1_cv = cv2.projectPoints(world_pts_cv, rvec1, tvec1, K1, kc1)[0]
        proj_pts2_cv = cv2.projectPoints(world_pts_cv, rvec2, tvec2, K2, kc2)[0]

        # Check error
        err1 = ((proj_pts1_cv - pts1_cv)[:, 0, :] ** 2).sum(axis=1)
        err2 = ((proj_pts2_cv - pts2_cv)[:, 0, :] ** 2).sum(axis=1)
        errors = np.hstack([err1, err2])

        # Get 3d points in each camera's reference frame
        # Note RT1 is the identity and RT are 3x4, so no need for `from_homog`
        # Return points in with shape (N,3)
        pts1_3d = RT1.dot(to_homog(world_pts)).T
        pts2_3d = RT2.dot(to_homog(world_pts)).T
        return pts1_3d, pts2_3d, errors, fishlen

    @staticmethod
    def minimum_weight_assignment(cost_errors):
        """
        Finds optimal assignment of left-camera to right-camera detections

        Doctest:
            >>> # Rows are detections in img1, cols are detections in img2
            >>> from viame.processes.camtrawl.algos import *
            >>> self = FishStereoMeasurments()
            >>> cost_errors = np.array([
            >>>     [9, 2, 1, 9],
            >>>     [4, 1, 5, 5],
            >>>     [9, 9, 2, 4],
            >>> ])
            >>> assign1 = self.minimum_weight_assignment(cost_errors)
            >>> assign2 = self.minimum_weight_assignment(cost_errors.T)
        """
        n1, n2 = cost_errors.shape
        n = max(n1, n2)
        # Embed the [n1 x n2] matrix in a padded (with inf) [n x n] matrix
        cost_matrix = np.full((n, n), fill_value=np.inf)
        cost_matrix[0:n1, 0:n2] = cost_errors

        # Find an effective infinite value for infeasible assignments
        is_infeasible = np.isinf(cost_matrix)
        is_positive = cost_matrix > 0
        feasible_vals = cost_matrix[~(is_infeasible & is_positive)]
        large_val = (n + feasible_vals.sum()) * 2
        # replace infinite values with effective infinite values
        cost_matrix[is_infeasible] = large_val

        # Solve munkres problem for minimum weight assignment
        indexes = list(zip(*scipy.optimize.linear_sum_assignment(cost_matrix)))
        # Return only the feasible assignments
        assignment = [(i, j) for (i, j) in indexes
                      if cost_matrix[i, j] < large_val]
        return assignment

    def find_matches(self, cal, detections1, detections2):
        """
        Match detections from the left camera to detections in the right camera

        Doctest:
            >>> # Rows are detections in img1, cols are detections in img2
            >>> from viame.processes.camtrawl.algos import *
            >>> detections1, detections2, cal = demodata_detections(target_step='triangulate', target_frame_num=6)
            >>> self = FishStereoMeasurments()
            >>> assignment, assign_data, cand_errors = self.find_matches(cal, detections1, detections2)
        """
        logger.debug('find matches')
        n_detect1, n_detect2 = len(detections1), len(detections2)
        cand_data = {}

        # Initialize matrix of reprojection errors
        cost_errors = np.full((n_detect1, n_detect2), fill_value=np.inf)
        cand_errors = np.full((n_detect1, n_detect2), fill_value=np.inf)

        # Find the liklihood that each pair of detections matches by
        # triangulating and then measuring the reprojection error.
        for (i, det1), (j, det2) in it.product(enumerate(detections1),
                                               enumerate(detections2)):
            # Triangulate assuming det1 and det2 match, but return the
            # reprojection error so we can check if this assumption holds
            pts1_3d, pts2_3d, errors, fishlen = self.triangulate(cal, det1, det2)
            error = errors.mean()

            # Mark the pair (i, j) as a potential candidate match
            cand_errors[i, j] = error

            # Record information about this candidate match
            cand_data[(i, j)] = {
                'world_pts': pts1_3d,
                'error': error,
                'fishlen': fishlen,
                'range': pts1_3d.T[2].mean(),
                'dz': np.abs(np.diff(pts1_3d.T[2]))[0],
                'box_pts1': det1.box_points(),
                'box_pts2': det2.box_points(),
            }

            # Check chirality
            # Both Z-coordinates must be positive (i.e. in front the cameras)
            z_coords1 = pts1_3d.T[2]
            z_coords2 = pts2_3d.T[2]
            both_in_front = np.all(z_coords1 > 0) and np.all(z_coords2 > 0)
            if not both_in_front:
                # Ignore out-of-view correspondences
                continue

            # Check if reprojection error is too high
            max_error = self.config['max_err']
            small_len = self.config['small_len']  # hardcoded to 15cm in matlab
            if len(max_error) == 2:
                error_thresh = max_error[0] if fishlen <= small_len else max_error[1]
            else:
                error_thresh = max_error[0]

            if error  >= error_thresh:
                # Ignore correspondences with high reprojection error
                continue

            cost_errors[i, j] = error

        # Find the matching with minimum reprojection error, such that each
        # detection in one camera can match at most one detection in the other.
        assignment = self.minimum_weight_assignment(cost_errors)

        # get associated data with each assignment
        assign_data = []
        for i, j in assignment:
            data = {'ij': (i, j)}
            data.update(cand_data[(i, j)])
            assign_data.append(data)
        return assignment, assign_data, cand_errors

class StereoCalibration(object):
    """
    Helper class for reading / accessing stereo camera calibration params

    Doctest:
        >>> from viame.processes.camtrawl.algos import *
        >>> cal_fpath = '/home/joncrall/data/camtrawl_stereo_sample_data/201608_calibration_data/selected/Camtrawl_2016.npz'
        >>> cal = StereoCalibration.from_file(cal_fpath)
    """
    def __init__(cal, data=None):
        cal.data = data
        cal.unit = 'milimeters'

    def __str__(cal):
        return '{}({})'.format(cal.__class__.__name__, cal.data)

    def extrinsic_vecs(cal):
        rvec1 = cal.data['left']['extrinsic']['om']
        tvec1 = cal.data['right']['extrinsic']['om']

        rvec2 = cal.data['right']['extrinsic']['om']
        tvec2 = cal.data['right']['extrinsic']['T']
        return rvec1, tvec1, rvec2, tvec2

    def distortions(cal):
        kc1 = cal.data['right']['intrinsic']['kc']
        kc2 = cal.data['left']['intrinsic']['kc']
        return kc1, kc2

    def intrinsic_matrices(cal):
        K1 = cal._make_intrinsic_matrix(cal.data['left']['intrinsic'])
        K2 = cal._make_intrinsic_matrix(cal.data['right']['intrinsic'])
        return K1, K2

    @staticmethod
    def _make_intrinsic_matrix(intrin):
        """ convert intrinsic dict to matrix """
        fc = intrin['fc']
        cc = intrin['cc']
        alpha_c = intrin['alpha_c']
        KK = np.array([
            [fc[0], alpha_c * fc[0], cc[0]],
            [    0,           fc[1], cc[1]],
            [    0,               0,     1],
        ])
        return KK

    @staticmethod
    def _make_intrinsic_params(K):
        """ convert intrinsic matrix to dict """
        intrin = {}
        fc = intrin['fc'] = np.zeros(2)
        cc = intrin['cc'] = np.zeros(2)
        [[fc[0], alpha_c_fc0, cc[0]],
         [    _,       fc[1], cc[1]],
         [    _,           _,     _]] = K
        intrin['alpha_c'] = np.array([alpha_c_fc0 / fc[0]])
        return intrin

    @classmethod
    def from_file(StereoCalibration, cal_fpath):
        """
        Loads a camera calebration from a .mat or .npz file

        SeeAlso:
            from_npzfile
            from_matfile
        """
        ext = splitext(cal_fpath)[1].lower()
        if ext == '.mat':
            return StereoCalibration.from_matfile(cal_fpath)
        elif ext == '.npz':
            return StereoCalibration.from_npzfile(cal_fpath)
        else:
            raise ValueError('unknown extension {}'.format(ext))

    @classmethod
    def from_npzfile(StereoCalibration, cal_fpath):
        """
        For the npz file the root object should be a dict with the following
        keys and values:
            R: extrinsic rotation matrix
            T: extrinsic translation
            cameraMatrixL: dict of intrinsict parameters for the left camera
                fc: focal length
                cc: principle point
                alpha_c: skew
            cameraMatrixR: dict of intrinsict parameters for the right camera
                fc: focal length
                cc: principle point
                alpha_c: skew
            distCoeffsL: distortion coefficients for the left camera
            distCoeffsR: distortion coefficients for the right camera
        """
        logger.debug('Loading npzfile {}'.format(cal_fpath))
        data = dict(np.load(cal_fpath))
        flat_dict = {}
        flat_dict['om'] = cv2.Rodrigues(data['R'])[0].ravel()
        flat_dict['T'] = data['T'].ravel()

        K1 = data['cameraMatrixL']
        intrin1 = StereoCalibration._make_intrinsic_params(K1)
        flat_dict['fc_left'] = intrin1['fc']
        flat_dict['cc_left'] = intrin1['cc']
        flat_dict['alpha_c_left'] = intrin1['alpha_c']
        flat_dict['kc_left'] = data['distCoeffsL'].ravel()

        K2 = data['cameraMatrixR']
        intrin2 = StereoCalibration._make_intrinsic_params(K2)
        flat_dict['fc_right'] = intrin2['fc']
        flat_dict['cc_right'] = intrin2['cc']
        flat_dict['alpha_c_right'] = intrin2['alpha_c']
        flat_dict['kc_right'] = data['distCoeffsR'].ravel()
        return StereoCalibration._from_flat_dict(flat_dict)

    def from_cameras(StereoCalibration, camera1, camera2):
        pass

    @classmethod
    def from_matfile(StereoCalibration, cal_fpath):
        """
        Loads a matlab camera calibration file from disk

        For the mat file, the root structure should be a dict with the key
        `Cal` whose value is a dict with the following items:
            om: extrinsic rotation vector (note rotation matrix is rodrigues(om))
            T: extrinsic translation
            fc_left: focal length of the left camera
            cc_left: principle point
            alpha_c_left: skew
            kc_left: distortion coefficients for the left camera
            fc_right: focal length of the right camera
            cc_right: principle point
            alpha_c_right: skew
            kc_right: distortion coefficients for the right camera

        References:
            http://www.vision.caltech.edu/bouguetj/calib_doc/htmls/parameters.html
            http://www.vision.caltech.edu/bouguetj/calib_doc/htmls/example5.html

        Doctest:
            >>> from viame.processes.camtrawl.algos import *
            >>> from viame.processes.camtrawl.demo import *
            >>> cal_fpath = '/home/joncrall/data/autoprocess_test_set/cal_201608.mat'
            >>> cal = StereoCalibration.from_matfile(cal_fpath)
            >>> print('cal = {}'.format(cal))
        """
        import scipy.io
        logger.debug('Loading matfile {}'.format(cal_fpath))
        cal_data = scipy.io.loadmat(cal_fpath)
        keys = ['om', 'T', 'fc_left', 'fc_right', 'cc_left', 'cc_right',
                'kc_left', 'kc_right', 'alpha_c_left', 'alpha_c_right']

        if isinstance(cal_data, dict) and 'Cal' in cal_data:
            vals = cal_data['Cal'][0][0]
            flat_dict = {k: v.ravel() for k, v in zip(keys, vals)}
        else:
            flat_dict = {key: cal_data[key].ravel() for key in keys}
        return StereoCalibration._from_flat_dict(flat_dict)

    @classmethod
    def _from_flat_dict(StereoCalibration, flat_dict):
        """ helper used by matlab and numpy readers """
        data = {
            'left': {
                'extrinsic': {
                    # Center wold on the left camera
                    'om': np.zeros(3),  # rotation vector
                    'T': np.zeros(3),  # translation vector
                },
                'intrinsic': {
                    'fc': flat_dict['fc_left'],  # focal point
                    'cc': flat_dict['cc_left'],  # principle point
                    'alpha_c': flat_dict['alpha_c_left'][0],  # skew
                    'kc': flat_dict['kc_left'],  # distortion
                }
            },

            'right': {
                'extrinsic': {
                    'om': flat_dict['om'],  # rotation vector
                    'T': flat_dict['T'],  # translation vector
                },
                'intrinsic': {
                    'fc': flat_dict['fc_right'],  # focal point
                    'cc': flat_dict['cc_right'],  # principle point
                    'alpha_c': flat_dict['alpha_c_right'][0],  # skew
                    'kc': flat_dict['kc_right'],  # distortion
                }
            },
        }
        cal = StereoCalibration()
        cal.data = data
        return cal

# if __name__ == '__main__':
#     from . import demo
#     logging.basicConfig()
#     demo.demo()
#     # import camtrawl_demo
#     # camtrawl_demo.demo()
#     # import utool as ut
#     # ut.dump_profile_text()
#     # import pytest
#     # pytest.main([__file__, '--doctest-modules'])
if __name__ == '__main__':
    r"""
    CommandLine:
        python -m viame.processes.camtrawl.algos

    Ignore:
        workon_py2
        source ~/code/VIAME/build/install/setup_viame.sh
        export SPROKIT_PYTHON_MODULES=camtrawl_processes:kwiver.processes:viame.processes
        export PYTHONPATH=$(pwd):$PYTHONPATH
        python ~/code/VIAME/plugins/camtrawl/python/camtrawl_demo.py
        ffmpeg -y -f image2 -i out_haul83/%*.png -vcodec mpeg4 -vf "setpts=10*PTS" haul83-results.avi
    """
    import xdoctest
    xdoctest.doctest_module(__file__)

reason: import numpy,import scipy


text:
import torch
import numpy as np
import pandas as pd
from torch.autograd import Variable
from sklearn.model_selection import train_test_split
from torch.utils.data import DataLoader,TensorDataset
import matplotlib.pyplot as plt
import torch.nn as nn

path="mnist/"

#first column label , other columns image pixels

#PREPARE THE DATA

train = pd.read_csv(path+"train.csv",dtype="float32")
train_y=train.label.values
train_x=train.drop('label',axis=1).values/255

#splitting data
x_train,x_test,y_train,y_test=train_test_split(train_x,train_y,test_size=0.2,random_state=42)

#array to tensor
x_train=torch.from_numpy(x_train)
x_test=torch.from_numpy(x_test)

y_train=torch.from_numpy(y_train).type(torch.LongTensor)
y_test=torch.from_numpy(y_test).type(torch.LongTensor)

#hyperparameters
BATCH_SIZE=100
n_iters=10001
num_epoch=n_iters/ (len(x_train)/BATCH_SIZE)
num_epoch=int(num_epoch)

#tensor dataset and dataloader

train=TensorDataset(x_train,y_train)
test=TensorDataset(x_test,y_test)

train_loader=DataLoader(train,batch_size=BATCH_SIZE,shuffle=False)
test_loader=DataLoader(test,batch_size=BATCH_SIZE,shuffle=False)

#data visualization
"""
plt.imshow(train_x[10].reshape(28,28))
plt.axis("off")
plt.title(str(train_y[10]))
plt.savefig("graph.png")
plt.show()"""

class ANNModel(nn.Module):
    #init ve forward: initte sadece tanımlanır
    # 784-150 relu-150-tanh 150-elu 10 initilaze etti
    #forward'a gidip birbirine baglanacak

    def __init__(self,input_dim,hidden_dim,output_dim):
        super(ANNModel,self).__init__()
        #linear func 1 : input_dim -> hidden , 784-150
        self.fc1=nn.Linear(input_dim, hidden_dim)
        #nonlinearty 1
        self.relu1=nn.ReLU()

        #linear func 2 : hidden -> hidden , 150-150
        self.fc2=nn.Linear(hidden_dim, hidden_dim)
        #nnreality2
        self.tanh2=nn.Tanh()

        #linear func 3 : hidden -> hidden , 150-150
        self.fc3=nn.Linear(hidden_dim, hidden_dim)
        #nnreality3
        self.elu3=nn.ELU()

        #linear func 3 : hidden -> out , 150-10
        self.fc4=nn.Linear(hidden_dim, output_dim)

    def forward(self, x):

        out=self.fc1(x)
        out=self.relu1(out)
        out=self.fc2(out)
        out=self.tanh2(out)
        out=self.fc3(out)
        out=self.elu3(out)
        out=self.fc4(out)

        return out

input_dim=28*28
hidden_dim=150
output_dim=10

model=ANNModel((input_dim), hidden_dim, output_dim)
error=nn.CrossEntropyLoss()
lr=0.002
optimizer=torch.optim.SGD(params=model.parameters(), lr=lr)

count=0
loss_list=[]
iter_list=[]
acc_list=[]

for epoch in range(num_epoch):
    for i, (images,labels) in enumerate(train_loader):

        train= Variable(images.view(-1,28*28))
        labels=Variable(labels)

        #clear grads
        optimizer.zero_grad()

        #forward prop
        outputs=model(train)

        #loss
        loss=error(outputs,labels)

        #calculate grads
        loss.backward()

        #update params
        optimizer.step()

        count+=1

        if count%50==0:
            correct=0
            total=0

            for images,labels in test_loader:

                test=Variable(images.view(-1,28*28))
                outputs=model(test)
                #print(outputs)
                predicted=torch.max(outputs.data,1)[1]#max prob value
                total+=len(labels)

                correct+=(predicted ==labels).sum()

            acc=100*correct/float(total)

            loss_list.append(loss.data)
            iter_list.append(count)
            acc_list.append(acc)

        if count %100==0:
            print("Iteration {} Loss {} Accuracy {}".format(count,loss.data,acc))

#loss & acc visualization

plt.plot(iter_list,loss_list)
plt.xlabel("Iteratıon")
plt.ylabel("Loss")
plt.title("ANN: Loss vs Number of Iteration")
plt.show()

plt.plot(iter_list,acc_list,color="red")
plt.xlabel("Iteratıon")
plt.ylabel("Accuracy")
plt.title("ANN: Accuracy vs Number of Iteration")

plt.show()

ext: py

max_stars_repo_name: mftnakrsu/basic-tutorial-on-pytorch

max_issues_repo_name: mftnakrsu/basic-tutorial-on-pytorch

max_forks_repo_name: mftnakrsu/basic-tutorial-on-pytorch

reason: import numpy


text:
from flask import Flask
from flask_restful import Resource
from flask_restful import Api

import numpy as np
import cv2

import werkzeug
from flask_restful import reqparse

parser = reqparse.RequestParser()
parser.add_argument("file", type = werkzeug.datastructures.FileStorage, location = "files")

app = Flask(__name__)
api = Api(app)

import base64

class ImageServer(Resource):

    def _decode(self, image):
        image = np.fromstring(image, np.uint8)
        image = cv2.imdecode(image, cv2.IMREAD_UNCHANGED)

        return image

    def _encode(self, image):
        _, encoded_image = cv2.imencode(".jpg", image)
        image_data = base64.b64encode(encoded_image).decode("utf-8")

        return f"data:image/jpeg;base64,{image_data}"

    def post(self):
        data = parser.parse_args()
        image = data["file"].read()
        image = self._decode(image = image)
        image = self._encode(image = image)

        return {"image" : image}

api.add_resource(ImageServer, "/")

if __name__ == "__main__":
    app.run(debug = True, port = 5000)


text:
from datasets.MOT.dataset import MultipleObjectTrackingDatasetSequence_MemoryMapped
from ._algo import do_siamfc_pair_sampling_positive_only, do_siamfc_pair_sampling_negative_only, do_siamfc_pair_sampling
import numpy as np
from data.tracking.sampler._sampling_algos.sequence_sampling.common._dummy_bbox import generate_dummy_bbox_xyxy
from data.operator.bbox.validity import bbox_is_valid
from data.operator.bbox.spatial.utility.half_pixel_offset.image import bounding_box_is_intersect_with_image
from data.tracking.sampler.SiamFC.type import SiamesePairSamplingMethod

def _get_bbox_and_generate_dummy_bbox_if_empty(frame, track_id, rng_engine, reference_bbox):
    if frame.has_object(track_id):
        x_obj_info = frame.get_object_by_id(track_id)
        x_bbox_validity = x_obj_info.get_bounding_box_validity_flag()
        if x_bbox_validity is not None and not x_bbox_validity:
            x_bbox = generate_dummy_bbox_xyxy(frame.get_image_size(), rng_engine, reference_bbox)
        else:
            x_bbox = x_obj_info.get_bounding_box()
    else:
        x_bbox = generate_dummy_bbox_xyxy(frame.get_image_size(), rng_engine, reference_bbox)
    return x_bbox

def _get_positive_frame_image_and_bbox(sequence, index, track_id):
    frame = sequence.get_frame(index)
    frame_image = frame.get_image_path()
    frame_bbox = frame.get_object_by_id(track_id).get_bounding_box()

    assert any(v > 0 for v in frame.get_image_size())
    assert bbox_is_valid(frame_bbox) and bounding_box_is_intersect_with_image(frame_bbox, frame_image.get_image_size())
    return frame_image, frame_bbox

def _get_frame_image_and_bbox(sequence, index, track_id, rng_engine, reference_bbox):
    frame = sequence.get_frame(index)
    frame_image = frame.get_image_path()
    frame_bbox = _get_bbox_and_generate_dummy_bbox_if_empty(frame, track_id, rng_engine, reference_bbox)
    assert any(v > 0 for v in frame.get_image_size())
    assert bbox_is_valid(frame_bbox) and bounding_box_is_intersect_with_image(frame_bbox, frame.get_image_size())
    return frame_image, frame_bbox

def _data_getter(sequence: MultipleObjectTrackingDatasetSequence_MemoryMapped, track_id, index_of_frames, rng_engine: np.random.Generator):
    z_image, z_bbox = _get_positive_frame_image_and_bbox(sequence, index_of_frames[0], track_id)

    if len(index_of_frames) == 1:
        return ((z_image, z_bbox), )

    x_image, x_bbox = _get_frame_image_and_bbox(sequence, index_of_frames[1], track_id, rng_engine, z_bbox)
    return ((z_image, z_bbox), (x_image, x_bbox))

def _sampling_one_track_in_sequence_and_generate_object_visible_mask(sequence: MultipleObjectTrackingDatasetSequence_MemoryMapped, rng_engine: np.random.Generator):
    index_of_track = rng_engine.integers(0, sequence.get_number_of_objects())
    track = sequence.get_object(index_of_track)

    mask = np.zeros(len(sequence), dtype=np.bool_)
    if track.get_all_bounding_box_validity_flag() is not None:
        ind = track.get_all_frame_index()[track.get_all_bounding_box_validity_flag()]
        mask[ind] = True
    else:
        mask[track.get_all_frame_index()] = True
    return mask, track.get_id()

def do_positive_sampling_in_multiple_object_tracking_dataset_sequence(sequence: MultipleObjectTrackingDatasetSequence_MemoryMapped, frame_range: int, sampling_method: SiamesePairSamplingMethod, rng_engine: np.random.Generator):
    mask, track_id = _sampling_one_track_in_sequence_and_generate_object_visible_mask(sequence, rng_engine)

    return _data_getter(sequence, track_id, do_siamfc_pair_sampling_positive_only(len(sequence), frame_range, mask, sampling_method, rng_engine), rng_engine)

def do_negative_sampling_in_multiple_object_tracking_dataset_sequence(sequence: MultipleObjectTrackingDatasetSequence_MemoryMapped, frame_range: int, rng_engine: np.random.Generator):
    mask, track_id = _sampling_one_track_in_sequence_and_generate_object_visible_mask(sequence, rng_engine)

    return _data_getter(sequence, track_id, do_siamfc_pair_sampling_negative_only(len(sequence), frame_range, mask, rng_engine), rng_engine)

def do_sampling_in_multiple_object_tracking_dataset_sequence(sequence: MultipleObjectTrackingDatasetSequence_MemoryMapped, frame_range: int, sampling_method: SiamesePairSamplingMethod, rng_engine: np.random.Generator):
    mask, track_id = _sampling_one_track_in_sequence_and_generate_object_visible_mask(sequence, rng_engine)

    indices, is_positive = do_siamfc_pair_sampling(len(sequence), frame_range, mask, sampling_method, rng_engine)
    return _data_getter(sequence, track_id, indices, rng_engine), is_positive


text:
import time
import unittest
import logging
from unittest.mock import patch, Mock
import sys
import numpy as np
import pandas as pd
import pprint
from google.protobuf.json_format import MessageToDict

import graphsignal
from graphsignal import statistics
from graphsignal.windows import *
from graphsignal import metrics_pb2

logger = logging.getLogger('graphsignal')

class StatisticsTest(unittest.TestCase):
    def setUp(self):
        self.maxDiff = None
        if len(logger.handlers) == 0:
            logger.addHandler(logging.StreamHandler(sys.stdout))
        graphsignal.configure(api_key='k1', debug_mode=True)

    def tearDown(self):
        graphsignal.shutdown()

    def test_data_size(self):
        self.assertEqual(statistics.estimate_size([[1], [2]]), 2)
        self.assertEqual(statistics.estimate_size(np.array([[1], [2]])), 2)
        self.assertEqual(statistics.estimate_size(
            {'a': [1, 2, 3], 'b': [4, 5, 6]}), 3)
        self.assertEqual(statistics.estimate_size(
            {'a': np.asarray([1, 2, 3]), 'b': np.asarray([4, 5, 6])}), 3)
        self.assertEqual(
            statistics.estimate_size(
                pd.DataFrame(
                    data=[
                        [1],
                        [2]])),
            2)

    def test_convert_to_numpy_dict_of_lists(self):
        d1 = {'c1': [1, 2], 'c2': [3.1, 4.1]}
        d2 = {'c1': [3], 'c2': [5.1]}
        columns, data = statistics._convert_to_numpy([d1, d2])
        self.assertIsNotNone(data)
        self.assertEqual(columns, ['c1', 'c2'])
        self.assertTrue(np.array_equal(
            data[0], np.asarray([1, 2, 3])))
        self.assertTrue(np.array_equal(
            data[1], np.asarray([3.1, 4.1, 5.1])))
        self.assertEqual(data[0].dtype, 'int64')
        self.assertEqual(data[1].dtype, 'float64')

    def test_convert_to_numpy_list_of_list(self):
        d1 = [[1, 2], [3.1, 4.1]]
        d2 = [[3, 4], [5.1, 6.1]]
        columns, data = statistics._convert_to_numpy(
            [d1, d2], columns=['A', 'B'])
        self.assertIsNotNone(data)
        self.assertEqual(columns, ['A', 'B'])
        self.assertTrue(np.array_equal(
            data[0], np.asarray([1, 3.1, 3, 5.1])))
        self.assertTrue(np.array_equal(
            data[1], np.asarray([2, 4.1, 4, 6.1])))

    def test_convert_to_numpy_dataframe_with_columns(self):
        d1 = pd.DataFrame(data={'c1': [1, 2], 'c2': [3.1, 4.1]})
        d2 = pd.DataFrame(data={'c1': [3], 'c2': [5.1]})
        columns, data = statistics._convert_to_numpy([d1, d2])
        self.assertIsNotNone(data)
        self.assertEqual(columns, ['c1', 'c2'])
        self.assertTrue(np.array_equal(
            data[0], np.asarray([1, 2, 3])))
        self.assertTrue(np.array_equal(
            data[1], np.asarray([3.1, 4.1, 5.1])))
        self.assertEqual(data[0].dtype, 'int64')
        self.assertEqual(data[1].dtype, 'float64')

    def test_convert_to_numpy_dataframe_wo_columns(self):
        d = pd.DataFrame(data=[[1, 3.1], [2, 4.1]])
        columns, data = statistics._convert_to_numpy([d])
        self.assertIsNotNone(data)
        self.assertEqual(columns, ['0', '1'])
        self.assertTrue(np.array_equal(
            data[0], np.asarray([1, 2])))
        self.assertTrue(np.array_equal(
            data[1], np.asarray([3.1, 4.1])))

    def test_convert_to_numpy_ndarray(self):
        d1 = np.asarray([[1, 3.1], [2, 4.1]])
        d2 = np.asarray([[8, 12]])
        columns, data = statistics._convert_to_numpy([d1, d2])
        self.assertIsNotNone(data)
        self.assertEqual(columns, ['0', '1'])
        self.assertTrue(np.array_equal(
            data[0], np.asarray([1, 2, 8])))
        self.assertTrue(np.array_equal(
            data[1], np.asarray([3.1, 4.1, 12])))
        self.assertEqual(data[0].dtype, 'float64')
        self.assertEqual(data[1].dtype, 'float64')

    def test_update_data_metrics_empty(self):
        window = metrics_pb2.MetricWindow()
        statistics.update_data_metrics({}, window, [])

        self.assertEqual(len(window.data_streams[str(
            metrics_pb2.DataStream.DataSource.INPUT_FEATURES)].metrics), 0)
        self.assertEqual(len(window.data_streams[str(
            metrics_pb2.DataStream.DataSource.OUTPUTS)].metrics), 0)

    def test_update_data_metrics(self):
        window = metrics_pb2.MetricWindow()
        metric_updaters = {}

        features = {
            'f1': [1, 1, 2, 0],
            'f2': [3.5, 4.5, 5.5, 0],
            'f3': ['a', 'b', 'c', 'c'],
            'f4': [0, float('nan'), float('inf'), 2],
            'f5': [True, True, False, None]}
        outputs = [[0.1], [0.2], [0.1], [0.4]]

        statistics.update_data_metrics(
            metric_updaters, window, [
                PredictionRecord(
                    features=features, outputs=outputs)])

        for metric_updater in metric_updaters.values():
            metric_updater.finalize()

        window_dict = MessageToDict(window)

        input_metrics_json = window_dict['dataStreams'][str(
            metrics_pb2.DataStream.DataSource.INPUT_FEATURES)]['metrics']
        output_metrics_json = window_dict['dataStreams'][str(
            metrics_pb2.DataStream.DataSource.OUTPUTS)]['metrics']

        for metric in input_metrics_json.values():
            if 'distributionValue' in metric:
                del metric['distributionValue']['sketchKll10']

        for metric in output_metrics_json.values():
            if 'distributionValue' in metric:
                del metric['distributionValue']['sketchKll10']

        #pp = pprint.PrettyPrinter()
        # pp.pprint(input_metrics_json)
        # pp.pprint(output_metrics_json)

        self.assertEqual(input_metrics_json,
                         {'02ada935d60a': {'dimensions': {'feature_name': 'f5'},
                                           'distributionValue': {'sketchImpl': 'KLL10'},
                                           'name': 'distribution',
                                           'type': 'DISTRIBUTION'},
                             '076aa29f76cf': {'dimensions': {'feature_name': 'f4'},
                                              'name': 'integer_values',
                                              'ratioValue': {'counter': 2.0, 'total': 4.0},
                                              'type': 'RATIO'},
                             '222d0b1dfa42': {'dimensions': {'feature_name': 'f2'},
                                              'name': 'missing_values',
                                              'ratioValue': {'total': 4.0},
                                              'type': 'RATIO'},
                             '23777f8d6c5b': {'dimensions': {'feature_name': 'f3'},
                                              'name': 'string_values',
                                              'ratioValue': {'counter': 4.0, 'total': 4.0},
                                              'type': 'RATIO'},
                             '2df60f25840a': {'dimensions': {'feature_name': 'f3'},
                                              'distributionValue': {'sketchImpl': 'KLL10'},
                                              'name': 'distribution',
                                              'type': 'DISTRIBUTION'},
                             '31cdc4047da9': {'dimensions': {'feature_name': 'f5'},
                                              'name': 'missing_values',
                                              'ratioValue': {'counter': 1.0, 'total': 4.0},
                                              'type': 'RATIO'},
                             '342246abc789': {'dimensions': {'feature_name': 'f1'},
                                              'name': 'missing_values',
                                              'ratioValue': {'total': 4.0},
                                              'type': 'RATIO'},
                             '38b5e7971c3f': {'dimensions': {'feature_name': 'f1'},
                                              'distributionValue': {'sketchImpl': 'KLL10'},
                                              'name': 'distribution',
                                              'type': 'DISTRIBUTION'},
                             '3aa8ce8a745e': {'counterValue': {'counter': 4.0},
                                              'name': 'instance_count',
                                              'type': 'COUNTER'},
                             '4283d6f1fd2b': {'dimensions': {'feature_name': 'f5'},
                                              'name': 'boolean_values',
                                              'ratioValue': {'counter': 3.0, 'total': 4.0},
                                              'type': 'RATIO'},
                             '4a9116752484': {'dimensions': {'feature_name': 'f4'},
                                              'name': 'float_values',
                                              'ratioValue': {'total': 4.0},
                                              'type': 'RATIO'},
                             '5f5ec237ea4f': {'dimensions': {'feature_name': 'f1'},
                                              'name': 'integer_values',
                                              'ratioValue': {'counter': 4.0, 'total': 4.0},
                                              'type': 'RATIO'},
                             '62480aeb7da2': {'dimensions': {'feature_name': 'f1'},
                                              'name': 'float_values',
                                              'ratioValue': {'total': 4.0},
                                              'type': 'RATIO'},
                             '6add4609813c': {'gaugeValue': {'gauge': 5.0},
                                              'name': 'column_count',
                                              'type': 'GAUGE'},
                             '709339fc41ee': {'dimensions': {'feature_name': 'f1'},
                                              'name': 'zero_values',
                                              'ratioValue': {'counter': 1.0, 'total': 4.0},
                                              'type': 'RATIO'},
                             '72dc8927bd35': {'dimensions': {'feature_name': 'f2'},
                                              'name': 'zero_values',
                                              'ratioValue': {'counter': 1.0, 'total': 4.0},
                                              'type': 'RATIO'},
                             '8bcee9cf32e3': {'dimensions': {'feature_name': 'f4'},
                                              'distributionValue': {'sketchImpl': 'KLL10'},
                                              'name': 'distribution',
                                              'type': 'DISTRIBUTION'},
                             'c236202ee147': {'dimensions': {'feature_name': 'f4'},
                                              'name': 'zero_values',
                                              'ratioValue': {'counter': 1.0, 'total': 4.0},
                                              'type': 'RATIO'},
                             'd00eadf9e4cf': {'dimensions': {'feature_name': 'f2'},
                                              'name': 'integer_values',
                                              'ratioValue': {'counter': 1.0, 'total': 4.0},
                                              'type': 'RATIO'},
                             'd082b0be0bf4': {'dimensions': {'feature_name': 'f3'},
                                              'name': 'missing_values',
                                              'ratioValue': {'total': 4.0},
                                              'type': 'RATIO'},
                             'df1ea8f7032a': {'dimensions': {'feature_name': 'f4'},
                                              'name': 'missing_values',
                                              'ratioValue': {'counter': 2.0, 'total': 4.0},
                                              'type': 'RATIO'},
                             'e2844d940cc2': {'dimensions': {'feature_name': 'f2'},
                                              'name': 'float_values',
                                              'ratioValue': {'counter': 3.0, 'total': 4.0},
                                              'type': 'RATIO'},
                             'e675c4eab4f7': {'dimensions': {'feature_name': 'f3'},
                                              'name': 'empty_values',
                                              'ratioValue': {'total': 4.0},
                                              'type': 'RATIO'},
                             'f1d8c84f680a': {'dimensions': {'feature_name': 'f2'},
                                              'distributionValue': {'sketchImpl': 'KLL10'},
                                              'name': 'distribution',
                                              'type': 'DISTRIBUTION'}})
        self.assertEqual(output_metrics_json,
                         {'094ce9a486c0': {'gaugeValue': {'gauge': 1.0},
                                           'name': 'column_count',
                                           'type': 'GAUGE'},
                          '15f3585b5a2b': {'dimensions': {'output_name': '0'},
                                           'name': 'missing_values',
                                           'ratioValue': {'total': 4.0},
                                           'type': 'RATIO'},
                          '2d4c320b727b': {'dimensions': {'output_name': '0'},
                                           'name': 'integer_values',
                                           'ratioValue': {'total': 4.0},
                                           'type': 'RATIO'},
                          '350dc33d58cf': {'dimensions': {'output_name': '0'},
                                           'name': 'zero_values',
                                           'ratioValue': {'total': 4.0},
                                           'type': 'RATIO'},
                          '3dea0616b753': {'dimensions': {'output_name': '0'},
                                           'name': 'float_values',
                                           'ratioValue': {'counter': 4.0, 'total': 4.0},
                                           'type': 'RATIO'},
                          '59abb90a51e4': {'dimensions': {'output_name': '0'},
                                           'distributionValue': {'sketchImpl': 'KLL10'},
                                           'name': 'distribution',
                                           'type': 'DISTRIBUTION'},
                          '8e1bc8013472': {'counterValue': {'counter': 4.0},
                                           'name': 'instance_count',
                                           'type': 'COUNTER'}})

    def test_update_data_metrics_perf(self):
        window = metrics_pb2.MetricWindow()
        metric_updaters = {}

        prediction_buffer = []
        for i in range(100):
            features = {}
            for j in range(50):
                features['f1_' + str(j)] = [i]
                features['f2_' + str(j)] = [i + 3.53242342]
                features['f3_' + str(j)] = [100 * 'abc' + str(i)]
                features['f4_' + str(j)] = [float('nan')]
                features['f4_' + str(j)] = [float('+inf')]
            outputs = [[0.1], [0.2], [0.1], [0.4]]
            prediction_buffer.append(
                PredictionRecord(
                    features=features,
                    outputs=outputs))

        #import cProfile
        #from pstats import Stats, SortKey

        start_ts = time.time()

        # with cProfile.Profile() as pr:
        for k in range(10):
            statistics.update_data_metrics(
                metric_updaters, window, prediction_buffer)
        #stats = Stats(pr)
        # stats.sort_stats(SortKey.CUMULATIVE).print_stats(25)

        for metric_updater in metric_updaters.values():
            metric_updater.finalize()

        took = time.time() - start_ts

        print('update_data_metrics took (sec): ', round(took, 6))
        self.assertTrue(len(window.data_streams[str(
            metrics_pb2.DataStream.DataSource.INPUT_FEATURES)].metrics) > 0)
        self.assertTrue(len(window.data_streams[str(
            metrics_pb2.DataStream.DataSource.OUTPUTS)].metrics) > 0)
        self.assertTrue(took < 1)

    def test_truncate_strings(self):
        truncated = statistics._truncate_strings(
            ['abc', '1234567890abcdefgh', '1234567890abcdefghi'])
        self.assertEqual(
            truncated, [
                'abc', '1234567890abcdefgh', '1234567890...efghi'])

    def test_update_performance_metrics_binary(self):
        window = metrics_pb2.MetricWindow()
        metric_updaters = {}

        evaluation_buffer = [
            EvaluationRecord(
                label=True,
                prediction=True,
                segments=[
                    's1',
                    's2']),
            EvaluationRecord(
                label=True,
                prediction=True,
                segments=[
                    's1',
                    's3']),
            EvaluationRecord(
                label=True,
                prediction=False,
                segments=[
                    's1',
                    's3']),
            EvaluationRecord(
                label=False,
                prediction=True,
                segments=[
                    's1',
                    's3']),
            EvaluationRecord(
                label=False,
                prediction=False,
                segments=[
                    's1',
                    's2'])
        ]

        statistics.update_performance_metrics(
            metric_updaters, window, evaluation_buffer)

        for metric_updater in metric_updaters.values():
            metric_updater.finalize()

        window_dict = MessageToDict(window)
        metrics_json = window_dict['dataStreams'][str(
            metrics_pb2.DataStream.DataSource.EVALUATION_BINARY)]['metrics']

        #pp = pprint.PrettyPrinter()
        # pp.pprint(metrics_json)

        self.assertEqual(metrics_json, {'2a8a320fc710': {'name': 'accuracy',
                                                         'ratioValue': {'counter': 3.0, 'total': 5.0},
                                                         'type': 'RATIO'},
                                        '2a915536a163': {'counterValue': {'counter': 1.0},
                                                         'name': 'binary_false_positives',
                                                         'type': 'COUNTER'},
                                        '3e98594499f3': {'counterValue': {'counter': 1.0},
                                                         'name': 'binary_false_negatives',
                                                         'type': 'COUNTER'},
                                        '4769fae29cea': {'counterValue': {'counter': 1.0},
                                                         'name': 'binary_true_negatives',
                                                         'type': 'COUNTER'},
                                        'a335baad03d4': {'distributionValue': {'sketchImpl': 'KLL10',
                                                                               'sketchKll10': {'H': '1',
                                                                                'c': 0.6666666666666666,
                                                                                'compactorsString': [{'items': ['s1',
                                                                                's2',
                                                                                's1',
                                                                                's3',
                                                                                's1',
                                                                                's2']}],
                                                                                'itemType': 'STRING',
                                                                                'k': '10',
                                                                                'maxSize': '11',
                                                                                'size': '6'}},
                                                         'name': 'segment_matches',
                                                         'type': 'DISTRIBUTION'},
                                        'e98d593706c6': {'distributionValue': {'sketchImpl': 'KLL10',
                                                                               'sketchKll10': {'H': '1',
                                                                                'c': 0.6666666666666666,
                                                                                'compactorsString': [{'items': ['s1',
                                                                                's2',
                                                                                's1',
                                                                                's3',
                                                                                's1',
                                                                                's3',
                                                                                's1',
                                                                                's3',
                                                                                's1',
                                                                                's2']}],
                                                                                'itemType': 'STRING',
                                                                                'k': '10',
                                                                                'maxSize': '11',
                                                                                'size': '10'}},
                                                         'name': 'segment_totals',
                                                         'type': 'DISTRIBUTION'},
                                        'ee427ea584e7': {'counterValue': {'counter': 2.0},
                                                         'name': 'binary_true_positives',
                                                         'type': 'COUNTER'}})

    def test_update_performance_metrics_categorical(self):
        window = metrics_pb2.MetricWindow()
        metric_updaters = {}

        evaluation_buffer = [
            EvaluationRecord(
                label='c1',
                prediction='c1',
                segments=[
                    's1',
                    's2']),
            EvaluationRecord(
                label='c1',
                prediction='c1',
                segments=[
                    's1',
                    's3']),
            EvaluationRecord(
                label='c1',
                prediction='c2',
                segments=[
                    's1',
                    's3']),
            EvaluationRecord(
                label='c1',
                prediction='c2',
                segments=[
                    's1',
                    's3']),
            EvaluationRecord(
                label='c2',
                prediction='c1',
                segments=[
                    's1',
                    's3']),
            EvaluationRecord(
                label='c2',
                prediction='c2',
                segments=[
                    's1',
                    's2']),
            EvaluationRecord(
                label='c3',
                prediction='c3',
                segments=[
                    's1',
                    's2']),
        ]

        statistics.update_performance_metrics(
            metric_updaters, window, evaluation_buffer)

        for metric_updater in metric_updaters.values():
            metric_updater.finalize()

        window_dict = MessageToDict(window)
        metrics_json = window_dict['dataStreams'][str(
            metrics_pb2.DataStream.DataSource.EVALUATION_CATEGORICAL)]['metrics']

        data_stream = window.data_streams[str(
            metrics_pb2.DataStream.DataSource.EVALUATION_CATEGORICAL)]

        total = None
        sketches = {}
        for metric in data_stream.metrics.values():
            if metric.name == 'total':
                total = metric.counter_value.counter
            elif metric.name == 'class_true_positives':
                sketches['class_true_positives'] = KLLSketch()
                sketches['class_true_positives'].from_proto(
                    metric.distribution_value.sketch_kll10)
            elif metric.name == 'class_false_positives':
                sketches['class_false_positives'] = KLLSketch()
                sketches['class_false_positives'].from_proto(
                    metric.distribution_value.sketch_kll10)
            elif metric.name == 'class_false_negatives':
                sketches['class_false_negatives'] = KLLSketch()
                sketches['class_false_negatives'].from_proto(
                    metric.distribution_value.sketch_kll10)

        # c1 2f22765d
        # c2 6b1f5330
        # c3 a625406f
        self.assertEqual(total, 7.0)
        self.assertEqual(sketches['class_true_positives'].distribution(),
                         [['2f22765d', 2], ['6b1f5330', 1], ['a625406f', 1]])
        self.assertEqual(sketches['class_false_positives'].distribution(),
                         [['2f22765d', 1], ['6b1f5330', 2]])
        self.assertEqual(sketches['class_false_negatives'].distribution(),
                         [['2f22765d', 2], ['6b1f5330', 1]])

    def test_update_performance_metrics_numeric(self):
        window = metrics_pb2.MetricWindow()
        metric_updaters = {}

        evaluation_buffer = [
            EvaluationRecord(
                label=1.0,
                prediction=3.4,
                segments=[
                    's1',
                    's2']),
            EvaluationRecord(
                label=2.4,
                prediction=0.2,
                segments=[
                    's1',
                    's3'])
        ]

        statistics.update_performance_metrics(
            metric_updaters, window, evaluation_buffer)

        for metric_updater in metric_updaters.values():
            metric_updater.finalize()

        window_dict = MessageToDict(window)
        metrics_json = window_dict['dataStreams'][str(
            metrics_pb2.DataStream.DataSource.EVALUATION_NUMERIC)]['metrics']

        #pp = pprint.PrettyPrinter()
        # pp.pprint(metrics_json)

        self.assertEqual(metrics_json, {'2789175b4787': {'counterValue': {'counter': 4.6},
                                                         'name': 'mae_sum',
                                                         'type': 'COUNTER'},
                                        '8c82f7b3ed54': {'counterValue': {'counter': 2.0},
                                                         'name': 'mse_n',
                                                         'type': 'COUNTER'},
                                        'fe4a651921fe': {'counterValue': {'counter': 10.599999999999998},
                                                         'name': 'mse_sum',
                                                         'type': 'COUNTER'}})

    def test_update_performance_metrics_segments(self):
        window = metrics_pb2.MetricWindow()
        metric_updaters = {}

        evaluation_buffer = [
            EvaluationRecord(
                label=True,
                prediction=True,
                segments=[
                    's1',
                    's2']),
            EvaluationRecord(
                label=True,
                prediction=True,
                segments=[
                    's1',
                    's3']),
            EvaluationRecord(
                label=True,
                prediction=False,
                segments=[
                    's1',
                    's3']),
            EvaluationRecord(
                label=False,
                prediction=True,
                segments=[
                    's1',
                    's3']),
            EvaluationRecord(
                label=False,
                prediction=False,
                segments=[
                    's1',
                    's2'])
        ]

        statistics.update_performance_metrics(
            metric_updaters, window, evaluation_buffer)

        for metric_updater in metric_updaters.values():
            metric_updater.finalize()

        data_stream = window.data_streams[str(
            metrics_pb2.DataStream.DataSource.EVALUATION_BINARY)]

        sketches = {}
        for metric in data_stream.metrics.values():
            if metric.name == 'segment_totals':
                sketches['segment_totals'] = KLLSketch()
                sketches['segment_totals'].from_proto(
                    metric.distribution_value.sketch_kll10)
            elif metric.name == 'segment_matches':
                sketches['segment_matches'] = KLLSketch()
                sketches['segment_matches'].from_proto(
                    metric.distribution_value.sketch_kll10)

        self.assertEqual(sketches['segment_totals'].distribution(),
                         [['s1', 5], ['s2', 2], ['s3', 3]])
        self.assertEqual(sketches['segment_matches'].distribution(), [
                         ['s1', 3], ['s2', 2], ['s3', 1]])

max_stars_repo_name: graphsignal/graphsignal

max_issues_repo_name: graphsignal/graphsignal-python

max_forks_repo_name: graphsignal/graphsignal-python


text:
import os
import sys
curDir = os.path.dirname(__file__)
sys.path.append('{0}/../scripts/'.format(curDir))
sys.path.append('{0}/../models/'.format(curDir))

import pandas as pd
import datetime
import numpy as np
from indicators import Indicators
# from auto_sklearn_model import AutoSklearnModel
from svm_model import SVMModel

START_DATE = '2007-01-01'
END_DATE = datetime.date.today()

def create_classify_labels(adj_close_prices, num_days):
    classified = adj_close_prices.rolling(window=num_days+1).\
                    apply(lambda t: 1 if t[num_days] >= t[0] else -1)[num_days:]
    return pd.DataFrame(classified.values, columns=['Labels'])

def get_indicators(stock):
    start = pd.to_datetime(START_DATE)
    end = pd.to_datetime(END_DATE)
    indicators = Indicators(stock, start, end)
    return indicators.calculate_all_indicators()

def buy_or_sell(stock, num_days):
    # Get indicators data and labels
    indicators_df_origin = get_indicators(stock)
    train_Y_df = create_classify_labels(indicators_df_origin['Adj Close Price'], num_days)
    train_Y = np.array(train_Y_df.transpose().values[0])
    train_X_df = indicators_df_origin[:len(train_Y_df)]
    test_X_df = indicators_df_origin[len(train_Y_df):]

    # AutoSklearnModel
    # auto_model = AutoSklearnModel()
    # (model, training_score) = auto_model.get_model(train_X, train_Y)
    # predicted = model.predict(test_X)

    # SVM Classification
    svm_model = SVMModel()
    training_score = svm_model.train(train_X_df, train_Y)
    predicted = svm_model.predict(test_X_df)

    dates = pd.to_datetime(test_X_df.index)

    for i in range(len(predicted)):
        print("    {0}".format(dates[i].strftime('%Y-%m-%d')))
        print("    next {0}-day prediction: {1}".format(num_days, predicted[i]))
        print("")

    print("Training data size was: {0}".format(len(train_X_df)))
    print("Training accuracy was: {0:.4f}".format(training_score))
    print("")

    return predicted[0]

if __name__ == "__main__":

    if len(sys.argv) != 3:
        print("Usage: python3 predict.py <num_days> <stock_symbol>")
        quit()

    buy_or_sell(sys.argv[2], int(sys.argv[1]))


text:
########################################################
# divide_spectra.py                                    #
# Matheus J. Castro                                    #
# Version 2.0                                          #
# Last Modification: 11/11/2021 (month/day/year)       #
# https://github.com/MatheusJCastro/spectra_comparator #
# Licensed under MIT License                           #
########################################################

from astropy.io import fits
import matplotlib.pyplot as plt
import numpy as np
import sys

def open_spec(fl_name):
    # Subroutine to open the .fits spectrum and read it

    hdul = fits.open(fl_name)  # open the file
    spec_data = hdul[0].data  # get the data
    spec_header = hdul[0].header  # get the header

    if spec_data.shape != (2048,):  # get only the actual spectrum (for multidimensional data)
        spec_data = spec_data[1][0]

    # Get the wavelength information from the header
    # CDELT1 or CD1_1
    wl = spec_header['CRVAL1'] + spec_header['CD1_1'] * np.arange(0, len(spec_data))

    hdul.close()  # close the file

    return wl, spec_data, spec_header

def plot_spectra(spec1, spec2, specdiv, show=False, sharey=False):
    # Subroutine to plot the spectra

    specs = [spec1, spec2, specdiv]  # spectra data
    names = ["New Lamp", "Old Lamp", "Ratio"]  # spectra names
    colors = ["blue", "orange", "green"]  # spectra colors

    # noinspection PyTypeChecker
    fig, axes = plt.subplots(figsize=(21, 9), nrows=len(specs), sharex=True, sharey=sharey)
    axes[0].set_title("Difference between lamps", fontsize=28)

    for i, spec in enumerate(specs):
        axes[i].plot(spec[0], spec[1], label=names[i], color=colors[i])
        # axes[i].set_yscale("log")
        # axes[i].set_yticklabels([])
        axes[i].tick_params(axis='y', which='both', labelsize=16)
        axes[i].yaxis.set_major_formatter(plt.FormatStrFormatter('%.2f'))
        axes[i].legend(loc="upper left", fontsize=20)
        axes[i].grid(which="both")
        if i == 2 and not sharey:
            axes[i].set_ylim(0.99*min(spec[1]), 0.43*max(spec[1]))

    axes[len(specs)-1].set_xlabel("Wavelength (\u212b)", fontsize=20)
    plt.xticks(fontsize=16)
    plt.tight_layout()
    plt.subplots_adjust(hspace=.0)  # no horizontal space between the plots

    shy = ""
    if sharey:
        shy = "_sharey"

    plt.savefig("Lamps_Comparison{}".format(shy))

    if show:
        plt.show()

    plt.close()

def main(arg):
    # Main subroutine to plot the spectra

    # If one of the arguments passed is "sharey" the y-axis will be shared among the subplots
    sharey = False
    if any(i == "sharey" for i in arg):
        sharey = True

    # File names. 3 must be given
    files = ["multi_tha_novo.fits",
             "tha_velha.fits",
             "tha_divided.fits"]

    # Open the spectra
    old_spec = open_spec(files[0])
    new_spec = open_spec(files[1])
    div_spec = open_spec(files[2])

    for i in [True, False]:  # Bypass of the sharey variable
        plot_spectra(old_spec, new_spec, div_spec, sharey=i)

if __name__ == '__main__':
    args = sys.argv[1:]  # get the arguments passed in the execution
    main(args)

reason: import numpy,from astropy