text:
"""

===============
Incremental PCA
===============

Incremental principal component analysis (IPCA) is typically used as a
replacement for principal component analysis (PCA) when the dataset to be
decomposed is too large to fit in memory. IPCA builds a low-rank approximation
for the input data using an amount of memory which is independent of the
number of input data samples. It is still dependent on the input data features,
but changing the batch size allows for control of memory usage.

This example serves as a visual check that IPCA is able to find a similar
projection of the data to PCA (to a sign flip), while only processing a
few samples at a time. This can be considered a "toy example", as IPCA is
intended for large datasets which do not fit in main memory, requiring
incremental approaches.

"""
print(__doc__)

# Authors: Kyle Kastner
# License: BSD 3 clause

import numpy as np
import matplotlib.pyplot as plt

from sklearn.datasets import load_iris
from sklearn.decomposition import PCA, IncrementalPCA

iris = load_iris()
X = iris.data
y = iris.target

n_components = 2
ipca = IncrementalPCA(n_components=n_components, batch_size=10)
X_ipca = ipca.fit_transform(X)

pca = PCA(n_components=n_components)
X_pca = pca.fit_transform(X)

colors = ['navy', 'turquoise', 'darkorange']

for X_transformed, title in [(X_ipca, "Incremental PCA"), (X_pca, "PCA")]:
    plt.figure(figsize=(8, 8))
    for color, i, target_name in zip(colors, [0, 1, 2], iris.target_names):
        plt.scatter(X_transformed[y == i, 0], X_transformed[y == i, 1],
                    color=color, lw=2, label=target_name)

    if "Incremental" in title:
        err = np.abs(np.abs(X_pca) - np.abs(X_ipca)).mean()
        plt.title(title + " of iris dataset\nMean absolute unsigned error "
                  "%.6f" % err)
    else:
        plt.title(title + " of iris dataset")
    plt.legend(loc="best", shadow=False, scatterpoints=1)
    plt.axis([-4, 4, -1.5, 1.5])

#plt.show()

ext: py

reason: import numpy


text:
import warnings

import mmcv
import numpy as np

from mmpose.core import imshow_keypoints, imshow_keypoints_3d
from .. import builder
from ..builder import POSENETS
from .base import BasePose

try:
    from mmcv.runner import auto_fp16
except ImportError:
    warnings.warn('auto_fp16 from mmpose will be deprecated from v0.15.0'
                  'Please install mmcv>=1.1.4')
    from mmpose.core import auto_fp16

@POSENETS.register_module()
class PoseLifter(BasePose):
    """Pose lifter that lifts 2D pose to 3D pose.

    The basic model is a pose model that predicts root-relative pose. If
    traj_head is not None, a trajectory model that predicts absolute root joint
    position is also built.

    Args:
        backbone (dict): Config for the backbone of pose model.
        neck (dict|None): Config for the neck of pose model.
        keypoint_head (dict|None): Config for the head of pose model.
        traj_backbone (dict|None): Config for the backbone of trajectory model.
            If traj_backbone is None and traj_head is not None, trajectory
            model will share backbone with pose model.
        traj_neck (dict|None): Config for the neck of trajectory model.
        traj_head (dict|None): Config for the head of trajectory model.
        loss_semi (dict|None): Config for semi-supervision loss.
        train_cfg (dict|None): Config for keypoint head during training.
        test_cfg (dict|None): Config for keypoint head during testing.
        pretrained (str|None): Path to pretrained weights.
    """

    def __init__(self,
                 backbone,
                 neck=None,
                 keypoint_head=None,
                 traj_backbone=None,
                 traj_neck=None,
                 traj_head=None,
                 loss_semi=None,
                 train_cfg=None,
                 test_cfg=None,
                 pretrained=None):
        super().__init__()
        self.fp16_enabled = False

        self.train_cfg = train_cfg
        self.test_cfg = test_cfg

        # pose model
        self.backbone = builder.build_backbone(backbone)

        if neck is not None:
            self.neck = builder.build_neck(neck)

        if keypoint_head is not None:
            keypoint_head['train_cfg'] = train_cfg
            keypoint_head['test_cfg'] = test_cfg
            self.keypoint_head = builder.build_head(keypoint_head)

        # trajectory model
        if traj_head is not None:
            self.traj_head = builder.build_head(traj_head)

            if traj_backbone is not None:
                self.traj_backbone = builder.build_backbone(traj_backbone)
            else:
                self.traj_backbone = self.backbone

            if traj_neck is not None:
                self.traj_neck = builder.build_neck(traj_neck)

        # semi-supervised learning
        self.semi = loss_semi is not None
        if self.semi:
            assert keypoint_head is not None and traj_head is not None
            self.loss_semi = builder.build_loss(loss_semi)

        self.init_weights(pretrained=pretrained)

    @property
    def with_neck(self):
        """Check if has keypoint_neck."""
        return hasattr(self, 'neck')

    @property
    def with_keypoint(self):
        """Check if has keypoint_head."""
        return hasattr(self, 'keypoint_head')

    @property
    def with_traj_backbone(self):
        """Check if has trajectory_backbone."""
        return hasattr(self, 'traj_backbone')

    @property
    def with_traj_neck(self):
        """Check if has trajectory_neck."""
        return hasattr(self, 'traj_neck')

    @property
    def with_traj(self):
        """Check if has trajectory_head."""
        return hasattr(self, 'traj_head')

    def init_weights(self, pretrained=None):
        """Weight initialization for model."""
        self.backbone.init_weights(pretrained)
        if self.with_neck:
            self.neck.init_weights()
        if self.with_keypoint:
            self.keypoint_head.init_weights()
        if self.with_traj_backbone:
            self.traj_backbone.init_weights(pretrained)
        if self.with_traj_neck:
            self.traj_neck.init_weights()
        if self.with_traj:
            self.traj_head.init_weights()

    @auto_fp16(apply_to=('input', ))
    def forward(self,
                input,
                target=None,
                target_weight=None,
                metas=None,
                return_loss=True,
                **kwargs):
        """Calls either forward_train or forward_test depending on whether
        return_loss=True.

        Note:
            Note:
            batch_size: N
            num_input_keypoints: Ki
            input_keypoint_dim: Ci
            input_sequence_len: Ti
            num_output_keypoints: Ko
            output_keypoint_dim: Co
            input_sequence_len: To

        Args:
            input (torch.Tensor[NxKixCixTi]): Input keypoint coordinates.
            target (torch.Tensor[NxKoxCoxTo]): Output keypoint coordinates.
                Defaults to None.
            target_weight (torch.Tensor[NxKox1]): Weights across different
                joint types. Defaults to None.
            metas (list(dict)): Information about data augmentation
            return_loss (bool): Option to `return loss`. `return loss=True`
                for training, `return loss=False` for validation & test.
        Returns:
            dict|Tensor: if `reutrn_loss` is true, return losses. Otherwise
                return predicted poses
        """
        if return_loss:
            return self.forward_train(input, target, target_weight, metas,
                                      **kwargs)
        else:
            return self.forward_test(input, metas, **kwargs)

    def forward_train(self, input, target, target_weight, metas, **kwargs):
        """Defines the computation performed at every call when training."""
        assert input.size(0) == len(metas)

        # supervised learning
        # pose model
        features = self.backbone(input)
        if self.with_neck:
            features = self.neck(features)
        if self.with_keypoint:
            output = self.keypoint_head(features)

        losses = dict()
        if self.with_keypoint:
            keypoint_losses = self.keypoint_head.get_loss(
                output, target, target_weight)
            keypoint_accuracy = self.keypoint_head.get_accuracy(
                output, target, target_weight, metas)
            losses.update(keypoint_losses)
            losses.update(keypoint_accuracy)

        # trajectory model
        if self.with_traj:
            traj_features = self.traj_backbone(input)
            if self.with_traj_neck:
                traj_features = self.traj_neck(traj_features)
            traj_output = self.traj_head(traj_features)

            traj_losses = self.traj_head.get_loss(traj_output,
                                                  kwargs['traj_target'], None)
            losses.update(traj_losses)

        # semi-supervised learning
        if self.semi:
            ul_input = kwargs['unlabeled_input']
            ul_features = self.backbone(ul_input)
            if self.with_neck:
                ul_features = self.neck(ul_features)
            ul_output = self.keypoint_head(ul_features)

            ul_traj_features = self.traj_backbone(ul_input)
            if self.with_traj_neck:
                ul_traj_features = self.traj_neck(ul_traj_features)
            ul_traj_output = self.traj_head(ul_traj_features)

            output_semi = dict(
                labeled_pose=output,
                unlabeled_pose=ul_output,
                unlabeled_traj=ul_traj_output)
            target_semi = dict(
                unlabeled_target_2d=kwargs['unlabeled_target_2d'],
                intrinsics=kwargs['intrinsics'])

            semi_losses = self.loss_semi(output_semi, target_semi)
            losses.update(semi_losses)

        return losses

    def forward_test(self, input, metas, **kwargs):
        """Defines the computation performed at every call when training."""
        assert input.size(0) == len(metas)

        results = {}

        features = self.backbone(input)
        if self.with_neck:
            features = self.neck(features)
        if self.with_keypoint:
            output = self.keypoint_head.inference_model(features)
            keypoint_result = self.keypoint_head.decode(metas, output)
            results.update(keypoint_result)

        if self.with_traj:
            traj_features = self.traj_backbone(input)
            if self.with_traj_neck:
                traj_features = self.traj_neck(traj_features)
            traj_output = self.traj_head.inference_model(traj_features)
            results['traj_preds'] = traj_output

        return results

    def forward_dummy(self, input):
        """Used for computing network FLOPs.

        See ``tools/get_flops.py``.

        Args:
            input (torch.Tensor): Input pose

        Returns:
            Tensor: Model output
        """
        output = self.backbone(input)
        if self.with_neck:
            output = self.neck(output)
        if self.with_keypoint:
            output = self.keypoint_head(output)

        if self.with_traj:
            traj_features = self.traj_backbone(input)
            if self.with_neck:
                traj_features = self.traj_neck(traj_features)
            traj_output = self.traj_head(traj_features)
            output = output + traj_output

        return output

    def show_result(self,
                    result,
                    img=None,
                    skeleton=None,
                    pose_kpt_color=None,
                    pose_limb_color=None,
                    radius=8,
                    thickness=2,
                    vis_height=400,
                    win_name='',
                    show=False,
                    wait_time=0,
                    out_file=None):
        """Visualize 3D pose estimation results.

        Args:
            result (list[dict]): The pose estimation results containing:
                - "keypoints_3d" ([K,4]): 3D keypoints
                - "keypoints" ([K,3] or [T,K,3]): Optional for visualizing
                    2D inputs. If a sequence is given, only the last frame
                    will be used for visualization
                - "bbox" ([4,] or [T,4]): Optional for visualizing 2D inputs
                - "title" (str): title for the subplot
            img (str or Tensor): Optional. The image to visualize 2D inputs on.
            skeleton (list of [idx_i,idx_j]): Skeleton described by a list of
                limbs, each is a pair of joint indices.
            pose_kpt_color (np.array[Nx3]`): Color of N keypoints.
                If None, do not draw keypoints.
            pose_limb_color (np.array[Mx3]): Color of M limbs.
                If None, do not draw limbs.
            radius (int): Radius of circles.
            thickness (int): Thickness of lines.
            vis_height (int): The image hight of the visualization. The width
                will be N*vis_height depending on the number of visualized
                items.
            win_name (str): The window name.
            wait_time (int): Value of waitKey param.
                Default: 0.
            out_file (str or None): The filename to write the image.
                Default: None.

        Returns:
            Tensor: Visualized img, only if not `show` or `out_file`.
        """

        assert len(result) > 0
        result = sorted(result, key=lambda x: x.get('track_id', 0))

        # draw image and input 2d poses
        if img is not None:
            img = mmcv.imread(img)

            bbox_result = []
            pose_input_2d = []
            for res in result:
                if 'bbox' in res:
                    bbox = np.array(res['bbox'])
                    if bbox.ndim != 1:
                        assert bbox.ndim == 2
                        bbox = bbox[-1]  # Get bbox from the last frame
                    bbox_result.append(bbox)
                if 'keypoints' in res:
                    kpts = np.array(res['keypoints'])
                    if kpts.ndim != 2:
                        assert kpts.ndim == 3
                        kpts = kpts[-1]  # Get 2D keypoints from the last frame
                    pose_input_2d.append(kpts)

            if len(bbox_result) > 0:
                bboxes = np.vstack(bbox_result)
                mmcv.imshow_bboxes(
                    img,
                    bboxes,
                    colors='green',
                    top_k=-1,
                    thickness=thickness,
                    show=False)
            if len(pose_input_2d) > 0:
                imshow_keypoints(
                    img,
                    pose_input_2d,
                    skeleton,
                    kpt_score_thr=0.3,
                    pose_kpt_color=pose_kpt_color,
                    pose_limb_color=pose_limb_color,
                    radius=radius,
                    thickness=thickness)
            img = mmcv.imrescale(img, scale=vis_height / img.shape[0])

        img_vis = imshow_keypoints_3d(result, img, skeleton, pose_kpt_color,
                                      pose_limb_color, vis_height)

        if show:
            mmcv.visualization.imshow(img_vis, win_name, wait_time)

        if out_file is not None:
            mmcv.imwrite(img_vis, out_file)

        return img_vis

max_stars_repo_name: SummerVideoAnalysis/mmpose

max_issues_repo_name: connor-john/mmpose

max_forks_repo_name: connor-john/mmpose


text:
import numpy as np

class DataPreprocessorModel(object):
    def __init__(self, X, labels):
        """

        训练数据集，将数据分为训练集、验证集、测试集，比例为2：1：1
        并将数据集预处理

        Parameters
        ----------
        X : {array-like, sparse matrix} of shape(n_samples, n_features)
            Training data

        labels : {array-like} of shape(n_samples, 1_label)

        """
        self.X = X
        self.labels = labels

    def __normalize(self):
        # (X-U)/delta
        mean = np.mean(self.X, axis = 0)
        X_norm = self.X - mean
        std = np.std(X_norm, axis = 0)
        X_norm /= std + 10 ** (-5)
        return (X_norm, mean, std)

    def __PCA_white(self, X_train):
        mean  = np.mean(X_train, axis = 0)
        X_norm = X_train - mean
        cov = np.dot(X_norm.T, X_norm)/X_norm.shape[0]
        U, S, V = np.linalg.svd(cov)
        X_norm = np.dot(X_norm, U)
        X_norm /= np.sqrt(S + 10**(-5))
        return (X_norm, mean, U, S)

    def __split_data(self):

        num_examples = self.X.shape[0]
        shuffle_no = list(range(num_examples))
        np.random.shuffle(shuffle_no)

        X_train = self.X[ shuffle_no[:num_examples//2] ]
        labels_train = labels[ shuffle_no[:num_examples//2] ]

        X_val = self.X[ shuffle_no[num_examples//2:num_examples//2+num_examples//4]]
        labels_val = labels[ shuffle_no[num_examples//2:num_examples//2+num_examples//4]]

        X_test = self.X[shuffle_no[-num_examples//4 :]]
        labels_test = labels[ shuffle_no[-num_examples//4 :]]

        return (X_train, labels_train, X_val, labels_val, X_test, labels_test)

    def __data_preprocess(self, X_train, X_val, X_test):
        (X_train_pca, mean, U, S) = self.__PCA_white(X_train)
        X_val_pca = np.dot(X_val - mean, U)
        X_val_pca /= np.sqrt(S + 10**(-5))
        X_test_pca = np.dot(X_test - mean, U)
        X_test_pca /= np.sqrt(S + 10**(-5))
        return (X_train_pca, X_val_pca, X_test_pca)

    def fit(self):
        """

        训练数据集，将数据分为训练集、验证集、测试集，比例为2：1：1
        并将数据集预处理(中心化,规范化,PCA处理)

        Returns
        -------
        X_train_pca : 训练集X

        labels_train : 训练集标记集合

        X_val_pca : 验证集X

        labels_val : 验证集标记集合

        X_test_pca : 测试集X

        labels_test : 测试集标记集合

        """
        X_train, labels_train, X_val, labels_val, X_test, labels_test = self.__split_data()
        X_train_pca, X_val_pca, X_test_pca = self.__data_preprocess(X_train=X_train, X_val=X_val, X_test=X_test)
        return (X_train_pca, labels_train, X_val_pca, labels_val, X_test_pca, labels_test)


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
FELPY

__author__ = "Trey Guest"
__credits__ = ["Trey Guest"]
__license__ = "EuXFEL"
__version__ = "0.2.1"
__maintainer__ = "Trey Guest"
__email__ = "trey.guest@xfel.eu"
__status__ = "Developement"
"""

import sys

import scipy

import numpy as np

from wpg import srwlib
from wpg.srw import srwlpy

from copy import deepcopy

from scipy.constants import c
from matplotlib import pyplot as plt

import wpg.srwlpy
from felpy.utils.opt_utils import ekev2wav
from felpy.model.wavefront import Wavefront
from wpg.wpg_uti_wf import look_at_q_space, calculate_fwhm, calc_pulse_energy
from wpg.wpg_uti_wf import plot_intensity_map as plotIntensity
from wpg.wpg_uti_wf import plot_intensity_qmap as plotKspace

from wpg.generators import build_gauss_wavefront_xy as build_gaussian
from wpg.generators import build_gauss_wavefront as build_gaussian_3D

from felpy.model.beamline import Beamline

from wpg.srwlib import SRWLOptL as thinLens

from wpg.optical_elements import Drift
from wpg.srwlib import SRWLOptD

from felpy.backend.wpg_converters import wavefront_from_array

fwhm2rms = np.sqrt(8*np.log(2)) ### FWHM = sqrt(8ln(2))*sigma
h = scipy.constants.physical_constants['Planck constant in eV s'][0]

def get_pulse_energy(wfr):
    wfr.set_electric_field_representation('t')
    energy = calc_pulse_energy(wfr)[0]
    wfr.set_electric_field_representation('f')
    return energy

def analytical_pulse_energy(q, ekev):
    """
    Estimate of analytical_pulse_energy from electron bunch charge and radiation energy

    :param q: electron bunch charge [nC]
    :param ekev: radiation energy [keV]

    :return P: pulse energy [J]
    """

    P = 19*q/ekev
    return P/1e3

def analytical_pulse_duration(q):
    """
    Estimate analytical_pulse_duration from electron bunch charge

    :param q: electron bunch charge [nC]

    :return t: Duration of pulse [s]
    """

    t = (q*1e3)/9.8
    return t*1e-15

def analytical_pulse_width(ekev):
    """
    Estimate analytical_pulse_width (FWHM) from radiation energy (assumes symmetrical beam)

    :param ekev: radiation energy [keV]

    :return sig: Radiation pulse width (FWHM) [m]
    """

    sig = np.log((7.4e03/ekev))*6
    return sig/1e6

def analytical_pulse_divergence(q, ekev):
    """
    Estimate of analytical_pulse_divergence (half-angle) from electron bunch charge and radiation energy

    :param q: electron bunch charge [nC]
    :param ekev: radiation energy [keV]

    :return dtheta: pulse divergence [rad]
    """

    dtheta = (17.2*6.4*np.sqrt(q))/ekev**(0.85)
    return dtheta/1e6

def modify_beam_divergence(wfr, sig, dtheta):
    """
    Construct a thin-lens with a focal length defined by a desired divergence,
    see: thin_lens_mod, and then propagate to 2f

    :param wfr: WPG wfr structure
    :params sig: beam fwhm [m]
    :params dtheta: pulse divergence [rad]
    """
    f = sig/(np.tan(dtheta))

    tl = thinLens(f,f)
    bl = Beamline()
    bl.append(tl, [0,0,0,0,0,1,1,1,1,0,0,0])
    bl.append(Drift(2*f), [0,0,0,0,0,1,1,1,1,0,0,0])
    bl.propagate(wfr)

    wfr.params.Mesh.zCoord = 0

def modify_beam_energy(wfr, desired):

    actual = get_pulse_energy(wfr)
    diff = actual/desired

    wfr.data.arrEhor *= diff

def construct_gaussian(nx, ny, ekev, extent, sigX, sigY, divergence, xoff = 0, yoff = 0, tau = 1, mx = 0, my = 0, tiltX = 0, tiltY = 0):
    gsnBm = build_gaussian(nx, ny, ekev, xMin = extent[0], xMax = extent[1], yMin = extent[2], yMax = extent[3],
                           sigX = sigX, sigY = sigY, d2waist = 1, xoff = xoff, yoff = yoff, pulseTau = tau, _mx = mx, _my = my, tiltX= tiltX, tiltY = tiltY)

    wfr = Wavefront(gsnBm)
    wfr.params.wavelength = ekev2wav(ekev)

    modify_beam_divergence(wfr,sigX, divergence)
    return wfr

def construct_SA1_wavefront(nx, ny, ekev, q, xoff = 0, yoff = 0, mx = 0, my = 0,
                            tiltX = 0, tiltY = 0,
                            xMin = -400e-06, xMax = 400e-06,
                            yMin = -400e-06, yMax = 400e-06):
    """
    Construct a fully-coherent Gaussian source with properties related to the
    energy of the radiation and beam charge.

    Important points: pulseTau (coherence length) is set to tau (pulse length)
    in the fully coherent case

    :param nx: number of horizontal pixels [int]
    :param ny: number of vertical pixels [int]
    :param nz: number of slices [int]
    :param ekev: radiation energy [keV]
    :param q: electron beam bunch charge [nC]
    :param xoff: horizontal offset of beam maximum [m]
    :param yoff: vertical offset of beam maximum [m]
    :param modify_beam_divergenceergence: Choose to set non-diffraction limited div [bool]
    """
    wavelength = (h*c)/(ekev*1e3)

    sigX, sigY = analytical_pulse_width(ekev)/fwhm2rms, analytical_pulse_width(ekev)/fwhm2rms
    pulseEn = analytical_pulse_energy(q, ekev)

    dtheta = analytical_pulse_divergence(q, ekev)
    tau = analytical_pulse_duration(q)

    gsnBm = build_gaussian(nx, ny, ekev, xMin, xMax, yMin, yMax, sigX, sigY, 1, xoff = xoff, yoff = yoff, pulseTau = tau, _mx = mx, _my = my, tiltX= tiltX, tiltY = tiltY)

    wfr = Wavefront(gsnBm)
    wfr.params.wavelength = wavelength
    srwlib.srwl.SetRepresElecField(wfr._srwl_wf, 'f')

    modify_beam_divergence(wfr,sigX, analytical_pulse_divergence(q,ekev)) ## note 300 is hardcoded fix, should not stay past tue 13/01/21

    return wfr

def construct_SA1_pulse(nx, ny, nz, ekev, q, mx = 0, my = 0):
    """
    Construct a fully-coherent Gaussian source with properties related to the
    energy of the radiation and beam charge.

    Important points: pulseTau (coherence length) is set to tau (pulse length)
    in the fully coherent case

    :param nx: number of horizontal pixels [int]
    :param ny: number of vertical pixels [int]
    :param nz: number of slices [int]
    :param ekev: radiation energy [keV]
    :param q: electron beam bunch charge [nC]
    :param xoff: horizontal offset of beam maximum [m]
    :param yoff: vertical offset of beam maximum [m]
    :param modify_beam_divergenceergence: Choose to set non-diffraction limited div [bool]
    """
    wavelength = (h*c)/(ekev*1e3)

    xMin, xMax = -400e-06, 400e-06 #based on fwhm of 1 nC, 3 keV beam
    yMin, yMax = -400e-06, 400e-06 #based on fwhm of 1 nC, 3 keV beam

    sigX, sigY = analytical_pulse_width(ekev)/fwhm2rms, analytical_pulse_width(ekev)/fwhm2rms
    pulseEn = analytical_pulse_energy(q, ekev)

    dtheta = analytical_pulse_divergence(q, ekev)
    tau = analytical_pulse_duration(q)

    gsnBm = build_gaussian_3D(nx = nx,
                              ny = ny,
                              nz = nz,
                              ekev = ekev,
                              xMin = xMin, xMax = xMax,
                              yMin = yMin, yMax = yMax,
                              sigX = sigX, sigY = sigY,
                              d2waist = 1,
                              tau = tau,
                              pulseRange = 1,
                              pulseEn = 2*pulseEn,
                              _mx = mx,
                              _my = my)

    wfr = Wavefront(gsnBm)
    srwlib.srwl.SetRepresElecField(wfr._srwl_wf, 'f')

    wfr.params.wavelength = wavelength

    #modify_beam_divergence(wfr,analytical_pulse_width(ekev),analytical_pulse_divergence(q,ekev))

    modify_beam_energy(wfr, desired = 2*pulseEn)
    return wfr

def test_coherent_pulse_divergence(outdir = None):
    """
    Test that the FWHM of the source matches the analytical prediction

    :param outdir: output directory [str]
    """

    fig = plt.figure(figsize = [12,8])
    ax = fig.add_subplot(111)
    ax.set_title("Source Divergence FWHM", fontsize = 22)
    ax.set_ylabel(r"d$\theta$ [$\mu$rad]", fontsize = 16)
    ax.set_xlabel("Radiation Energy [keV]", fontsize = 16)
    ax.set_xlim([2.5,17.5])

    nx, ny = 512, 512

    esteps = 5
    energyrange = np.linspace(3,16, esteps)

    for q in [0.01, 0.1, 0.3, 0.5, 0.8, 1.0]:

        divergences = []
        for energy in energyrange:
            wfr = construct_SA1_wavefront(nx, ny, energy, q)
            divergences.append(wfr.get_divergence()[0])
        ax.scatter(energyrange, np.array(divergences)*1e6, marker = 'x')

        analytical_data = [analytical_pulse_divergence(q, energy) for energy in energyrange]
        ax.plot(energyrange, np.array(analytical_data)*1e6, '--')

    ax.legend(["0.01 nC", "0.1 nC", "0.3 nC", "0.5 nC", "0.8 nC", "1.0 nC"])

    if outdir == None:
        pass
    else:
        fig.savefig(outdir + "SourceDivergence_Energy.png")

if __name__ == '__main__':

    pass

max_stars_repo_name: twguest/FELPy

max_issues_repo_name: twguest/FELPy

max_forks_repo_name: twguest/FELPy

reason: import numpy,import scipy,from scipy


text:
#!/usr/bin/env python
"""
Generic python script.
"""
__author__ = "Alex Drlica-Wagner"
import unittest
# Execute tests in order: https://stackoverflow.com/a/22317851/4075339
unittest.TestLoader.sortTestMethodsUsing = None

import numpy as np
import fitsio

import ugali.analysis.loglike
from ugali.utils.logger import logger
logger.setLevel(logger.WARN)

CONFIG='tests/config.yaml'
LON = RA = 53.92
LAT = DEC = -54.05
IDX = [1,2537,9000]

class TestLoglike(unittest.TestCase):
    """Test the loglikelihood"""

    def setUp(self):
        self.loglike = ugali.analysis.loglike.createLoglike(CONFIG,lon=LON,lat=LAT)
        self.source = self.loglike.source
        self.filename = 'test-membership.fits'

    def test_initial_config(self):
        # Likelihood configuration.
        np.testing.assert_equal(self.source.isochrone.name,'Bressan2012')
        np.testing.assert_equal(self.source.kernel.name,'RadialPlummer')
        np.testing.assert_equal(self.source.richness,1000.)

    def test_initial_probability(self):
        # Probability calculations

        np.testing.assert_allclose(self.loglike.f,0.08595560,rtol=1e-6)
        np.testing.assert_allclose(self.loglike.u[IDX],
                                   [5.29605173e-03, 1.80040569e-03, 5.52283081e-09],
                                   rtol=1e-6)
        np.testing.assert_allclose(self.loglike.b[IDX],
                                   [4215.31143651, 9149.29106545, 1698.22182173],
                                   rtol=1e-6)
        np.testing.assert_allclose(self.loglike.p[IDX],
                                   [1.25480793e-03, 1.96742181e-04, 3.25212568e-09],
                                   rtol=1e-6)

        np.testing.assert_allclose(self.loglike(),3948.1559048)
        np.testing.assert_allclose(self.loglike.ts(),7896.31181)
        np.testing.assert_allclose(self.loglike.nobs,85.9556015)

    def test_fit_richness(self):
        # Fit the richness
        interval = self.loglike.richness_interval()
        np.testing.assert_allclose(interval,(31596.21551, 32918.707276))

        lnl,rich,para = self.loglike.fit_richness()
        np.testing.assert_allclose(lnl,8449.77225)
        np.testing.assert_allclose(rich,32252.807226)
        np.testing.assert_allclose(self.loglike.source.richness,rich)

    def test_write_membership(self):
        # Write membership
        self.loglike.write_membership(self.filename)

        # Read membership and metadata
        mem,hdr = fitsio.read(self.filename,header=True)

        np.testing.assert_allclose(self.loglike.p,mem['PROB'])
        np.testing.assert_allclose(self.loglike.catalog.ra,mem['RA'])
        np.testing.assert_allclose(self.loglike.catalog.dec,mem['DEC'])
        np.testing.assert_allclose(self.loglike.catalog.color,mem['COLOR'])

        # Testing output metadata
        print (self.loglike.ts(),hdr['TS'])
        np.testing.assert_allclose(self.loglike.ts(),hdr['TS'])
        np.testing.assert_allclose(self.source.richness,hdr['RICHNESS'])
        np.testing.assert_allclose(self.source.age,hdr['AGE'])
        np.testing.assert_allclose(self.source.z,hdr['METALLICITY'])
        np.testing.assert_allclose(self.source.distance_modulus,hdr['DISTANCE_MODULUS'])
        np.testing.assert_allclose(self.source.lon,hdr['LON'])
        np.testing.assert_allclose(self.source.lat,hdr['LAT'])

if __name__ == "__main__":
    unittest.main()

max_stars_repo_name: mcnanna/ugali

max_issues_repo_name: mcnanna/ugali

max_forks_repo_name: kadrlica/ugali


text:
import keras.backend as K
import numpy as np
import logging
import os
from keras.models import Model
from keras import optimizers
from layers import (CNNContextLayer, Dot_Abs, ESIMAttention, NormAttention, HingeLoss)
from char_model import CharRNN
from keras.layers import *

logger = logging.getLogger(__name__)
local_file = os.path.split(__file__)[-1]
logging.basicConfig(
    format='%(asctime)s : %(filename)s : %(funcName)s : %(levelname)s : %(message)s',
    level=logging.INFO)

def create_model(word_embedding_matrix, word_list, char_list, params):
    logger.info(params)
    words_num = len(word_list)
    char_num = len(char_list)
    word_embedding_dim = np.shape(word_embedding_matrix)[-1]
    sentence_length = params.sentence_length
    filters = params.filters
    kernel_sizes = [int(x) for x in params.kernel_sizes.split(',')]
    char_rnn_dim = params.char_rnn_dim
    dropout = params.dropout
    character_length = params.character_length
    char_dim = params.char_dim
    lr = params.lr
    decay_rate = params.decay_rate
    hinge = params.hinge
    coherence_k = params.voting_k * 2
    context_sentence_length = params.context_sentence_length
    context_rnn_dim = params.context_rnn_dim

    # Build input
    mention = Input(batch_shape=(None, sentence_length), dtype='int32', name='mention')
    candidate = Input(batch_shape=(None, sentence_length), dtype='int32', name='candidate')
    men_char = Input(batch_shape=(None, sentence_length, character_length), dtype='int32', name='men_char')
    can_char = Input(batch_shape=(None, sentence_length, character_length), dtype='int32', name='can_char')
    can_prior = Input(batch_shape=(None, 1), dtype='float32', name='can_prior')
    can_context = Input(batch_shape=(None, context_sentence_length), dtype='int32', name='can_context')
    candidate_emb = Input(batch_shape=(None, 50), dtype='float32', name='candidate_emb')
    voting_candidates_emb = Input(batch_shape=(None, coherence_k, 50), dtype='float32', name='voting_candidates_emb')

    pos_candidate = Input(batch_shape=(None, sentence_length), dtype='int32', name='pos_candidate')
    pos_can_char = Input(batch_shape=(None, sentence_length, character_length), dtype='int32', name='pos_can_char')
    neg_candidate = Input(batch_shape=(None, sentence_length), dtype='int32', name='neg_candidate')
    neg_can_char = Input(batch_shape=(None, sentence_length, character_length), dtype='int32', name='neg_can_char')
    pos_can_prior = Input(batch_shape=(None, 1), dtype='float32', name='pos_can_prior')
    neg_can_prior = Input(batch_shape=(None, 1), dtype='float32', name='neg_can_prior')
    pos_candidate_emb = Input(batch_shape=(None, 50), dtype='float32', name='pos_candidate_emb')
    neg_candidate_emb = Input(batch_shape=(None, 50), dtype='float32', name='neg_candidate_emb')

    word_layer = Embedding(
        input_dim=words_num,
        output_dim=word_embedding_dim,
        weights=[word_embedding_matrix],
        trainable=True
    )
    w_res1 = word_layer(mention)
    w_res2 = word_layer(candidate)

    # char embedding
    char_layer = Embedding(char_num, char_dim, input_length=(sentence_length, character_length),
                           trainable=True, mask_zero=True,
                           name='char_embedding', embeddings_initializer='he_uniform')

    char_res1 = char_layer(men_char)
    char_res2 = char_layer(can_char)
    char_shape = K.int_shape(char_res1)

    #char_repre
    char_rnn = CharRNN(char_shape, dropout, char_rnn_dim)

    # batch * sequenth length * rnn dim
    ch_res1 = char_rnn(char_res1)
    ch_res2 = char_rnn(char_res2)

    w_res1 = concatenate([w_res1, ch_res1])
    w_res2 = concatenate([w_res2, ch_res2])

    #attention layer
    shape = K.int_shape(w_res1)
    attention1, attention2 = ESIMAttention(shape)(w_res1, w_res2)

    shape = K.int_shape(attention1)
    mention_att1 = Dot_Abs(shape, axis=2)(w_res1, attention1)
    repre1 = concatenate([w_res1, attention1], axis=-1)
    com_repre1 = concatenate([repre1, mention_att1], axis=-1)

    mention_att2 = Dot_Abs(shape, axis=2)(w_res2, attention2)
    repre2 = concatenate([w_res2, attention2], axis=-1)
    com_repre2 = concatenate([repre2, mention_att2], axis=-1)

    # Cnn context layer
    w_res1_shape = K.int_shape(com_repre1)
    cnns = [CNNContextLayer(filters=filters, kernel_size=kernel_size, input_shape=(w_res1_shape[1], w_res1_shape[2]),
                            dropout=dropout) for kernel_size in kernel_sizes]

    sequence1 = concatenate([cnn(com_repre1) for cnn in cnns])
    sequence2 = concatenate([cnn(com_repre2) for cnn in cnns])
    sequence = concatenate([sequence1, sequence2])

    #prior
    if params.add_prior:
        sequence = concatenate([sequence, can_prior])

    #context
    context_rnn1 = Bidirectional(GRU(context_rnn_dim))
    can_repre = word_layer(candidate)
    can_repre = context_rnn1(can_repre)
    context_rnn2 = Bidirectional(GRU(context_rnn_dim, return_sequences=True))
    context_repre = word_layer(can_context)
    context_repre = context_rnn2(context_repre)
    norm_att = NormAttention(K.int_shape(context_repre))
    weight = norm_att(can_repre, context_repre)
    context_repre = Dot(axes=1)([weight, context_repre])
    context_score = Dot(axes=1, normalize=True)([can_repre, context_repre])
    if params.add_context:
        sequence = concatenate([sequence, context_score])

    #coherence
    candidate_embs = RepeatVector(K.int_shape(voting_candidates_emb)[1])(candidate_emb)
    coherence_feature = Dot(normalize=True, axes=-1)([candidate_embs, voting_candidates_emb])
    pool = GlobalAveragePooling1D()(coherence_feature)
    coherence_score = Lambda(K.mean, output_shape=(1,), arguments={'axis': -1, 'keepdims':True})(pool)
    if params.add_coherence:
        sequence = concatenate([sequence, coherence_score])

    #output
    sequence = Dropout(dropout)(sequence)
    dense = Dense(64, activation='relu', input_shape=K.int_shape(sequence))(sequence)
    dense = Dropout(dropout)(dense)
    pred = Dense(1, activation='sigmoid', name='final_dense')(dense)

    #hinge loss
    score_model = Model(
        inputs=[mention, candidate, men_char, can_char, can_prior, candidate_emb, voting_candidates_emb, can_context], outputs=pred)

    pos_score = score_model(
        [mention, pos_candidate, men_char, pos_can_char, pos_can_prior, pos_candidate_emb, voting_candidates_emb, can_context])
    neg_score = score_model(
        [mention, neg_candidate, men_char, neg_can_char, neg_can_prior, neg_candidate_emb, voting_candidates_emb, can_context])

    result = HingeLoss(alpha=hinge)(pos_score, neg_score)

    train_inputs = (
    mention, pos_candidate, men_char, pos_can_char, neg_candidate, neg_can_char, pos_can_prior, neg_can_prior,
    pos_candidate_emb, neg_candidate_emb, voting_candidates_emb, can_context)
    train_model = Model(inputs=train_inputs, outputs=result)

    # Compile model
    nadam = optimizers.nadam(lr=lr, schedule_decay=decay_rate)
    train_model.compile(loss=lambda y_true, y_pred: y_pred, optimizer=nadam)
    pre_inputs = (mention, candidate, men_char, can_char, can_prior, candidate_emb, voting_candidates_emb, can_context)
    predict_model = Model(inputs=pre_inputs, outputs=pred)
    predict_model.compile(loss=lambda y_true, y_pred: y_pred, optimizer=nadam)

    return train_model, predict_model

if __name__ == '__main__':
    pass


text:
from path import path_code_dir
import sys
sys.path.insert(0, path_code_dir)
from amftrack.util import get_dates_datetime, get_dirname
import scipy.io as sio
from pymatreader import read_mat
import numpy as np
from amftrack.pipeline.functions.image_processing.extract_graph import (
    sparse_to_doc,
    get_degree3_nodes,
)
import open3d as o3d
from cycpd import rigid_registration
from amftrack.pipeline.paths.directory import directory_scratch
import pandas as pd

i = int(sys.argv[-1])
op_id = int(sys.argv[-2])
thresh = int(sys.argv[1])
directory = str(sys.argv[2])

run_info = pd.read_json(f'{directory_scratch}temp/{op_id}.json')
folder_list = list(run_info['folder'])
folder_list.sort()

dates_datetime_chosen = folder_list[i : i + 2]
print("========")
print(f"Matching plate {dates_datetime_chosen[0]} at dates {dates_datetime_chosen}")
print("========")
dates = dates_datetime_chosen

dilateds = []
skels = []
skel_docs = []
for directory_name in dates:
    path_snap = directory + directory_name
    skel_info = read_mat(path_snap + "/Analysis/skeleton_pruned.mat")
    skel = skel_info["skeleton"]
    skels.append(skel)
    skel_doc = sparse_to_doc(skel)
    skel_docs.append(skel_doc)
isnan = True
for order in [(0,1),(1,0)]:
    skeleton1, skeleton2 = skel_docs[order[0]], skel_docs[order[1]]
    skelet_pos = np.array(list(skeleton1.keys()))
    samples = np.random.choice(skelet_pos.shape[0], len(skeleton2.keys()) // 100)
    X = np.transpose(skelet_pos[samples, :])
    skelet_pos = np.array(list(skeleton2.keys()))
    samples = np.random.choice(skelet_pos.shape[0], len(skeleton2.keys()) // 100)
    Y = np.transpose(skelet_pos[samples, :])
    reg = rigid_registration(
        **{
            "X": np.transpose(X.astype(float)),
            "Y": np.transpose(Y.astype(float)),
            "scale": False, 'tolerance' : 1e-5, 'w' : 1e-5
        }
    )
    out = reg.register()
    Rfound = reg.R[0:2, 0:2]
    tfound = np.dot(Rfound, reg.t[0:2])
    if order == (0,1):
        t_init = -tfound
        Rot_init = Rfound
    else:
        Rot_init,t_init = np.linalg.inv(Rfound), np.dot(np.linalg.inv(Rfound),tfound)
    sigma2 = reg.sigma2
    if sigma2>=thresh:
        print("========")
        print(f"Failed to match plate {dates_datetime_chosen[0]} at dates {dates_datetime_chosen}")
        print("========")
        continue
    isnan = np.isnan(tfound[0])
    if isnan:
        continue
#     X = np.transpose(
#         np.array([pos1[node] for node in pruned1 if pruned1.degree(node) == 3])
#     )
#     Y = np.transpose(
#         np.array([pos2[node] for node in pruned2 if pruned2.degree(node) == 3])
#     )
    skeleton1, skeleton2 = skel_docs[0], skel_docs[1]
    X = np.transpose(
        np.array(get_degree3_nodes(skeleton1))
    )
    Y = np.transpose(
        np.array(get_degree3_nodes(skeleton2))
    )
    # fig=plt.figure(figsize=(10,9))
    # ax = fig.add_subplot(111)
    # ax.scatter(X[0,:],X[1,:])
    # ax.scatter(Y[0,:],Y[1,:])
#     Xex = np.transpose(np.transpose(np.dot(Rot_init, X)) + t_init)
    # fig=plt.figure(figsize=(10,9))
    # ax = fig.add_subplot(111)
    # ax.scatter(Xex[0,:],Xex[1,:])
    # ax.scatter(Y[0,:],Y[1,:])
    X = np.insert(X, 2, values=0, axis=0)
    Y = np.insert(Y, 2, values=0, axis=0)
    print(X.shape, Y.shape)
    vectorX = o3d.utility.Vector3dVector(np.transpose(X))
    vectorY = o3d.utility.Vector3dVector(np.transpose(Y))
    source = o3d.geometry.PointCloud(vectorX)
    target = o3d.geometry.PointCloud(vectorY)
    threshold = 200
    trans_init = np.asarray(
        [
            [Rot_init[0, 0], Rot_init[0, 1], 0, t_init[0]],
            [Rot_init[1, 0], Rot_init[1, 1], 0, t_init[1]],
            [0, 0, 1, 0],
            [0.0, 0.0, 0.0, 1.0],
        ]
    )
    reg_p2p = o3d.registration.registration_icp(
        source,
        target,
        threshold,
        trans_init,
        o3d.registration.TransformationEstimationPointToPoint(),
    )
    print(reg_p2p)
    Rfound = reg_p2p.transformation[0:2, 0:2]
    tfound = reg_p2p.transformation[0:2, 3]
    print(Rfound, tfound)
    X, Y = X[0:2, :], Y[0:2, :]
    Yrep = np.transpose(np.transpose(np.dot(Rfound, X)) + tfound)
    # fig=plt.figure(figsize=(10,9))
    # ax = fig.add_subplot(111)
    # ax.scatter(np.transpose(Yrep)[:,0],np.transpose(Yrep)[:,1])
    # ax.scatter(np.transpose(Y)[:,0],np.transpose(Y)[:,1])
    break

if not isnan:
    sio.savemat(path_snap + "/Analysis/transform.mat", {"R": Rfound, "t": tfound})
else :
    sio.savemat(path_snap + "/Analysis/transform_corrupt.mat", {"R": np.array([[1,0],[0,1]]), "t": np.array([0,0])})

max_stars_repo_name: Cocopyth/MscThesis

max_issues_repo_name: Cocopyth/MscThesis

max_forks_repo_name: Cocopyth/MscThesis

reason: import numpy,import scipy


text:
# Copyright (c) 2017, Apple Inc. All rights reserved.
#
# Use of this source code is governed by a BSD-3-clause license that can be
# found in the LICENSE.txt file or at https://opensource.org/licenses/BSD-3-Clause

import numpy as _np

from . import datatypes
from ._feature_management import process_or_validate_features
from ._feature_management import is_valid_feature_list
from . import _feature_management as _fm

def set_classifier_interface_params(spec, features, class_labels,
        model_accessor_for_class_labels, output_features = None, training_features=None):
    """
    Common utilities to set the regression interface params.
    """
    # Normalize the features list.
    features = _fm.process_or_validate_features(features)

    if class_labels is None:
        raise ValueError("List of class labels must be provided.")

    n_classes = len(class_labels)

    output_features = _fm.process_or_validate_classifier_output_features(output_features, class_labels)

    if len(output_features) == 1:
        predicted_class_output, pred_cl_type = output_features[0]
        score_output = None
    elif len(output_features) == 2:
        predicted_class_output, pred_cl_type = output_features[0]
        score_output, score_output_type = output_features[1]
    else:
        raise ValueError("Provided output classes for a classifier must be "
                "a list of features, predicted class and (optionally) class_score.")

    spec.description.predictedFeatureName = predicted_class_output

    # Are they out of order?
    if not (pred_cl_type == datatypes.Int64() or pred_cl_type == datatypes.String()):
        raise ValueError("Provided predicted class output type not Int64 or String (%s)."
                % repr(pred_cl_type))

    if score_output is not None:
        if not isinstance(score_output_type, datatypes.Dictionary):
            raise ValueError("Provided class score output type not a Dictionary (%s)."
                    % repr(score_output_type))

        if score_output_type.key_type != pred_cl_type:
            raise ValueError(("Provided class score output (%s) key_type (%s) does not "
                    "match type of class prediction (%s).")
                    % (score_output, repr(score_output_type.key_type), repr(pred_cl_type)))

        spec.description.predictedProbabilitiesName = score_output

    # add input
    for index, (cur_input_name, input_type) in enumerate(features):
        input_ = spec.description.input.add()
        input_.name = cur_input_name
        datatypes._set_datatype(input_.type, input_type)

    # add output
    for index, (cur_output_name, output_type) in enumerate(output_features):
        output_ = spec.description.output.add()
        output_.name = cur_output_name
        datatypes._set_datatype(output_.type, output_type)

    # Add training features
    if training_features is not None:
        spec = set_training_features(spec, training_features)

    # Worry about the class labels
    if pred_cl_type == datatypes.String():
        try:
            for c in class_labels:
                getattr(spec, model_accessor_for_class_labels).stringClassLabels.vector.append(str(c))
        # Not all the classifiers have class labels; in particular the pipeline
        # classifier.  Thus it's not an error if we can't actually set them.
        except AttributeError:
            pass

    else:
        for c in class_labels:
            conv_error = False
            try:
                if not (int(c) == c):
                    conv_error = True
            except:
                conv_error = True

            if conv_error:
                raise TypeError(("Cannot cast '%s' class to an int type " % str(c))
                + "(class type determined by type of first class).")

            try:
                getattr(spec, model_accessor_for_class_labels).int64ClassLabels.vector.append(int(c))
            # Not all the classifiers have class labels; in particular the pipeline
            # classifier.  Thus it's not an error if we can't actually set them.
            except AttributeError:
                break

    # And we are done!
    return spec

def set_regressor_interface_params(spec, features, output_features, training_features=None):
    """ Common utilities to set the regressor interface params.
    """
    if output_features is None:
        output_features = [("predicted_class", datatypes.Double())]
    else:
        output_features = _fm.process_or_validate_features(output_features, 1)

    if len(output_features) != 1:
        raise ValueError("Provided output features for a regressor must be "
                    "one Double feature.")

    if output_features[0][1] != datatypes.Double():
        raise ValueError("Output type of a regressor must be a Double.")

    prediction_name = output_features[0][0]
    spec.description.predictedFeatureName = prediction_name

    # Normalize the features list.
    features = _fm.process_or_validate_features(features)

    # add input and output features
    for cur_input_name, feature_type in features:
        input_ = spec.description.input.add()
        input_.name = cur_input_name
        datatypes._set_datatype(input_.type, feature_type)

    # Add training features
    if training_features is not None:
        spec = set_training_features(spec, training_features)

    output_ = spec.description.output.add()
    output_.name = prediction_name
    datatypes._set_datatype(output_.type, 'Double')
    return spec

def set_transform_interface_params(spec, input_features, output_features, are_optional=False, training_features=None):
    """ Common utilities to set transform interface params.
    """
    input_features = _fm.process_or_validate_features(input_features)
    output_features = _fm.process_or_validate_features(output_features)

    # Add input and output features
    for (fname, ftype) in input_features:
        input_ = spec.description.input.add()
        input_.name = fname
        datatypes._set_datatype(input_.type, ftype)
        if are_optional:
            input_.type.isOptional = are_optional

    for (fname, ftype) in output_features:
        output_ = spec.description.output.add()
        output_.name = fname
        datatypes._set_datatype(output_.type, ftype)

    # Add training features
    if training_features is not None:
        spec = set_training_features(spec, training_features)

    return spec

def set_training_features(spec, training_features):

    for (fname, ftype) in training_features:
        training_input_ = spec.description.trainingInput.add()
        training_input_.name = fname
        if ftype:
            datatypes._set_datatype(training_input_.type, ftype)

    return spec

max_stars_repo_name: Gerzer/coremltools

max_issues_repo_name: Gerzer/coremltools

max_forks_repo_name: Gerzer/coremltools


text:
"""
Main file for training multi-camera pose
"""

import os
import optparse
import time
import cPickle as pickle
import numpy as np
from skimage import color

from pyKinectTools.dataset_readers.KinectPlayer import KinectPlayer, display_help
from pyKinectTools.utils.DepthUtils import world2depth, depthIm2XYZ, skel2depth, depth2world
from pyKinectTools.utils.SkeletonUtils import display_skeletons, transform_skels, kinect_to_msr_skel
# from pyKinectTools.algs.GeodesicSkeleton import generateKeypoints, distance_map
from pyKinectTools.algs.RandomForestPose import RFPose

""" Debugging """
from IPython import embed
np.seterr(all='ignore')

# -------------------------MAIN------------------------------------------

def main(visualize=False, learn=False):
	# Init both cameras
	# fill = True
	fill = False
	get_color = True
	cam = KinectPlayer(base_dir='./', device=1, bg_subtraction=True, get_depth=True, get_color=get_color, get_skeleton=True, fill_images=fill)
	cam2 = KinectPlayer(base_dir='./', device=2, bg_subtraction=True, get_depth=True, get_color=get_color, get_skeleton=True, fill_images=fill)
	# Transformation matrix from first to second camera
	data = pickle.load(open("Registration.dat", 'r'))
	transform_c1_to_c2 = data['transform']

	# Get random forest parameters
	if learn:
		rf = RFPose(offset_max=100, feature_count=300)
	else:
		rf = RFPose()
		rf.load_forest()

	ii = 0
	# cam.next(200)
	all_joint_ims_z = []
	all_joint_ims_c = []
	while cam.next() and ii < 200:
		# Update frames
		cam2.sync_cameras(cam)
		if ii%10 != 0:
			ii += 1
			continue

		# Transform skels from cam1 to cam2
		cam_skels = [np.array(cam.users[s]['jointPositions'].values()) for s in cam.users.keys()]
		# Get rid of bad skels
		cam_skels = [s for s in cam_skels if np.all(s[0] != -1)]

		if len(cam_skels) == 0:
			continue
		ii+=1

		# Save images
		if 1:
			joint_ims_z = []
			joint_ims_c = []
			dx = 10
			skel_tmp = skel2depth(cam_skels[0], [240,320])
			for j_pos in skel_tmp:
				# embed()
				joint_ims_z += [cam.depthIm[j_pos[0]-dx:j_pos[0]+dx, j_pos[1]-dx:j_pos[1]+dx]]
				joint_ims_c += [cam.colorIm[j_pos[0]-dx:j_pos[0]+dx, j_pos[1]-dx:j_pos[1]+dx]]
			if len(joint_ims_z) > 0:
				all_joint_ims_z += [joint_ims_z]
				all_joint_ims_c += [joint_ims_c]

		if 0:
			cam2_skels = transform_skels(cam_skels, transform_c1_to_c2, 'image')

			try:
				depth = cam2.get_person()
				if learn:
					rf.add_frame(depth, cam2_skels[0])
				else:
					rf.infer_pose(depth)

				if visualize:
					cam2.depthIm = display_skeletons(cam2.depthIm, cam2_skels[0], (5000,), skel_type='Low')
					skel1 = kinect_to_msr_skel(skel2depth(cam_skels[0], [240,320]))
					cam.depthIm = display_skeletons(cam.depthIm, skel1, (5000,), skel_type='Low')
					cam.visualize()
					cam2.visualize()
			except:
				pass

	embed()
	if learn:
		print "Starting forest"
		rf.learn_forest()

	print 'Done'

if __name__=="__main__":

	parser = optparse.OptionParser()
	parser.add_option('-v', '--visualize', dest='viz', action="store_true", default=False, help='Enable visualization')
	parser.add_option('-l', '--learn', dest='learn', action="store_true", default=False, help='Training phase')
	parser.add_option('-i', '--infer', dest='infer', action="store_true", default=False, help='Training phase')
	(opt, args) = parser.parse_args()

	if len(args) > 0:
		print "Wrong input argument"
	elif opt.infer==False and opt.learn==False:
		print "You must supply the program with either -l (learn parameters) or -i (infer pose)."
	else:
		if opt.infer:
			main(visualize=opt.viz, learn=False)
		else:
			main(visualize=opt.viz, learn=True)

	# '''Profiling'''
	# cProfile.runctx('main()', globals(), locals(), filename="ShowSkeletons.profile")

max_stars_repo_name: colincsl/pyKinectTools

max_issues_repo_name: colincsl/pyKinectTools

max_forks_repo_name: colincsl/pyKinectTools


text:
from neural_deprojection.models.identify_medium.generate_data import decode_examples
from neural_deprojection.graph_net_utils import vanilla_training_loop, TrainOneEpoch, AbstractModule, \
    get_distribution_strategy, build_log_dir, build_checkpoint_dir
import glob, os
import tensorflow as tf
from functools import partial
import pylab as plt

from graph_nets import blocks
import sonnet as snt
from graph_nets.graphs import GraphsTuple
from graph_nets.utils_tf import fully_connect_graph_dynamic
import json

class RelationNetwork(AbstractModule):
    """Implementation of a Relation Network.

    See https://arxiv.org/abs/1706.01427 for more details.

    The global and edges features of the input graph are not used, and are
    allowed to be `None` (the receivers and senders properties must be present).
    The output graph has updated, non-`None`, globals.
    """

    def __init__(self,
                 edge_model_fn,
                 global_model_fn,
                 reducer=tf.math.unsorted_segment_mean,
                 use_globals=False,
                 name="relation_network"):
        """Initializes the RelationNetwork module.

        Args:
          edge_model_fn: A callable that will be passed to EdgeBlock to perform
            per-edge computations. The callable must return a Sonnet module (or
            equivalent; see EdgeBlock for details).
          global_model_fn: A callable that will be passed to GlobalBlock to perform
            per-global computations. The callable must return a Sonnet module (or
            equivalent; see GlobalBlock for details).
          reducer: Reducer to be used by GlobalBlock to aggregate edges. Defaults to
            tf.math.unsorted_segment_sum.
          name: The module name.
        """
        super(RelationNetwork, self).__init__(name=name)

        self._edge_block = blocks.EdgeBlock(
            edge_model_fn=edge_model_fn,
            use_edges=False,
            use_receiver_nodes=True,
            use_sender_nodes=True,
            use_globals=use_globals)

        self._global_block = blocks.GlobalBlock(
            global_model_fn=global_model_fn,
            use_edges=True,
            use_nodes=False,
            use_globals=use_globals,
            edges_reducer=reducer)

    def _build(self, graph:GraphsTuple) -> GraphsTuple:
        """Connects the RelationNetwork.

        Args:
          graph: A `graphs.GraphsTuple` containing `Tensor`s, except for the edges
            and global properties which may be `None`.

        Returns:
          A `graphs.GraphsTuple` with updated globals.

        Raises:
          ValueError: If any of `graph.nodes`, `graph.receivers` or `graph.senders`
            is `None`.
        """
        output_graph = self._global_block(self._edge_block(graph))
        return output_graph  # graph.replace(globals=output_graph.globals)

class Model(AbstractModule):
    def __init__(self, image_feature_size=16, kernel_size=3, name=None, **unused_kwargs):
        super(Model, self).__init__(name=name)
        self.encoder_graph = RelationNetwork(lambda: snt.nets.MLP([32,16], activate_final=True),
                                       lambda: snt.nets.MLP([32, 16], activate_final=True))
        self.encoder_image = RelationNetwork(lambda: snt.nets.MLP([32, 16], activate_final=True),
                                       lambda: snt.nets.MLP([32, 16], activate_final=True))
        self.image_cnn = snt.Sequential([snt.Conv2D(16,kernel_size, stride=2), tf.nn.relu,
                                         snt.Conv2D(image_feature_size, kernel_size, stride=2), tf.nn.relu])
        self.compare = snt.nets.MLP([32, 1])
        self.image_feature_size=image_feature_size

        self._step = None

    @property
    def step(self):
        if self._step is None:
            raise ValueError("Need to set step idx variable. model.step = epoch")
        return self._step

    @step.setter
    def step(self, value):
        self._step = value

    def _build(self, batch, *args, **kwargs):
        (graph, img, c) = batch
        del c
        encoded_graph = self.encoder_graph(graph)
        tf.summary.image(f'img_before_cnn', img[None,...], step=self.step)
        img = self.image_cnn(img[None,...])
        for channel in range(img.shape[-1]):
            tf.summary.image(f'img_after_cnn[{channel}]', img[...,channel:channel+1], step=self.step)
        #1, w,h,c -> w*h, c
        nodes = tf.reshape(img, (-1,self.image_feature_size))
        img_graph = GraphsTuple(nodes=nodes,
                            edges=None,
                            globals=None,
                            receivers=None,
                            senders=None,
                            n_node=tf.shape(nodes)[0:1],
                            n_edge=tf.constant([0]))
        connected_graph = fully_connect_graph_dynamic(img_graph)
        encoded_img = self.encoder_image(connected_graph)
        return self.compare(tf.concat([encoded_graph.globals, encoded_img.globals], axis=1))#[1]

MODEL_MAP = dict(model1=Model)

def build_training(model_type, model_parameters, optimizer_parameters, loss_parameters, strategy=None) -> TrainOneEpoch:
    model_cls = MODEL_MAP[model_type]

    model = model_cls(**model_parameters)

    def build_opt(**kwargs):
        opt_type = kwargs.get('opt_type')
        if opt_type == 'adam':
            learning_rate = kwargs.get('learning_rate', 1e-4)
            opt = snt.optimizers.Adam(learning_rate)
        else:
            raise ValueError('Opt {} invalid'.format(opt_type))
        return opt

    def build_loss(**loss_parameters):
        def loss(model_outputs, batch):
            (graph, img, c) = batch
            # loss =  mean(-sum_k^2 true[k] * log(pred[k]/true[k]))
            return tf.reduce_mean(tf.losses.binary_crossentropy(c[None,None],model_outputs, from_logits=True))
        return loss

    loss = build_loss(**loss_parameters)
    opt = build_opt(**optimizer_parameters)

    training = TrainOneEpoch(model, loss, opt, strategy=strategy)

    return training

def build_dataset(data_dir):
    """
    Build data set from a directory of tfrecords.

    Args:
        data_dir: str, path to *.tfrecords

    Returns: Dataset obj.
    """
    tfrecords = glob.glob(os.path.join(data_dir, '*.tfrecords'))
    dataset = tf.data.TFRecordDataset(tfrecords).map(partial(decode_examples,
                                                             node_shape=(4,),
                                                             image_shape=(18, 18, 1)))  # (graph_data_dict, image, idx)
    _graphs = dataset.map(lambda graph_data_dict, img, idx: (graph_data_dict, idx)).shuffle(buffer_size=50)
    _images = dataset.map(lambda graph_data_dict, img, idx: (img, idx)).shuffle(buffer_size=50)
    shuffled_dataset = tf.data.Dataset.zip((_graphs, _images))  # ((graph_data_dict, idx1), (img, idx2))
    shuffled_dataset = shuffled_dataset.map(lambda ds1, ds2: (ds1[0], ds2[0], ds1[1] == ds2[1]))  # (graph, img, yes/no)
    shuffled_dataset = shuffled_dataset.filter(lambda graph, img, c: ~c)
    shuffled_dataset = shuffled_dataset.map(lambda graph_data_dict, img, c: (graph_data_dict, img, tf.cast(c, tf.int32)))
    nonshuffeled_dataset = dataset.map(
        lambda graph_data_dict, img, idx: (graph_data_dict, img, tf.constant(1, dtype=tf.int32)))  # (graph, img, yes)
    dataset = tf.data.experimental.sample_from_datasets([shuffled_dataset, nonshuffeled_dataset])
    dataset = dataset.map(lambda graph_data_dict, img, c: (GraphsTuple(globals=None, edges=None, **graph_data_dict), img, c))
    return dataset

def main(data_dir, config):
    # Make strategy at the start of your main before any other tf code is run.
    strategy = get_distribution_strategy(use_cpus=True, logical_per_physical_factor=1)

    train_dataset = build_dataset(os.path.join(data_dir,'train'))
    test_dataset = build_dataset(os.path.join(data_dir,'test'))

    for (graph, img, c) in iter(test_dataset):
        print(graph)
        break

    with strategy.scope():
        train_one_epoch = build_training(**config)

    log_dir = build_log_dir('test_log_dir', config)
    checkpoint_dir = build_checkpoint_dir('test_checkpointing', config)
    with open(os.path.join(checkpoint_dir,'config.json'), 'w') as f:
        json.dump(config, f)

    vanilla_training_loop(train_one_epoch=train_one_epoch,
                          training_dataset=train_dataset,
                          test_dataset=test_dataset,
                          num_epochs=1000,
                          early_stop_patience=3,
                          checkpoint_dir=checkpoint_dir,
                          log_dir=log_dir,
                          debug=False)

if __name__ == '__main__':
    import numpy as np
    for kernel_size in [3,4,5]:
        for learning_rate in 10**np.linspace(-5, -3,5):
            for image_feature_size in [4,8,16]:
                config = dict(model_type='model1',
                              model_parameters=dict(kernel_size=kernel_size, image_feature_size=image_feature_size),
                              optimizer_parameters=dict(learning_rate=learning_rate, opt_type='adam'),
                              loss_parameters=dict())
                main('data', config)


text:
# -*- coding: utf-8 -*-
"""
Created on 2021-2022

@author: Rob
"""

import os
import pickle
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime
from pyproj import Proj, Geod
from keras.models import Model, load_model, Sequential
from keras.layers import Dense, Dropout, LSTM, Input, RepeatVector, GRU, Bidirectional
from keras.callbacks import EarlyStopping, ModelCheckpoint, ReduceLROnPlateau
from maeh_extension import MAEH_Extension as ext

"""
    LOADING DATA AND PREPARE STRUCTURE (FLAT)
"""
# load preprocessed data by experiments
X_train = np.load("X_train.npy")
X_val = np.load("X_val.npy")
X_test = np.load("X_test.npy")
y_train = np.load("y_train.npy")
y_val = np.load("y_val.npy")
y_test = np.load("y_test.npy")

# depends by experiment
X_train = X_train[:,:,:5] # drop d/Lat and d/Lon - FOR EXP-A ONLY
X_val = X_val[:,:,:5] # drop d/Lat and d/Lon - FOR EXP-A ONLY
X_test = X_test[:,:,:5] # drop d/Lat and d/Lon - FOR EXP-A ONLY
y_train = y_train[:,:,:2] # drop all features except lat and long
y_val = y_val[:,:,:2] # drop all features except lat and long
y_test = y_test[:,:,:2] # drop all features except lat and long

# flatten output structure
y_train_flat = np.copy(y_train)
y_val_flat = np.copy(y_val)
y_test_flat = np.copy(y_test)
y_train_flat = y_train_flat.reshape(y_train_flat.shape[0], y_train_flat.shape[1]*y_train_flat.shape[2])
y_val_flat = y_val_flat.reshape(y_val_flat.shape[0], y_val_flat.shape[1]*y_val_flat.shape[2])
y_test_flat = y_test_flat.reshape(y_test_flat.shape[0], y_test_flat.shape[1]*y_test_flat.shape[2])

"""
    HYPERPARAMETERS
"""
geod = Geod(ellps='WGS84', proj="utm", zone=33)
p = Proj (
    proj="utm",
    zone=33,
    ellps="WGS84"
)

# experiment parameters
lstm_cells = [25,50,75,100,125,150,175,200,225,250,275,300]
model_repeat_time = 10

seq_input_length = 30;
seq_output_length = 20;

# network parameters
n_features = X_train.shape[2]
lstm_layer_size = 25
num_epochs = 100
batch_size = 128

"""
    DEEP LEARNING ARCHITECTURES
"""
# define LSTM AE model (Sequence-to-sequence autoencoder)
def lstm_autoencoder_flat(rewrite_model, index):
    model_file_name = f'LSTM-AE_cells-{lstm_layer_size}_epochs-{num_epochs}_batch-{batch_size}_Test-{index}.h5'
    history_file_name = f'LSTM-AE_cells-{lstm_layer_size}_epochs-{num_epochs}_batch-{batch_size}_Test-{index}_history.npy'
    if(os.path.isfile(model_file_name) and not rewrite_model):
        autoencoder = load_model(model_file_name)
        history = np.load(history_file_name, allow_pickle=True).item()
    else:
        # define encoder
        inputs = Input(shape=(seq_input_length, n_features))
        encoder = LSTM(lstm_layer_size, activation='relu', return_sequences=False)(inputs)
        encoder = Dropout(0.01)(encoder)
        # define reconstruct decoder
        decoder = RepeatVector(seq_output_length)(encoder)
        decoder = LSTM(lstm_layer_size, activation='relu', return_sequences=False)(decoder)
        decoder = Dropout(0.01)(decoder)
        decoder = Dense(y_train_flat.shape[1], activation='relu')(decoder)
        # create autoencoder
        if(os.path.isfile(model_file_name)):
           autoencoder = load_model(model_file_name)
        else:
            autoencoder = Model(inputs, decoder)
        autoencoder.compile(optimizer='adam', loss='mse')

        reduce_lr = ReduceLROnPlateau(monitor="val_loss", factor=0.1, patience=8, min_lr=0.0001)
        earlyStop = EarlyStopping(monitor="val_loss", verbose=1, mode='min', patience=8)
        checkpoint = ModelCheckpoint(model_file_name, monitor='val_loss', verbose=1, save_best_only=True, mode='min')
        # fit model
        history = autoencoder.fit(X_train, y_train_flat, epochs=num_epochs, batch_size=batch_size, verbose=1, validation_data=(X_val, y_val_flat),
                                  callbacks=[reduce_lr,earlyStop,checkpoint])
        np.save(history_file_name, history.history)
    return autoencoder, history

def basic_gru_flat(rewrite_model, index):
    model_file_name = f'Basic-GRU_cells-{lstm_layer_size}_epochs-{num_epochs}_batch-{batch_size}_Test-{index}.h5'
    history_file_name = f'Basic-GRU_cells-{lstm_layer_size}_epochs-{num_epochs}_batch-{batch_size}_Test-{index}_history.npy'
    model = Sequential()
    if(os.path.isfile(model_file_name) and not rewrite_model):
        model = load_model(model_file_name)
        history = np.load(history_file_name, allow_pickle=True).item()
    else:
        if(os.path.isfile(model_file_name)):
           model = load_model(model_file_name)
        else:
            model = Sequential()

        model.add(GRU(lstm_layer_size, activation='relu', input_shape=(seq_input_length, n_features)))
        model.add(Dropout(0.01))

        model.add(Dense(y_train_flat.shape[1], activation='relu'))
        model.compile(optimizer='adam', loss='mse')

        reduce_lr = ReduceLROnPlateau(monitor="val_loss", factor=0.1, patience=8, min_lr=0.0001)
        earlyStop = EarlyStopping(monitor="val_loss", verbose=1, mode='min', patience=8)
        checkpoint = ModelCheckpoint(model_file_name, monitor='val_loss', verbose=1, save_best_only=True, mode='min')
        # fit model
        history = model.fit(X_train, y_train_flat, epochs=num_epochs, batch_size=batch_size, verbose=1, validation_data=(X_val, y_val_flat),
                            callbacks=[reduce_lr,earlyStop,checkpoint])
        np.save(history_file_name, history.history)
    return model, history

def basic_bidirectional_lstm_flat(rewrite_model, index):
    model_file_name = f'Basic-Bidirectional-LSTM_cells-{lstm_layer_size}_epochs-{num_epochs}_batch-{batch_size}_Test-{index}.h5'
    history_file_name = f'Basic-Bidirectional-LSTM_cells-{lstm_layer_size}_epochs-{num_epochs}_batch-{batch_size}_Test-{index}_history.npy'
    model = Sequential()
    if(os.path.isfile(model_file_name) and not rewrite_model):
        model = load_model(model_file_name)
        history = np.load(history_file_name, allow_pickle=True).item()
    else:
        if(os.path.isfile(model_file_name)):
           model = load_model(model_file_name)
        else:
            model = Sequential()

        model.add(Bidirectional(LSTM(lstm_layer_size, activation='relu'), input_shape=(seq_input_length, n_features)))
        model.add(Dropout(0.01))

        model.add(Dense(y_train_flat.shape[1], activation='relu'))
        model.compile(optimizer='adam', loss='mse')

        reduce_lr = ReduceLROnPlateau(monitor="val_loss", factor=0.1, patience=8, min_lr=0.0001)
        earlyStop = EarlyStopping(monitor="val_loss", verbose=1, mode='min', patience=8)
        checkpoint = ModelCheckpoint(model_file_name, monitor='val_loss', verbose=1, save_best_only=True, mode='min')
        # fit model
        history = model.fit(X_train, y_train_flat, epochs=num_epochs, batch_size=batch_size, verbose=1, validation_data=(X_val, y_val_flat),
                            callbacks=[reduce_lr,earlyStop,checkpoint])
        np.save(history_file_name, history.history)
    return model, history

"""
    MODELS TRAINING
"""
def write_logs(modelName):
    end_time = datetime.now()
    duration = end_time - begin_time
    m,s = divmod(duration.total_seconds(), 60)
    f = open('logs.txt', "a")
    f.write(modelName + '\n')
    f.write('Begin datetime: {} \t\t End datetime: {}'.format(begin_time, end_time) + '\n')
    f.write('Duration H:M:S is {}:{}:{}'.format(m//60,m%60,s) + '\n')
    model.summary(print_fn=lambda log: f.write(log + '\n'))
    f.write(model.to_json() + '\n\n')
    f.close()

# training by cell size
for lstm_layer_size in lstm_cells:
    for i in range(model_repeat_time):
        begin_time = datetime.now()
        model, history = lstm_autoencoder_flat(True, i) # training Autoencoder models
        # model, history = basic_gru_flat(True, i)  # training GRU models
        # model, history = basic_bidirectional_lstm_flat(True, i)  # training LSTM Bidirectional models
        write_logs('LSTM-AE_cells-{}_epochs-{}_batch-{}_Test-{}'.format(lstm_layer_size, num_epochs, batch_size, i))

"""
    ADDITIONAL METHODS
"""
def transform_and_scale(x):
    i = 7 - x.shape[2]
    if(i != 0):
        x = np.append(x, np.zeros((x.shape[0], x.shape[1], i)), axis=2)
    i = x.shape[1]
    x = x.reshape(x.shape[0]*x.shape[1], x.shape[2])
    x = scaler.inverse_transform(x)
    x = x.reshape(x.shape[0]//i, i, x.shape[1])
    return x

def transform_and_scale_utm(x):
    i = x.shape[1]
    x = x.reshape(x.shape[0]*x.shape[1], x.shape[2])
    x = scaler_utm.inverse_transform(x)
    x = x.reshape(x.shape[0]//i, i, x.shape[1])
    return x

def inverse_pol_to_coordinates(x, y_pred):
    x = transform_and_scale(x)
    lat = x[:,seq_input_length-1,0]
    long = x[:,seq_input_length-1,1]

    prediction = list()
    for j in range(len(x)):
        temp_lat = lat[j]
        temp_long = long[j]
        for i in range(seq_output_length):
            endLon,endLat,backAzimuth = geod.fwd(temp_long,temp_lat,y_pred[j,i,0],y_pred[j,i,1])
            prediction.append([endLat, endLon])
            temp_lat = endLat
            temp_long = endLon
    prediction = np.array(prediction)
    prediction = prediction.reshape(len(prediction)//seq_output_length,seq_output_length,2)
    return prediction

def inverse_utm_to_coordinates(X_val, y_pred):
    X_val = transform_and_scale(X_val)
    lat = X_val[:,seq_input_length-1,0]
    long = X_val[:,seq_input_length-1,1]

    y_pred = y_pred.reshape(len(y_pred)*seq_output_length, 2)
    y_pred = scaler_utm.inverse_transform(y_pred)
    y_pred = y_pred.reshape(len(y_pred)//seq_output_length,seq_output_length, 2)

    prediction = list()

    for j in range(len(X_val)):
        temp_lat = lat[j]
        temp_long = long[j]
        for i in range(seq_output_length):
            utm_x, utm_y = p(temp_lat, temp_long, inverse=False)
            temp_x = utm_x + y_pred[j,i,0]
            temp_y = utm_y + y_pred[j,i,1]

            endLat, endLon = p(temp_x, temp_y, inverse=True)
            prediction.append([endLat, endLon])
            temp_lat = endLat
            temp_long = endLon
    prediction = np.array(prediction)
    prediction = prediction.reshape(len(prediction)//seq_output_length,seq_output_length,2)

    return prediction

def plt_dynamic_training(x, vy, ty, ax, title, colors=['b']):
    ax.plot(x, vy, 'r', label="Validation Loss")
    ax.plot(x, ty, 'b', label="Train Loss")
    ax.legend()
    plt.title(title)
    plt.grid()
    plt.savefig(f'image_{lstm_layer_size}_{i}.jpg')
    fig.canvas.draw()

"""
    MODELS EVALUATION
"""
scaler = pickle.load(open("scaler.pkl", "rb"))
scaler_pol = pickle.load(open("scaler_pol.pkl", "rb"))
scaler_utm = pickle.load(open("scaler_utm.pkl", "rb"))

ytrue = np.load("y_test.npy")
ytrue = transform_and_scale(ytrue)

loss_avg = list()
distance_loss = np.zeros((model_repeat_time, len(lstm_cells)))
ii = 0
for lstm_layer_size in lstm_cells:
    score = 0
    for i in range(model_repeat_time):
        fig, ax = plt.subplots(1,1)
        ax.set_xlabel('Epoch')
        ax.set_ylabel('MAEH (Loss)')
        # load model and history without training new model
        model, history = lstm_autoencoder_flat(False, i)
        # list of each epoch numbers
        x = list(range(1, len(history['val_loss'])+1))
        # model evaluation (forecast/prediction)
        yhat = model.predict(X_test, verbose=0, batch_size=batch_size)

        """FOR LAT LON"""
        # EXP-A and EXP-B
        yhat = yhat.reshape(len(yhat), seq_output_length, 2)
        yhat = transform_and_scale(yhat)

        """FOR POL"""
        # EXP-C
        """yhat = yhat.reshape(len(yhat), seq_output_length, 2) #flat
        yhat = yhat.reshape(len(yhat)*seq_output_length, 2)
        yhat = scaler_pol.inverse_transform(yhat)
        yhat = yhat.reshape(len(yhat)//seq_output_length,seq_output_length, 2)
        yhat = inverse_pol_to_coordinates(X_test, yhat)"""

        """FOR UTM"""
        # EXP-D
        """yhat = inverse_utm_to_coordinates(X_test,yhat)"""

        error, dist = ext.true_and_predicted_distance(ext, ytrue, yhat, 'k')
        score += error
        distance_loss[i,ii] = error

        # training plot
        vy = history['val_loss']
        ty = history['loss']
        plt_dynamic_training(x[3:], vy[3:], ty[3:], ax, f'Cross Validation - AE-LSTM-Exp_A_cells-{lstm_layer_size}_Experiment-{i}')
    ii = ii + 1
    loss_avg.append(score/model_repeat_time)
    print('AVERAGE TEST SCORE: ', score/model_repeat_time)


text:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import csv
from plot import to_float
from plot import to_label

"""
    some queries may fail, so it worth specical treatment.
"""

bar_width = 0.6
xstick_offset = 0.4

font = {
    'family': 'serif',
    'size': 21}

matplotlib.rc('font', **font)
dpi = plt.figure().dpi
matplotlib.rcParams.update({'figure.autolayout': True})
# http://www.mulinblog.com/a-color-palette-optimized-for-data-visualization/
colors = ['#5da5da', '#faa43a', '#60bd68', '#f17cb0', '#b2912f', '#b276b2']
path = "/Users/chushumo/Project/papers/2015-multiwayjoin-revision/images"
query = "fb_q5"

def plot_runtime(algebras, data, output_name):
    # plot wall clock time
    runtime = []    # query runtime
    std = []        # standard deviation
    for row in data:
        if row[0] == query:
            runtime.append(to_float(row[2]))
            std.append(to_float(row[3]))

    ind = np.arange(len(runtime))  # the x locations for the groups
    fig, ax = plt.subplots()
    rects = ax.bar(ind+0.1, runtime, width=bar_width, color=colors,  yerr=std)
    ax.set_ylabel('Time (sec)')
    ax.set_xticks(ind+xstick_offset)
    ax.set_xticklabels(algebras)
    if query == "fb_q5":
        ax.set_ylim((0, 14000))
    # label the bars
    for rect in rects:
        height = rect.get_height()
        if int(height) == 0:
            ax.text(
                rect.get_x()+rect.get_width()/2.,
                1.05*height, "FAIL",
                ha='center', va='bottom', color='red')
        else:
            ax.text(
                rect.get_x()+rect.get_width()/2.,
                1.05*height, to_label(height),
                ha='center', va='bottom')
    print "outputing {}".format(output_name)
    plt.savefig(output_name, format='pdf', dpi=dpi)

def plot_cpu(algebras, data, output_name):
    # plot wall clock time
    cputime = []    # query cputime
    for row in data:
        if row[0] == query:
            cputime.append(to_float(row[8]))

    ind = np.arange(len(cputime))  # the x locations for the groups
    fig, ax = plt.subplots()
    rects = ax.bar(ind+0.1, cputime, width=bar_width, color=colors)
    ax.set_ylabel('CPU Time (sec)')
    ax.set_xticks(ind+xstick_offset)
    ax.set_xticklabels(algebras)
    if query == "fb_q5":
        ax.set_ylim((0, 280000))
    # label the bars
    for rect in rects:
        height = rect.get_height()
        if int(height) == 0:
            ax.text(
                rect.get_x()+rect.get_width()/2.,
                1.05*height, "FAIL",
                ha='center', va='bottom', color='red')
        else:
            ax.text(
                rect.get_x()+rect.get_width()/2.,
                1.05*height, to_label(height),
                ha='center', va='bottom')
    print "outputing {}".format(output_name)
    plt.savefig(output_name, format='pdf', dpi=dpi)

def plot_shuffle(algebras, data, output_name):
    shuffle_size = []
    for row in data:
        if row[0] == query:
            shuffle_size.append(to_float(row[10]))
    ind = np.arange(len(shuffle_size))  # the x locations for the groups
    fig, ax = plt.subplots()
    rects = ax.bar(ind+0.1, shuffle_size, width=bar_width, color=colors)
    ax.set_ylabel('Tuples shuffled (million)')
    # ax.set_xlabel('Physical Algebra')
    ax.set_xticks(ind+xstick_offset)
    ax.set_xticklabels(algebras)
    if query == "fb_q5":
        ax.set_ylim((0, 17000))
    # label the bar
    for rect in rects:
        height = rect.get_height()
        if int(height) == 0:
            ax.text(
                rect.get_x()+rect.get_width()/2.,
                1.05*height, "FAIL",
                ha='center', va='bottom', color='red')
        else:
            ax.text(
                rect.get_x()+rect.get_width()/2.,
                1.05*height, to_label(height),
                ha='center', va='bottom')
    print "outputing {}".format(output_name)
    plt.savefig(output_name, format='pdf', dpi=dpi)

def plot():
    fname = "csvs/SIGMOD Experiment - summary.csv"
    agbrs = ('RS_HJ', 'RS_TJ', 'BR_HJ', 'BR_TJ', 'HC_HJ', 'HC_TJ')
    with open(fname, "rU") as f:
        csvreader = csv.reader(f)
        data = [r for r in csvreader]
    plot_runtime(agbrs, data, "{}/{}_wall_time.pdf".format(
        path, query))
    plot_cpu(agbrs, data, "{}/{}_cpu_time.pdf".format(
        path, query))
    plot_shuffle(agbrs, data, "{}/{}_shuffle_size.pdf".format(
        path, query))

if __name__ == '__main__':
    plot()

max_stars_repo_name: stechu/myria-commandline

max_issues_repo_name: stechu/myria-commandline

max_forks_repo_name: stechu/myria-commandline


text:
import os
from datetime import datetime
import pprint
import glob
import logging
import pickle
from pprint import pformat
import sys
import yaml
import pathpy
import pandas as pd
import dotenv
import numpy as np
from matplotlib import pyplot as plt
from scipy import stats

import sacred

from src.data_processing import generate_temporal_network, get_runs
from src.preprocess_experiment import create_train_test_split
from src.utils import config_adapt

def analyze(config):
    def report(i):
        i % 10 == 0 and print(f"Processed {i} logs of {total}")

    logger = logging.getLogger("hids.analyze")

    train, _ = create_train_test_split(config["data"]["runs"], config["model"]["train_examples"])
    runs = get_runs(config["data"]["runs"], train)

    total = len(runs)

    normal_graphs = [
        report(i) or generate_temporal_network(scenario, None)
        for i, scenario in enumerate(runs["path"])
    ]

    save_path = os.path.join(config["analyze"]["figures"], f'{config["dataset"]}_time_analyze.png')

    inter_event_times_all = analyze_time(normal_graphs, save_path)

    statistics = stats.describe(inter_event_times_all)

    logger.debug("runs for training")
    logger.debug(runs)

    logger.info("Summary time statistics")
    logger.info("")
    logger.info(statistics)

def analyze_time(temporal_nets, save_path=None):

    multi = [x.inter_event_times().flatten() for x in temporal_nets]
    means = [x.mean() / 1000.0 for x in multi]
    max_times = [x.max() / 1000000.0 for x in multi]
    multi = np.concatenate(multi)

    observation_lengths = [x.observation_length() for x in temporal_nets]

    plt.tight_layout()
    # these are my measurements, unsorted
    fig = plt.figure(figsize=(8, 20))
    fig.subplots_adjust(hspace=0.3)

    plt.subplot(411)
    plt.hist(
        multi,
        bins=np.logspace(np.log10(1), np.log10(np.percentile(multi, 99)), 50),
        weights=100 * np.ones(len(multi)) / len(multi),
        color="blue",
        edgecolor="black",
        alpha=0.5,
        range=(0, np.percentile(multi, 99)),
    )
    plt.gca().set_xscale("log")
    plt.title("Histogram Time Differences (Log)")
    plt.xlabel("Time differences [$\mu s$]")
    plt.ylabel("Percentage [%]")

    plt.subplot(412)
    plt.hist(
        multi,
        bins=50,
        weights=100 * np.ones(len(multi)) / len(multi),
        color="blue",
        edgecolor="black",
        alpha=0.5,
        range=(0, np.percentile(multi, 90)),
    )
    plt.title("Histogram Time Differences of 90 percentile")
    plt.xlabel("Time differences [$\mu s$]")
    plt.ylabel("Percentage [%]")

    plt.subplot(413)
    plt.hist(
        means,
        bins=50,
        weights=100 * np.ones(len(means)) / len(means),
        color="blue",
        edgecolor="black",
        alpha=0.5,
        range=(0, np.percentile(means, 99)),
    )
    plt.title("Histogram Mean Time Difference")
    plt.xlabel("Mean time differences [ms]")
    plt.ylabel("Percentage [%]")

    plt.subplot(414)
    plt.hist(
        max_times,
        bins=50,
        weights=100 * np.ones(len(means)) / len(means),
        color="blue",
        edgecolor="black",
        alpha=0.5,
        range=(0, np.percentile(max_times, 99)),
    )
    plt.title("Histogram Maximum Time Difference")
    plt.xlabel("Maximum time differences [s]")
    plt.ylabel("Percentage [%]")

    if save_path:
        os.makedirs(os.path.dirname(save_path), exist_ok=True)
        plt.savefig(save_path, bbox_inches="tight")

    return multi

reason: import numpy,from scipy


text:
import zmq
import signal
import sys
from time import sleep
from Logger import Logger
import numpy as np

class gpsPoint:
    def __init__(self, x, y, e=10):
        self.x = x
        self.y = y
        self.Lat = x*np.pi/180.0
        self.Lon = y*np.pi/180.0
        self.Err = e

class gpsNavigation:
    def __init__(self,lenOfFilter):
        self.lengthOfFilter = lenOfFilter
        self.ItemKey = {"ID":0,"TIME":1,"LON":2,"LAT":3,"P_ERR":4,"DIR":5}
        self.IDKey = {"GPS1":0, "GPS2":1}

        p0 = gpsPoint(0,0)
        self.gpsPoints = [[],[]]
        for i in range(lenOfFilter):
            self.gpsPoints[0].append(p0)
            self.gpsPoints[1].append(p0)

        self.gpsWayPoint = gpsPoint(0,0)

        self.dir = np.zeros((2,self.lengthOfFilter),np.float)
        self.Time = [0,0]

    def shift(self,sensorID):
        for j in range(self.lengthOfFilter-2,-1,-1):
            self.gpsPoints[sensorID][j+1] = self.gpsPoints[sensorID][j]
            self.dir[sensorID,j+1] = self.dir[sensorID,j]

        self.gpsPoints[sensorID][0] = gpsPoint(0,0)
        self.dir[sensorID,0] = 0

    def update(self,sensorID,datastring):
        itemid,itemVal = datastring.split(":")
        if(itemid == "LON") or (itemid == "LAT"):
            pos,sign = itemVal.split(",")
            if(sign == "N") or (sign == "E"):
                ## Position on globe is from -180° to +180°
                itemVal = float(pos)
            else:
                itemVal = -float(pos)
        if(itemid == "LON"):
            self.gpsPoints[sensorID][0].Lon = itemVal
            return 0
        elif(itemid == "LAT"):
            self.gpsPoints[sensorID][0].Lat = itemVal
            return 0
        elif(itemid == "P_ERR"):
            self.gpsPoints[sensorID][0].Err = float(itemVal)
            return 0
        elif(itemid == "DIR"):
            self.dir[sensorID,0] = int(itemVal)
            return 0
        elif(itemid == "TIME"):
            self.Time[sensorID] = float(itemVal)
            return 0
        elif(itemid == "ID"):
            return 0
        else:
            return -1

    def getDataString(self,sensorID):
        data  = "SensorID:" + str(sensorID)
        data += ";LON:" + str(self.gpsPoints[sensorID][0].Lon)
        data += ";LAT:" + str(self.gpsPoints[sensorID][0].Lat)
        data += ";P_ERR:" + str(self.gpsPoints[sensorID][0].PErr)
        data += ";DIR:" + str(self.dir[sensorID,0])
        data += ";TIME:" + str(self.Time[sensorID])
        return data

    def consistent(self):
        return (self.Time[0] == self.Time[1]) and (self.Time[0] > 1)

    def getCorrectedPosition(self):
        # get realdistance from points
        #   if(distance < dilutions:
        #       compute new position with higher precision
        #   else:
        #       compute new position with lower precision
        corX = (self.gpsPoints[0][0].x + self.gpsPoints[1][0].x) / 2.0
        corY = (self.gpsPoints[0][0].y + self.gpsPoints[1][0].y) / 2.0
        corErr = (self.gpsPoints[0][0].Err + self.gpsPoints[1][0].Err) / 2.0
        curPos = gpsPoint(corX,corY,corErr)
        return curPos

    def getDirAndDist(self,A,B):
        # You are at point A and want to go to point B
        dLat = B.Lat - A.Lat
        dLon = B.Lon - A.Lon

        R = 6378137.0

        distance = None
        angle = None

        ## Based on Orthrodromic navigation - see wiki
        alf = 2*np.arcsin(np.sqrt((np.sin((dLat)/2)**2
              +np.cos(B.Lat)*np.cos(A.Lat)*(dLon/2)**2)))

        angle = (np.arcsin((np.cos(B.Lat)*np.sin(dLon))/np.sin(alf))
                 /np.pi*180)%360

        if(dLat < 0):
            if(dLon < 0):
                angle -= 90
            elif(dLon>0):
                angle += 90
            else:
                angle += 180

        distance = np.round(alf*R)
        angle = np.round(angle)

        return distance,angle

class gpsModule:
    def __init__(self):

        self.logger = Logger("gpsNavigation")
        self.InPorts = ["10101","10102","10110"]
        self.OutPorts = ["10105"]
        self.enabled = True
        signal.signal(signal.SIGINT, self.sigINT_Handler)
        zMQC = zmq.Context()
        self.subscriber = zMQC.socket(zmq.SUB)
        self.publisher = zMQC.socket(zmq.PUB)

        self.GPSData = gpsNavigation(5)
        self.WayPoint = gpsPoint(0,0,0)
        self.waypointSet = False

        for port in self.InPorts:
            self.subscriber.connect('tcp://127.0.0.1:'+port)
            self.logger.save_line("Binded to port: " + port)
        self.subscriber.setsockopt(zmq.SUBSCRIBE, b"")

        for port in self.OutPorts:
            self.publisher.bind('tcp://127.0.0.1:'+port)
            self.logger.save_line("Connected to port: " + port)

        sleep(0.5)

    def sigINT_Handler(self, signal, frame):
        print("\ngpsNavigation detected SigINT signal")
        self.logger.save_line("Signal SigINT detected")
        self.enabled = False
        for port in self.InPorts:
            self.subscriber.disconnect('tcp://127.0.0.1:'+port)
            self.logger.save_line("SUB Disconnected from local port: " + port)
        for port in self.OutPorts:
            self.publisher.disconnect('tcp://127.0.0.1:'+port)
            self.logger.save_line("PUB Disconnected from local port: " + port)
        sys.exit(0)

    def updateData(self,line):
        lineSplit = line.split(";")
        gpsid = lineSplit[0].split(":")[1]
        if(gpsid in self.GPSData.IDKey):
            gpid = self.GPSData.IDKey[gpsid]
            self.GPSData.shift(gpid)
            for item in lineSplit:
                if(self.GPSData.update(gpid,item) == -1):
                   self.logger.save_line("Unsupported sentence : <"
                                         + item + ">")

    def updateWayPoint(self,line):
        lineSplit = line.split(";")
        self.WayPoint = gpsPoint(
                          float(lineSplit[1].split(":")[1])
                        , float(lineSplit[2].split(":")[1])
                        , float(lineSplit[3].split(":")[1]) )
        self.waypointSet = True

    def run(self):
        while(self.enabled):
            message = self.subscriber.recv_string()
            if(message.find("ID:GPS")>=0):
                self.updateData(message)
                if(self.GPSData.consistent()):
                    # watchdog = 0
                    if(self.waypointSet):
                        A = self.GPSData.getCorrectedPosition()
                        dist,azim = self.GPSData.getDirAndDist(A,self.WayPoint)
                        self.publisher.send_string("ID:NAV;DIST:"+str(dist)
                                                   +";AZIM:"+str(azim))
                    else:
                        A = self.GPSData.getCorrectedPosition()
                        self.logger.save_line("New position: LON: "
                                              + str(A.Lon) + " LAT: "
                                              + str(A.Lat))
            elif(message.find("ID:MAIN")>=0):
                self.updateWayPoint(message)
            else:
                self.logger.save_line("Ignoring sentense: <" + message + ">")
                # watchdog ++
                # if watchdog >= FiltLen:
                #   GPSdata.useOneSourceOnly == True
                # else:
                #   pass (wait for next update)
                sleep(0.0001)

max_stars_repo_name: TeaPackCZ/RobotZed

max_issues_repo_name: TeaPackCZ/RobotZed

max_forks_repo_name: TeaPackCZ/RobotZed


text:
import plotly.plotly as py
from plotly.graph_objs import *
import plotly
import numpy as np
import pandas as pd
import os.path
import plotly.graph_objs as go

path = "/ros_ws/src/data"

def plot_box_plot(data, title, y_key, image_file):
	keys = data.keys()
	traces = []
	for k in keys:
		traces.append(Box(y=data[k][y_key], name=k, boxpoints='all', jitter=0.3, pointpos=-1.8))
	layout = go.Layout(
	    title=title,
	    yaxis=dict(
	        title='Time in milliseconds',
	        titlefont=dict(
	            family='Courier New, monospace',
	            size=18,
	            color='#7f7f7f'
	        )
	    )
	)
	fig=dict(data=traces, layout=layout)
	plotly.offline.plot(
    fig,
    filename= './'+image_file+'.html',
    image='png')

def plot_volume_exploration_vs_time(data, max_time):
	keys = data.keys()
	traces = []
	x_side = 12 - (-15)
	y_side = 20 - (-20)
	z_side = 14 - 1
	max_volume_m = x_side*y_side*z_side
	trace_max_volume = go.Scatter(
	    x=[0, max_time],
	    y=[max_volume_m, max_volume_m],
	    name='Geofence volume',
	    mode='lines',
	)
	traces.append(trace_max_volume)
	for k in keys:
		traces.append(Scatter(x = data[k]['time ellapsed millis'], y=data[k]['volume'], name=k, mode = 'lines'))
	layout = go.Layout(
	    title='Progression of volume of explored space cases analyzed',
	    xaxis=dict(
	        title='Time in milli seconds',
	        titlefont=dict(
	            family='Courier New, monospace',
	            size=18,
	            color='#7f7f7f'
	        )
	    ),
	    yaxis=dict(
	        title='Explored volume in cubic meters',
	        titlefont=dict(
	            family='Courier New, monospace',
	            size=18,
	            color='#7f7f7f'
	        )
	    )
	)
	fig=dict(data=traces, layout=layout)
	plotly.offline.plot(fig, filename='./volume_exploration_vs_time.html', image='png')

def plot_ram_vs_time (data):
	keys = data.keys()
	traces = []
	for k in keys:
		traces.append(Scatter(y=data[k]['RAM'], name=k, mode = 'lines'))

	layout = go.Layout(
	    title='Progression of RAM used',
	    xaxis=dict(
	        title='Time in seconds',
	        titlefont=dict(
	            family='Courier New, monospace',
	            size=18,
	            color='#7f7f7f'
	        )
	    ),
	    yaxis=dict(
	        title='RAM used in Bytes',
	        exponentformat = 'none',
	        titlefont=dict(
	            family='Courier New, monospace',
	            size=18,
	            color='#7f7f7f'
	        )
	    )
	)
	fig=dict(data=traces, layout=layout)
	plotly.offline.plot(fig, filename='./ram_vs_time.html', image='png')

def extract_ram_data():
	csv_filepath = "./../data/"+str(file)+"/ram.csv"
	ram_data = pd.read_csv(csv_filepath, delim_whitespace=True)
	# print ram_data
	# print ram_data['free']
	data = {}
	data['RAM'] = ram_data['used']
	return data

def extract_volume_data(csv_filepath):
	volume_data = pd.read_csv(csv_filepath)
	data = {}
	data['time ellapsed millis'] = volume_data['time ellapsed millis'];
	data['volume'] = volume_data['volume'];
	return data

def extract_frontier_data(csv_filepath):
    frontiers_computation_time = pd.read_csv(csv_filepath)
    data = {}
    data['computation_time_millis'] = frontiers_computation_time['computation_time_millis']
    return data

def extract_lazy_theta_star_data(csv_filepath):
    frontiers_computation_time = pd.read_csv(csv_filepath)
    data = {}
    data['computation_time_millis'] = frontiers_computation_time['computation_time_millis']
    return data

def run_computation_time(csv_filepath):
	data = extract_computation_time_path(csv_filepath)
	plot_box_plot(data, "frontiers_computation_time_millis_box")

# exploration_time.csv

frontier_computation_time_data = {}
lazyThetaStar_computation_data = {}
volume_data = {}
ram_data = {}
max_time = 0
for file in os.walk('./../data').next()[1]:
	frontier_filename = "./../data/"+str(file)+"/current/frontiers_computation_time.csv"
	frontier_computation_time_data[file] = extract_frontier_data( frontier_filename )
	lazyThetaStar_filename = "./../data/"+str(file)+"/lazyThetaStar_computation_time.csv"
	lazyThetaStar_computation_data[file] = extract_lazy_theta_star_data( lazyThetaStar_filename )
	volume_filename = "./../data/"+str(file)+"/volume_explored.csv"
	volume_data[file] = extract_volume_data( volume_filename )
	ram_data[file]    = extract_ram_data()
	max_time = max(volume_data[file]['time ellapsed millis'])
plot_box_plot(frontier_computation_time_data, 'Amount of time used for generating frontiers', 'computation_time_millis', 'frontiers_time_box')
plot_box_plot(lazyThetaStar_computation_data, 'Amount of time used for path planning with Lazy Theta Star', 'computation_time_millis', 'lazyThetaStar_time_box')
plot_ram_vs_time(ram_data)
plot_volume_exploration_vs_time(volume_data, max_time)

# variables= {}
# execfile( "compare_box_graphs.py", variables )


text:
import game_map,player,sys,render
import numpy as np
import pygame as pyg
import time as t
from collections import namedtuple
from minimap import Minimap
from primitives import *

pyg.init()
internal_resolution = Resolution(256, 150)
screen_w,screen_h = resolution = Resolution(1024, 600)

pyg.display.set_caption("Raycast Engine v0.01")
screen = pyg.display.set_mode(resolution)

#velocity
v = 0.1
m = game_map.Map('map.txt')

p = player.Player(8, 4, v, 0, (np.pi/2))

mini = Minimap(screen, m, resolution.scale(0.35), Point(0,0))
while True:

    startTime = t.perf_counter()
    #key detect stuff#
    for event in pyg.event.get():
        if event.type == pyg.QUIT:
            sys.exit()

    #key detect#
    pressed = pyg.key.get_pressed()
    prevpos = p.x,p.y
    p.handle_keys(pressed)

    #check collition in a ghetto way#
    try:
        if (n:=m.map_array[int(p.y)][int(p.x)]) != 0:
            p.updatePosition(*prevpos)
    except IndexError:
        p.updatePosition(prevpos[0],prevpos[1])

    #draw bg#
    screen.fill((0,0,0))
    pyg.draw.rect(screen,(63,63,63),(0,int(screen_h/2),screen_w,int(screen_h/2)),0)

    coords = render.render(screen,p,resolution,internal_resolution,m.map())
    mini.draw()

    pyg.display.flip()

    t.sleep(max(0,1/60 - t.perf_counter() + startTime))