text:
# Copyright 2018 Regents of the University of Colorado. All Rights Reserved.
# Released under the MIT license.
# This software was developed at the University of Colorado's Laboratory for Atmospheric and Space Physics.
# Verify current version before use at: https://github.com/MAVENSDC/PyTplot

from __future__ import division
import numpy as np
from bokeh.plotting.figure import Figure
from bokeh.models import (ColumnDataSource, DatetimeAxis, HoverTool,
                          Range1d, Span, Title, Legend, BoxAnnotation)
from bokeh.models.glyphs import Line
from bokeh.models.tools import BoxZoomTool

import pytplot
from pytplot import tplot_utilities
from bokeh.models.formatters import DatetimeTickFormatter
from bokeh.models.markers import X

dttf = DatetimeTickFormatter(microseconds=["%H:%M:%S"],
                             milliseconds=["%H:%M:%S"],
                             seconds=["%H: %M:%S"],
                             minsec=["%H:%M:%S"],
                             minutes=["%H:%M:%S"],
                             hourmin=["%H:%M:%S"],
                             hours=["%H:%M"],
                             days=["%F"],
                             months=["%F"],
                             years=["%F"])

class TVarFigure1D(object):

    def __init__(self, tvar_name, auto_color, show_xaxis=False, slice=False):
        self.tvar_name = tvar_name
        self.auto_color = auto_color
        self.show_xaxis = show_xaxis
        if 'show_all_axes' in pytplot.tplot_opt_glob:
            if pytplot.tplot_opt_glob['show_all_axes']:
                self.show_xaxis = True
        self.slice = slice

        # Variables needed across functions
        self.colors = ['black', 'red', 'green', 'navy', 'orange', 'firebrick', 'pink', 'blue', 'olive']
        self.lineglyphs = []
        self.linenum = 0
        self.interactive_plot = None

        self.fig = Figure(x_axis_type='datetime',
                          tools=pytplot.tplot_opt_glob['tools'],
                          y_axis_type=self._getyaxistype())
        self.fig.add_tools(BoxZoomTool(dimensions='width'))
        self._format()

    @staticmethod
    def get_axis_label_color():
        if pytplot.tplot_opt_glob['black_background']:
            text_color = '#000000'
        else:
            text_color = '#FFFFFF'
        return text_color

    @staticmethod
    def getaxistype():
        axis_type = 'time'
        link_y_axis = False
        return axis_type, link_y_axis

    def getfig(self):
        if self.slice:
            return [self.fig, self.interactive_plot]
        else:
            return [self.fig]

    def setsize(self, width, height):
        self.fig.plot_width = width
        if self.show_xaxis:
            self.fig.plot_height = height + 22
        else:
            self.fig.plot_height = height

    def add_title(self):
        if 'title_text' in pytplot.tplot_opt_glob:
            if pytplot.tplot_opt_glob['title_text'] != '':
                title1 = Title(text=pytplot.tplot_opt_glob['title_text'],
                               align=pytplot.tplot_opt_glob['title_align'],
                               text_font_size=pytplot.tplot_opt_glob['title_size'],
                               text_color=self.get_axis_label_color())
                self.fig.title = title1
                self.fig.plot_height += 22

    def buildfigure(self):
        self._setminborder()
        self._setxrange()
        self._setxaxis()
        self._setyrange()
        self._addtimebars()
        self._visdata()
        self._setxaxislabel()
        self._setyaxislabel()
        self._addhoverlines()
        self._addlegend()

    def _format(self):
        # Formatting stuff
        self.fig.grid.grid_line_color = None
        self.fig.axis.major_tick_line_color = None
        self.fig.axis.major_label_standoff = 0
        self.fig.xaxis.formatter = dttf
        self.fig.title = None
        self.fig.toolbar.active_drag = 'auto'
        if not self.show_xaxis:
            self.fig.xaxis.major_label_text_font_size = '0pt'
            self.fig.xaxis.visible = False

    def _setxrange(self):
        # Check if x range is not set, if not, set good ones
        if 'x_range' not in pytplot.tplot_opt_glob:
            datasets = [pytplot.data_quants[self.tvar_name]]
            x_min_list = []
            x_max_list = []
            for oplot_name in pytplot.data_quants[self.tvar_name].attrs['plot_options']['overplots']:
                datasets.append(pytplot.data_quants[oplot_name])
            for dataset in datasets:
                x_min_list.append(np.nanmin(dataset.coords['time']))
                x_max_list.append(np.nanmax(dataset.coords['time']))
            pytplot.tplot_opt_glob['x_range'] = [np.nanmin(x_min_list), np.nanmax(x_max_list)]
            tplot_x_range = [np.nanmin(x_min_list), np.nanmax(x_max_list)]
            if self.show_xaxis:
                pytplot.lim_info['xfull'] = tplot_x_range
                pytplot.lim_info['xlast'] = tplot_x_range

        # Bokeh uses milliseconds since epoch for some reason
        x_range = Range1d(int(pytplot.tplot_opt_glob['x_range'][0]) * 1000.0,
                          int(pytplot.tplot_opt_glob['x_range'][1]) * 1000.0)
        self.fig.x_range = x_range

    def _setyrange(self):
        if self._getyaxistype() == 'log':
            if pytplot.data_quants[self.tvar_name].attrs['plot_options']['yaxis_opt']['y_range'][0] <= 0 or \
                    pytplot.data_quants[self.tvar_name].attrs['plot_options']['yaxis_opt']['y_range'][1] <= 0:
                return
        y_range = Range1d(pytplot.data_quants[self.tvar_name].attrs['plot_options']['yaxis_opt']['y_range'][0],
                          pytplot.data_quants[self.tvar_name].attrs['plot_options']['yaxis_opt']['y_range'][1])
        self.fig.y_range = y_range

    def _setminborder(self):
        self.fig.min_border_bottom = pytplot.tplot_opt_glob['min_border_bottom']
        self.fig.min_border_top = pytplot.tplot_opt_glob['min_border_top']
        if 'vertical_spacing' in pytplot.tplot_opt_glob:
            self.fig.min_border_bottom = int(pytplot.tplot_opt_glob['vertical_spacing'] / 2.0)
            self.fig.min_border_top = int(pytplot.tplot_opt_glob['vertical_spacing'] / 2.0)

    def _addtimebars(self):
        for time_bar in pytplot.data_quants[self.tvar_name].attrs['plot_options']['time_bar']:
            time_bar_line = Span(location=time_bar['location']*1000.0, dimension=time_bar['dimension'],
                                 line_color=time_bar['line_color'], line_width=time_bar['line_width'])
            self.fig.renderers.extend([time_bar_line])

    def _set_roi_lines(self, dataset):
        # Locating the two times between which there's a roi
        time = dataset.coords['time'].values
        roi_1 = pytplot.tplot_utilities.str_to_int(pytplot.tplot_opt_glob['roi_lines'][0])
        roi_2 = pytplot.tplot_utilities.str_to_int(pytplot.tplot_opt_glob['roi_lines'][1])
        # find closest time to user-requested time
        x = np.asarray(time)
        x_sub_1 = abs(x - roi_1 * np.ones(len(x)))
        x_sub_2 = abs(x - roi_2 * np.ones(len(x)))
        x_argmin_1 = np.nanargmin(x_sub_1)
        x_argmin_2 = np.nanargmin(x_sub_2)
        x_closest_1 = x[x_argmin_1]
        x_closest_2 = x[x_argmin_2]
        # Create roi box
        roi_box = BoxAnnotation(left=x_closest_1*1000.0, right=x_closest_2*1000.0, fill_alpha=0.2, fill_color='grey',
                                line_color='red', line_width=2.5)
        self.fig.renderers.extend([roi_box])

    def _setxaxis(self):
        xaxis1 = DatetimeAxis(major_label_text_font_size='0pt', formatter=dttf)
        xaxis1.visible = False
        self.fig.add_layout(xaxis1, 'above')

    def _getyaxistype(self):
        if 'y_axis_type' in pytplot.data_quants[self.tvar_name].attrs['plot_options']['yaxis_opt']:
            return pytplot.data_quants[self.tvar_name].attrs['plot_options']['yaxis_opt']['y_axis_type']
        else:
            return 'linear'

    def _setcolors(self):
        if 'line_color' in pytplot.data_quants[self.tvar_name].attrs['plot_options']['extras']:
            self.colors = pytplot.data_quants[self.tvar_name].attrs['plot_options']['extras']['line_color']

    def _setxaxislabel(self):
        self.fig.xaxis.axis_label = pytplot.data_quants[self.tvar_name].attrs['plot_options']['xaxis_opt']['axis_label']
        self.fig.xaxis.axis_label_text_font_size = str(pytplot.data_quants[self.tvar_name].attrs['plot_options']['extras']['char_size'])+'pt'
        self.fig.xaxis.axis_label_text_color = self.get_axis_label_color()

    def _setyaxislabel(self):
        self.fig.yaxis.axis_label = pytplot.data_quants[self.tvar_name].attrs['plot_options']['yaxis_opt']['axis_label']
        self.fig.yaxis.axis_label_text_font_size = str(pytplot.data_quants[self.tvar_name].attrs['plot_options']['extras']['char_size'])+'pt'
        self.fig.yaxis.axis_label_text_color = self.get_axis_label_color()

    def _visdata(self):
        self._setcolors()

        datasets = [pytplot.data_quants[self.tvar_name]]
        for oplot_name in pytplot.data_quants[self.tvar_name].attrs['plot_options']['overplots']:
            datasets.append(pytplot.data_quants[oplot_name])

        for dataset in datasets:
            # Get Linestyle
            line_style = None
            if 'linestyle' in pytplot.data_quants[self.tvar_name].attrs['plot_options']['extras']:
                line_style = pytplot.data_quants[self.tvar_name].attrs['plot_options']['extras']['linestyle']

            # Get a list of formatted times
            corrected_time = []
            for x in dataset.coords['time'].values:
                corrected_time.append(tplot_utilities.int_to_str(x))

            # Bokeh uses milliseconds since epoch for some reason
            x = dataset.coords['time'].values * 1000.0

            # Add region of interest (roi) lines if applicable
            if 'roi_lines' in pytplot.tplot_opt_glob.keys():
                self._set_roi_lines(dataset)

            plot_options = dataset.attrs['plot_options']
            df = pytplot.tplot_utilities.convert_tplotxarray_to_pandas_dataframe(dataset.name, no_spec_bins=True)

            # Create lines from each column in the dataframe
            for column_name in df.columns:
                y = df[column_name]

                # Account for log plotting
                if self._getyaxistype() == 'log':
                    y.loc[y <= 0] = np.NaN

                if 'line_style' in plot_options['line_opt']:
                    if plot_options['line_opt']['line_style'] == 'scatter':
                        Glyph = X
                    else:
                        Glyph = Line
                else:
                    Glyph = Line
                # Until what size of a data gap are we removing nan values from the dataset? Set by the user
                # (default is to plot as bokeh would normally plot w/o worrying about data gap handling).
                limit = pytplot.tplot_opt_glob['data_gap']
                if limit != 0:
                    # Grabbing the times associated with nan values (nan_values), and the associated "position" of those
                    # keys in the dataset list (nan_keys)
                    nan_values = y[y.isnull().values].index.tolist()
                    nan_keys = [y.index.tolist().index(j) for j in nan_values]

                    nans = dict(zip(nan_keys, nan_values))

                    count = 0   # Keeping a count of how big of a time gap we have
                    consec_list = list()  # List of consecutive nan values (composed of indices for gaps not bigger than
                    # the user-specified data gap)
                    for val in range(len(nan_keys)):
                        # Avoiding some weird issues with going to the last data point in the nan dictionary keys
                        if val != (len(nan_keys)-1):
                            # Difference between one index and another - if consecutive indices, the diff will be 1
                            diff = abs(nan_keys[val] - nan_keys[val+1])
                            # calculate time accumulated from one index to the next
                            t_now = nan_values[val]
                            t_next = nan_values[val + 1]
                            time_accum = abs(t_now - t_next)
                            # If we haven't reached the allowed data gap, just keep track of how big of a gap we're at,
                            # and the indices in the gap
                            if diff == 1 and count < limit:
                                count += time_accum
                                consec_list.append(nan_keys[val])
                            # This triggers when we initially exceed the allowed data gap
                            elif diff == 1 and count >= limit:
                                pass
                            # When we find that the previous index and the current one are not consecutive, stop adding to
                            # the consec_list/overall_list (if applicable), and start over the count of time accumulated
                            # in a gap, as well as the consecutive list of time values with nans
                            elif diff != 1:
                                # Restart the count and add the current val to the list of nan values to remove
                                count = 0
                                consec_list.append(nan_keys[val])

                    times = x.tolist()
                    for elem in consec_list:
                        # Unless the data gap was big enough, we need to remove nan values from the data,
                        # otherwise bokeh will automatically NOT interpolate (the exact opposite of behavior in
                        # pyqtgraph, which ALWAYS interpolates...).
                        times.remove(nans[elem]*1000.0)
                        del y[nans[elem]]
                        del corrected_time[corrected_time.index(tplot_utilities.int_to_str(nans[elem]))]

                    # Data to be plotted
                    line_source = ColumnDataSource(data=dict(x=times, y=y, corrected_time=corrected_time))
                else:
                    # Data to be plotted
                    line_source = ColumnDataSource(data=dict(x=x, y=y, corrected_time=corrected_time))
                if self.auto_color:
                    line = Glyph(x='x', y='y', line_color=self.colors[self.linenum % len(self.colors)])
                else:
                    line = Glyph(x='x', y='y')
                if Glyph == Line:
                    if 'line_style' not in plot_options['line_opt']:
                        if line_style is not None:
                            line.line_dash = line_style[self.linenum % len(line_style)]
                    else:
                        line.line_dash = plot_options['line_style']
                self.lineglyphs.append(self.fig.add_glyph(line_source, line))
                self.linenum += 1

    def _addhoverlines(self):
        # Add tools
        hover = HoverTool()
        hover.tooltips = [("Time", "@corrected_time"), ("Value", "@y")]
        self.fig.add_tools(hover)

    def _addlegend(self):
        # Add the Legend if applicable
        if 'legend_names' in pytplot.data_quants[self.tvar_name].attrs['plot_options']['yaxis_opt']:
            legend_names = pytplot.data_quants[self.tvar_name].attrs['plot_options']['yaxis_opt']['legend_names']
            if len(legend_names) != self.linenum:
                print("Number of lines do not match length of legend names")
            legend = Legend()
            legend.location = (0, 0)
            legend_items = []
            j = 0
            for legend_name in legend_names:
                legend_items.append((legend_name, [self.lineglyphs[j]]))
                j = j+1
                if j >= len(self.lineglyphs):
                    break
            legend.items = legend_items
            legend.label_text_font_size = str(pytplot.data_quants[self.tvar_name].attrs['plot_options']['extras']['char_size'])+'pt'
            legend.border_line_color = None
            legend.glyph_height = int(self.fig.plot_height / (len(legend_items) + 1))
            self.fig.add_layout(legend, 'right')

max_stars_repo_name: xnchu/PyTplot

max_issues_repo_name: xnchu/PyTplot

max_forks_repo_name: xnchu/PyTplot


text:
#!/usr/bin/env python3

from tqdm import trange
import numpy as np
from collections import defaultdict
import os
import os.path
import pandas as pd
import toml
from numpy import array as arr
from glob import glob
from scipy import optimize
import cv2

from .common import make_process_fun, find_calibration_folder, \
    get_video_name, get_cam_name, natural_keys

from aniposelib.cameras import CameraGroup

def proj(u, v):
    """Project u onto v"""
    return u * np.dot(v,u) / np.dot(u,u)

def ortho(u, v):
    """Orthagonalize u with respect to v"""
    return u - proj(v, u)

def get_median(all_points_3d, ix):
    pts = all_points_3d[:, ix]
    pts = pts[~np.isnan(pts[:, 0])]
    return np.median(pts, axis=0)

def correct_coordinate_frame(config, all_points_3d, bodyparts):
    """Given a config and a set of points and bodypart names, this function will rotate the coordinate frame to match the one in config"""
    bp_index = dict(zip(bodyparts, range(len(bodyparts))))
    axes_mapping = dict(zip('xyz', range(3)))

    ref_point = config['triangulation']['reference_point']
    axes_spec = config['triangulation']['axes']
    a_dirx, a_l, a_r = axes_spec[0]
    b_dirx, b_l, b_r = axes_spec[1]

    a_dir = axes_mapping[a_dirx]
    b_dir = axes_mapping[b_dirx]

    ## find the missing direction
    done = np.zeros(3, dtype='bool')
    done[a_dir] = True
    done[b_dir] = True
    c_dir = np.where(~done)[0][0]

    a_lv = get_median(all_points_3d, bp_index[a_l])
    a_rv = get_median(all_points_3d, bp_index[a_r])
    b_lv = get_median(all_points_3d, bp_index[b_l])
    b_rv = get_median(all_points_3d, bp_index[b_r])

    a_diff = a_rv - a_lv
    b_diff = ortho(b_rv - b_lv, a_diff)

    M = np.zeros((3,3))
    M[a_dir] = a_diff
    M[b_dir] = b_diff
    # form a right handed coordinate system
    if (a_dir,b_dir) in [(0,1), (2,0), (1,2)]:
        M[c_dir] = np.cross(a_diff, b_diff)
    else:
        M[c_dir] = np.cross(b_diff, a_diff)

    M /= np.linalg.norm(M, axis=1)[:,None]

    center = get_median(all_points_3d, bp_index[ref_point])

    all_points_3d_adj = all_points_3d.dot(M.T)
    center_new = get_median(all_points_3d_adj, bp_index[ref_point])
    all_points_3d_adj = all_points_3d_adj - center_new

    return all_points_3d_adj, M, center_new

def load_pose2d_fnames(fname_dict, offsets_dict=None, cam_names=None):
    if cam_names is None:
        cam_names = sorted(fname_dict.keys())
    pose_names = [fname_dict[cname] for cname in cam_names]

    if offsets_dict is None:
        offsets_dict = dict([(cname, (0,0)) for cname in cam_names])

    datas = []
    for ix_cam, (cam_name, pose_name) in \
            enumerate(zip(cam_names, pose_names)):
        dlabs = pd.read_hdf(pose_name)
        if len(dlabs.columns.levels) > 2:
            scorer = dlabs.columns.levels[0][0]
            dlabs = dlabs.loc[:, scorer]

        bp_index = dlabs.columns.names.index('bodyparts')
        joint_names = list(dlabs.columns.get_level_values(bp_index).unique())
        dx = offsets_dict[cam_name][0]
        dy = offsets_dict[cam_name][1]

        for joint in joint_names:
            dlabs.loc[:, (joint, 'x')] += dx
            dlabs.loc[:, (joint, 'y')] += dy

        datas.append(dlabs)

    n_cams = len(cam_names)
    n_joints = len(joint_names)
    n_frames = min([d.shape[0] for d in datas])

    # frame, camera, bodypart, xy
    points = np.full((n_cams, n_frames, n_joints, 2), np.nan, 'float')
    scores = np.full((n_cams, n_frames, n_joints), np.nan, 'float')

    for cam_ix, dlabs in enumerate(datas):
        for joint_ix, joint_name in enumerate(joint_names):
            points[cam_ix, :, joint_ix] = np.array(dlabs.loc[:, (joint_name, ('x', 'y'))])[:n_frames]
            try:
                scores[cam_ix, :, joint_ix] = np.array(dlabs.loc[:, (joint_name, ('likelihood'))])[:n_frames].ravel()
            except KeyError:
                pass

    return {
        'cam_names': cam_names,
        'points': points,
        'scores': scores,
        'bodyparts': joint_names
    }

def load_offsets_dict(config, cam_names, video_folder=None):
    ## TODO: make the recorder.toml file configurable
    # record_fname = os.path.join(video_folder, 'recorder.toml')

    # if os.path.exists(record_fname):
    #     record_dict = toml.load(record_fname)
    # else:
    #     record_dict = None
    #     # if 'cameras' not in config:
    #     # ## TODO: more detailed error?
    #     #     print("-- no crop windows found")
    #     #     return

    offsets_dict = dict()
    for cname in cam_names:
        # if record_dict is None:
        if 'cameras' not in config or cname not in config['cameras']:
            # print("W: no crop window found for camera {}, assuming no crop".format(cname))
            offsets_dict[cname] = (0, 0)
        else:
            offsets_dict[cname] = tuple(config['cameras'][cname]['offset'])
        # else:
        #     offsets_dict[cname] = record_dict['cameras'][cname]['video']['ROIPosition']

    return offsets_dict

def load_constraints(config, bodyparts, key='constraints'):
    constraints_names = config['triangulation'].get(key, [])
    bp_index = dict(zip(bodyparts, range(len(bodyparts))))
    constraints = []
    for a, b in constraints_names:
        assert a in bp_index, 'Bodypart {} from constraints not found in list of bodyparts'.format(a)
        assert b in bp_index, 'Bodypart {} from constraints not found in list of bodyparts'.format(b)
        con = [bp_index[a], bp_index[b]]
        constraints.append(con)
    return constraints

def triangulate(config,
                calib_folder, video_folder, pose_folder,
                fname_dict, output_fname):

    cam_names = sorted(fname_dict.keys())

    calib_fname = os.path.join(calib_folder, 'calibration.toml')
    cgroup = CameraGroup.load(calib_fname)

    offsets_dict = load_offsets_dict(config, cam_names, video_folder)

    out = load_pose2d_fnames(fname_dict, offsets_dict, cam_names)
    all_points_raw = out['points']
    all_scores = out['scores']
    bodyparts = out['bodyparts']

    cgroup = cgroup.subset_cameras_names(cam_names)

    n_cams, n_frames, n_joints, _ = all_points_raw.shape

    bad = all_scores < config['triangulation']['score_threshold']
    all_points_raw[bad] = np.nan

    if config['triangulation']['optim']:
        constraints = load_constraints(config, bodyparts)
        constraints_weak = load_constraints(config, bodyparts, 'constraints_weak')

        points_2d = all_points_raw
        scores_2d = all_scores

        points_shaped = points_2d.reshape(n_cams, n_frames*n_joints, 2)
        if config['triangulation']['ransac']:
            points_3d_init, _, _, _ = cgroup.triangulate_ransac(points_shaped, progress=True)
        else:
            points_3d_init = cgroup.triangulate(points_shaped, progress=True)
        points_3d_init = points_3d_init.reshape((n_frames, n_joints, 3))

        c = np.isfinite(points_3d_init[:, :, 0])
        if np.sum(c) < 20:
            print("warning: not enough 3D points to run optimization")
            points_3d = points_3d_init
        else:
            points_3d = cgroup.optim_points(
                points_2d, points_3d_init,
                constraints=constraints,
                constraints_weak=constraints_weak,
                # scores=scores_2d,
                scale_smooth=config['triangulation']['scale_smooth'],
                scale_length=config['triangulation']['scale_length'],
                scale_length_weak=config['triangulation']['scale_length_weak'],
                n_deriv_smooth=config['triangulation']['n_deriv_smooth'],
                reproj_error_threshold=config['triangulation']['reproj_error_threshold'],
                verbose=True)

        points_2d_flat = points_2d.reshape(n_cams, -1, 2)
        points_3d_flat = points_3d.reshape(-1, 3)

        errors = cgroup.reprojection_error(
            points_3d_flat, points_2d_flat, mean=True)
        good_points = ~np.isnan(all_points_raw[:, :, :, 0])
        num_cams = np.sum(good_points, axis=0).astype('float')

        all_points_3d = points_3d
        all_errors = errors.reshape(n_frames, n_joints)

        all_scores[~good_points] = 2
        scores_3d = np.min(all_scores, axis=0)

        scores_3d[num_cams < 1] = np.nan
        all_errors[num_cams < 1] = np.nan

    else:
        points_2d = all_points_raw.reshape(n_cams, n_frames*n_joints, 2)
        if config['triangulation']['ransac']:
            points_3d, picked, p2ds, errors = cgroup.triangulate_ransac(
                points_2d, min_cams=3, progress=True)

            all_points_picked = p2ds.reshape(n_cams, n_frames, n_joints, 2)
            good_points = ~np.isnan(all_points_picked[:, :, :, 0])

            num_cams = np.sum(np.sum(picked, axis=0), axis=1)\
                         .reshape(n_frames, n_joints)\
                         .astype('float')
        else:
            points_3d = cgroup.triangulate(points_2d, progress=True)
            errors = cgroup.reprojection_error(points_3d, points_2d, mean=True)
            good_points = ~np.isnan(all_points_raw[:, :, :, 0])
            num_cams = np.sum(good_points, axis=0).astype('float')

        all_points_3d = points_3d.reshape(n_frames, n_joints, 3)
        all_errors = errors.reshape(n_frames, n_joints)

        all_scores[~good_points] = 2
        scores_3d = np.min(all_scores, axis=0)

        scores_3d[num_cams < 2] = np.nan
        all_errors[num_cams < 2] = np.nan
        num_cams[num_cams < 2] = np.nan

    if 'reference_point' in config['triangulation'] and 'axes' in config['triangulation']:
        all_points_3d_adj, M, center = correct_coordinate_frame(config, all_points_3d, bodyparts)
    else:
        all_points_3d_adj = all_points_3d
        M = np.identity(3)
        center = np.zeros(3)

    dout = pd.DataFrame()
    for bp_num, bp in enumerate(bodyparts):
        for ax_num, axis in enumerate(['x','y','z']):
            dout[bp + '_' + axis] = all_points_3d_adj[:, bp_num, ax_num]
        dout[bp + '_error'] = all_errors[:, bp_num]
        dout[bp + '_ncams'] = num_cams[:, bp_num]
        dout[bp + '_score'] = scores_3d[:, bp_num]

    for i in range(3):
        for j in range(3):
            dout['M_{}{}'.format(i, j)] = M[i, j]

    for i in range(3):
        dout['center_{}'.format(i)] = center[i]

    dout['fnum'] = np.arange(n_frames)

    dout.to_csv(output_fname, index=False)

def process_session(config, session_path):
    pipeline_videos_raw = config['pipeline']['videos_raw']
    pipeline_calibration_results = config['pipeline']['calibration_results']
    pipeline_pose = config['pipeline']['pose_2d']
    pipeline_pose_filter = config['pipeline']['pose_2d_filter']
    pipeline_3d = config['pipeline']['pose_3d']

    calibration_path = find_calibration_folder(config, session_path)
    if calibration_path is None:
        return

    if config['filter']['enabled']:
        pose_folder = os.path.join(session_path, pipeline_pose_filter)
    else:
        pose_folder = os.path.join(session_path, pipeline_pose)

    calib_folder = os.path.join(calibration_path, pipeline_calibration_results)
    video_folder = os.path.join(session_path, pipeline_videos_raw)
    output_folder = os.path.join(session_path, pipeline_3d)

    pose_files = glob(os.path.join(pose_folder, '*.h5'))

    cam_videos = defaultdict(list)

    for pf in pose_files:
        name = get_video_name(config, pf)
        cam_videos[name].append(pf)

    vid_names = cam_videos.keys()
    vid_names = sorted(vid_names, key=natural_keys)

    if len(vid_names) > 0:
        os.makedirs(output_folder, exist_ok=True)

    for name in vid_names:
        fnames = cam_videos[name]
        cam_names = [get_cam_name(config, f) for f in fnames]
        fname_dict = dict(zip(cam_names, fnames))

        output_fname = os.path.join(output_folder, name + '.csv')

        print(output_fname)

        if os.path.exists(output_fname):
            continue

        try:
            triangulate(config,
                        calib_folder, video_folder, pose_folder,
                        fname_dict, output_fname)
        except ValueError:
            import traceback, sys
            traceback.print_exc(file=sys.stdout)

triangulate_all = make_process_fun(process_session)

ext: py

max_stars_repo_name: SjoerdBruijn/anipose

max_issues_repo_name: SjoerdBruijn/anipose

max_forks_repo_name: SjoerdBruijn/anipose

reason: import numpy,from numpy,from scipy


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Jun  2 15:33:35 2020

@author: matthew
"""

import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt
import xarray as xr
import sys

sys.path.append("/home/matthew/university_work/01_blind_signal_separation_python/13_ICASAR/ICASAR_GitHub")                  # location of ICASAR functions
from ICASAR_functions import ICASAR

sys.path.append("/home/matthew/university_work/python_stuff/python_scripts")
#from small_plot_functions import matrix_show #, col_to_ma, r2_arrays_to_googleEarth
# from synth_ts import *
# from small_plot_functions import low_resolution_ifgs
# from small_plot_functions import *
# from ifg_plotting_funcs import *
from insar_tools import files_to_daisy_chain, get_daisy_chain, mask_nans, r3_to_r2, baseline_from_names
from S1VP_auxilliary_functions import r2_arrays_to_googleEarth
#from S1VP_auxilliary_functions import r2_arrays_to_googleEarth, get_daisy_chain, mask_nans, r3_to_r2, baseline_from_names

def xarray_to_numpy(xarray_file):
    """ Open one of Milan's xarray files and extract ifgs (or coherence)
    """

    import numpy as np
    import xarray as xr

    dUNW = xr.open_dataset(xarray_file)

    # 0 get the data
    n_ifgs = len(dUNW.data_vars)
    nx = dUNW.dims['lon']
    ny = dUNW.dims['lat']

    phUnw_r3 = np.zeros((n_ifgs, ny, nx))                                              # initaite rank 3
    phUnw_r2 = np.zeros((n_ifgs, (ny*nx)))                                              # initaite and simpler rank2, with ifgs as rows
    ifg_names = []

    for ifg_n, data_variable in enumerate(dUNW.data_vars):
        phUnw_r3[ifg_n,] = dUNW[data_variable]
        phUnw_r2[ifg_n,] = np.ravel(dUNW[data_variable])
        ifg_names.append(data_variable)

    # 1: get the lats and lons
    lons = np.asarray(dUNW.coords['lon'])
    lats = np.asarray(dUNW.coords['lat'])

    return phUnw_r3, phUnw_r2, ifg_names, lons, lats

def tcs_plot(tcs, baselines_cs):
    """Pot the time courses to be easier to see (each gets its own row)
    """

    n_tcs = tcs.shape[1]                                                             # each new tc is a colum

    f, axes = plt.subplots(n_tcs, 1)
    f.canvas.set_window_title('Time courses plot')

    for n_tc in range(n_tcs):
        axes[n_tc].scatter(baselines_cs, tcs[:,n_tc])
        axes[n_tc].plot(baselines_cs, tcs[:,n_tc])
        axes[n_tc].axhline(0, c = 'k')

def baselines_cs_with_gaps(ifg_names):
    """Quick and dirty function to get the baslines for a time sereis with gaps
    """
    import numpy as np
    from datetime import datetime

    baselines_cs = []
    master_ifg1 = datetime.strptime(ifg_names[0][:8], '%Y%m%d')
    for ifg_name in ifg_names:
        slave_ifgn = datetime.strptime(ifg_name[9:], '%Y%m%d')
        baselines_cs.append((slave_ifgn - master_ifg1).days)

    baselines_cs = np.asarray(baselines_cs)
    return baselines_cs

def col_to_ma(col, pixel_mask):
    """ A function to take a column vector and a 2d pixel mask and reshape the column into a masked array.
    Useful when converting between vectors used by BSS methods results that are to be plotted
    Inputs:
        col | rank 1 array |
        pixel_mask | array mask (rank 2)
    Outputs:
        source | rank 2 masked array | colun as a masked 2d array
    2017/10/04 | collected from various functions and placed here.
    """
    import numpy.ma as ma
    import numpy as np

    source = ma.array(np.zeros(pixel_mask.shape), mask = pixel_mask )
    source.unshare_mask()
    source[~source.mask] = col.ravel()
    return source

#%% Things to set

ref_region = (10,10)                                                # set a reference region, in x then y, from top left (ie like a matrix)
point_interest = (42, 30)                                           # x then y, from top left?

ICASAR_settings = {"n_comp" : 6,                                    # number of components to recover with ICA (ie the number of PCA sources to keep)
                    "bootstrapping_param" : (200, 0),
                    "hdbscan_param" : (35, 10),                        # (min_cluster_size, min_samples)
                    "tsne_param" : (30, 12),                       # (perplexity, early_exaggeration)
                    "ica_param" : (1e-2, 150),                     # (tolerance, max iterations)
                    "figures" : 'png+window',
                    "ge_kmz"    :  True}                            # make a google earth .kmz of the ICs

#%% Open the xarray files for unwrapped phase and coherence.
#data = xr.open_dataset('unw.nc')

phUnw_r3, _, ifg_names, lons, lats = xarray_to_numpy('unw.nc')
coh_r3, _, _, _, _ = xarray_to_numpy('coh.nc')

print('Revesring the order of the latitudes (bit of a fudge) ')
lats = lats[::-1]

mean_coh = np.mean(coh_r3, axis = 0)
r2_arrays_to_googleEarth(mean_coh[np.newaxis,], lons, lats, layer_name_prefix = 'layer', kmz_filename = 'mean_coherence')

#%% Set a reference region

_, ny, nx = phUnw_r3.shape                                                                  # get size of ifg, in pixels
ref_pixel_value = phUnw_r3[:,ref_region[1], ref_region[0]][:, np.newaxis, np.newaxis]
ref_pixel_r3 = np.repeat(np.repeat(ref_pixel_value, ny, 1), nx, 2)                          # duplicate to be the same size as all the ifgs.
phUnw_r3 -= ref_pixel_r3                                                                    # remove from all ifgs so ref pixel is always 0

#%% Do ICASAR with all ifgs

# phUnw_r3_ma, n_nan = mask_nans(phUnw_r3, figures = True, threshold = 2)                          # mask any nans, and remove ifgs with more than threshold % of pixels as nans
# phUnw_r2 = r3_to_r2(phUnw_r3_ma)                                                                        # convert to rank 2 (ifgs as row vectors), ready for ICASAR

# sources, tcs, residual, Iq, n_clusters, S_all_info, phUnw_mean = ICASAR(phUnw_r2['ifgs'], mask = phUnw_r2['mask'],
#                                                                         lons = lons, lats = lats, **ICASAR_settings,
#                                                                         out_folder = './ICASAR_outputs_all_ifgs/')

#%% Do ICASAR with only the short baseline ifgs

phUnw_r3, ifg_names = get_daisy_chain(phUnw_r3, ifg_names)                                        # get the daisy chain ifgs.

# Milan - 2020/06/30 addition  PhUnwr3 is 83x110x120, as expected, but there are some pixels that are nan occasioanly in the off interferogram, and have to be removed from the entire time series.

phUnw_r3_ma, n_nan = mask_nans(phUnw_r3 , threshold = 2)                                          # mask any nans (and drop any ifgs that more than 2% of them are nans)

# now PhUnw_r3_ma is a masked array version of phUnw_r3 - same size (83x110x120), but some pixels in that are masked

phUnw_r2 = r3_to_r2(phUnw_r3_ma )                                                                  # convert to rank 2 (ie row vectors, which are needed for ICA)

# if we want to see the pixels I had to drop (ie the mask):
f, ax =plt.subplots()
ax.imshow(phUnw_r2['mask'])
ax.set_title('nans mask')

# and having a look at the size of things
n_pixels_masked = np.sum(phUnw_r2['mask'])                  # should be 60
# 60 pixels are masked.  (110x120) = 13200
#       |->                 -60    = 13140
# phUnw_r2['ifgs'] shape is     83 x 13140   (so the samenumber of pixels)

# to get back to a rank 2 masked array, use a little fuctnion called col_to_ma.
test_ifg = col_to_ma(phUnw_r2['ifgs'][0,], phUnw_r2['mask'])
# and continuing as before.

baselines_cs  = baselines_cs_with_gaps(ifg_names)                                                  # get the cumulative baselines (ie 6 12 18 24 etc if all 6 day)

sources, _,  tcs, residual, Iq, n_clusters, S_all_info, phUnw_mean  = ICASAR(phUnw_r2 ['ifgs'], mask = phUnw_r2 ['mask'],
                                                                                lons = lons, lats = lats, **ICASAR_settings,
                                                                                out_folder = './ICASAR_outputs/')

tcs_plot(np.cumsum(tcs , axis = 0), baselines_cs )                                                  # plot the cumulative time coursses

#%% Point of interest

poi_ts = phUnw_r3_ma [:, point_interest[1], point_interest[0]]
poi_ts_cs = np.cumsum(poi_ts)

f, ax = plt.subplots(1,1)
ax.scatter(np.arange(poi_ts_cs.shape[0]), poi_ts_cs)
ax.axhline(0, c = 'k')
plt.suptitle(f'Time series for one point (x:{point_interest[0]} y:{point_interest[1]}) ')
f.canvas.set_window_title('Time series for one point')

#%% Save the ifgs as pngs.

# from small_plot_functions import matrix_show

# for ifg_n, ifg_ma in enumerate(phUnw_r3_ma):
#     matrix_show(ifg_ma, ifg_names[ifg_n], save_path = "./ifgs/")

reason: import numpy


text:
"""
use cross validation to plot mean ROC curve, show std

ref:

https://scikit-learn.org/stable/auto_examples/model_selection/plot_roc_crossval.html#sphx-glr-auto-examples-model-selection-plot-roc-crossval-py

Note that you have to tune the parameters yourself

"""
from Bio.Seq import Seq
from Bio import SeqIO
from scipy import interp
import argparse
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
# import xgboost as xgb
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.colors as clr
import numpy as np
from matplotlib.colors import ListedColormap
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.colors import ListedColormap, LinearSegmentedColormap
import pandas as pd
import matplotlib.pylab as plt
import numpy as np
import scipy
import seaborn as sns
import glob
from sklearn.model_selection import KFold,StratifiedKFold
from sklearn import model_selection
from sklearn.linear_model import LogisticRegression,RidgeClassifier,SGDClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.ensemble import RandomForestClassifier
# from mlxtend.classifier import StackingCVClassifier
import os
import warnings
from sklearn.metrics import roc_curve,roc_auc_score,average_precision_score,precision_recall_curve
from sklearn.datasets import load_iris
# from mlxtend.classifier import StackingCVClassifier
# from mlxtend.feature_selection import ColumnSelector
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LogisticRegression
import warnings
warnings.filterwarnings('ignore')
from sklearn.exceptions import ConvergenceWarning
warnings.simplefilter(action='ignore', category=FutureWarning)
warnings.simplefilter(action='ignore', category=ConvergenceWarning)
from sklearn.ensemble import RandomForestRegressor,GradientBoostingRegressor,RandomForestClassifier
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import GridSearchCV
from sklearn.metrics.scorer import make_scorer
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
import uuid
from sklearn.base import TransformerMixin
from sklearn.datasets import make_regression
from sklearn.pipeline import Pipeline, FeatureUnion
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestRegressor,GradientBoostingClassifier
from sklearn.neighbors import KNeighborsRegressor
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.linear_model import LinearRegression, Ridge
import scipy
import numpy as np
from sklearn.model_selection import RandomizedSearchCV
from sklearn.model_selection import LeaveOneOut
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_absolute_error
from sklearn import linear_model
from sklearn.kernel_ridge import KernelRidge
from sklearn.svm import SVR,LinearSVC
from sklearn.neighbors import KNeighborsRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import Lasso
from sklearn.linear_model import Ridge,Lars,BayesianRidge
from copy import deepcopy as dp
from sklearn.datasets import make_moons, make_circles, make_classification
from sklearn.neural_network import MLPClassifier
from sklearn.neighbors import KNeighborsClassifier,RadiusNeighborsClassifier
from sklearn.svm import SVC
from sklearn.gaussian_process import GaussianProcessClassifier
from sklearn.gaussian_process.kernels import RBF
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier, AdaBoostClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis
from sklearn.gaussian_process import GaussianProcessClassifier
# from xgboost import XGBClassifier

def read_fasta(f):
	my_dict = {}
	for r in SeqIO.parse(f, "fasta"):
		my_dict[r.id] = str(r.seq).upper()
	return my_dict

def sklearn_RF(par=False):

	est = RandomForestClassifier(n_estimators=1000,random_state=0,warm_start=False,n_jobs=-1,class_weight={1:10,0:1})
	if par:
		est = RandomForestClassifier(**par)
	myDict = {}
	return est, myDict

def plot_top_features(reg,X,y,output):

	current_feature_df = pd.DataFrame()
	current_feature_df['features'] = X.columns.tolist()

	reg.fit(X,y)
	try:
		current_feature_df['score'] = list(reg.feature_importances_)
	except:
		try:
			current_feature_df['score'] = list(reg.coef_)
		except:
			current_feature_df['score'] = list(reg.coef_[0])

	current_feature_df = current_feature_df.sort_values('score',ascending=False)

	plt.figure(figsize=(len(current_feature_df['features']*2),8))
	sns.barplot(x=current_feature_df['features'],y=current_feature_df['score'] )
	plt.xticks(rotation=90)
	plt.xlabel("")
	plt.ylabel("Feature importance")
	plt.savefig("%s_feature_importance.pdf"%(output), bbox_inches='tight')

def to_true_label(pos,neg):
	pos_fasta = read_fasta(pos)
	neg_fasta = read_fasta(neg)
	df = pd.DataFrame()
	df['true'] = [1]*len(pos_fasta.keys())+[-1]*len(neg_fasta.keys())
	df.index = pos_fasta.keys()+neg_fasta.keys()
	return df

def gkm_SVM_fit_transform(pos_train,neg_train,pos_test,neg_test,d=3):
	"""From HemTools"""

	addon_string = str(uuid.uuid4()).split("-")[-1]
	train_out = addon_string
	classify_out = addon_string + ".classify.out"
	train_command = "gkmtrain -w 10 %s %s %s >/dev/null 2>&1"%(pos_train,neg_train, train_out)

	classify_input = addon_string + ".classify_input.fa"
	combine_pos_neg_command = "cat %s %s > %s"%(pos_test,neg_test,classify_input)
	classify_command = "gkmpredict %s %s.model.txt %s >/dev/null 2>&1"%(classify_input,train_out,classify_out)
	## run
	os.system(train_command)
	os.system(combine_pos_neg_command)
	os.system(classify_command)

	y_pred = pd.read_csv(classify_out,sep="\t",header=None,index_col=0)
	y_pred.columns = ['pred']
	os.system("rm %s*"%(addon_string))
	return y_pred

def write_fasta(file_name,myDict):
	out = open(file_name,"wt")
	for k in myDict:
		out.write(">"+k+"\n")
		out.write(myDict[k]+"\n")
	out.close()
def get_fasta_file(X_train,y_train,X_test,y_test):
	addon_string = str(uuid.uuid4()).split("-")[-1]
	pos_train_file = addon_string+"_pos_train_file.fa"
	neg_train_file = addon_string+"_neg_train_file.fa"
	pos_test_file = addon_string+"_pos_test_file.fa"
	neg_test_file = addon_string+"_neg_test_file.fa"

	write_fasta(pos_train_file,X_train.loc[y_train[y_train==1].index.tolist()].to_dict())
	write_fasta(pos_test_file,X_test.loc[y_test[y_test==1].index.tolist()].to_dict())
	write_fasta(neg_train_file,X_train.loc[y_train[y_train==0].index.tolist()].to_dict())
	write_fasta(neg_test_file,X_test.loc[y_test[y_test==0].index.tolist()].to_dict())

	# Y[Y==1]
	return pos_train_file,neg_train_file,pos_test_file,neg_test_file,addon_string

def find_neighbor(x,y):
	"""if y is within +-10bp of x"""
	# print (x,y)
	chr_x = x[:-1].split(":")[0]
	start_x = int(x[:-1].split(":")[-1].split("-")[0])
	chr_y = y[:-1].split(":")[0]
	start_y = int(y[:-1].split(":")[-1].split("-")[0])
	# print (chr_x,chr_y,start_x,start_y)
	if chr_x != chr_y:
		return False
	if abs(start_y-start_x)<=10:
		return True
	return False

def random_sample2(y):
	out = []
	cv=1
	used_seed=[]
	for i in range(cv):
		seed = y[~y.index.isin(used_seed)].sample(frac=0.1).index.tolist()
		print ("seed length",len(seed))
		unseed = list(set(y.index.tolist())-set(seed))
		for u in unseed:
			for s in seed:
				if find_neighbor(s,u):
					seed.append(u)
					break
		train = list(set(y.index.tolist())-set(seed))
		used_seed+=seed
		out.append([train,seed])
		print ("train size: %s "%(len(train)))
	return out

def random_sample(y,random_state):
	out = []
	cv=1
	for i in range(cv):
		seed = y.sample(frac=0.1,random_state=random_state).index.tolist()
		while True:
			flag_list = []
			unseed = list(set(y.index.tolist())-set(seed))
			for u in unseed:
				flag = 0
				for s in seed:
					if find_neighbor(s,u):
						seed.append(u)
						flag = 1
						break
				flag_list.append(flag)
			if sum(flag_list) == 0:
				break

		train = list(set(y.index.tolist())-set(seed))
		out.append([train,seed])
		# print ("train size: %s "%(len(train)))
	return out

def simple_CV_evaluation(model,X,y,k,random_state):

	# outer = StratifiedKFold(n_splits=5,shuffle=False)
	## total random sampling gave over-estimated scores because neighbor As (with almost the same features) are random splited.
	my_pred=[]
	my_true=[]

	auPRC_list = []
	auROC_list = []
	for train_index, test_index in random_sample(y,random_state):
		print ("train size: %s test size: %s"%(len(train_index),len(test_index)))
		X_train, X_test = X.loc[train_index], X.loc[test_index]
		y_train, y_test = y.loc[train_index], y.loc[test_index]
		if k == "RF-motif-ATAC":
			current_model = dp(model)
			current_model.fit(X_train,y_train)
			pred_y = current_model.predict_proba(X_test)
			pred_y = [x[1] for x in pred_y]
		elif k == "LS-GKM":

			pos_train_file,neg_train_file,pos_test_file,neg_test_file,addon_string = get_fasta_file(X_train,y_train,X_test,y_test)

			pred_y = gkm_SVM_fit_transform(pos_train_file,neg_train_file,pos_test_file,neg_test_file)
			pred_y = pred_y.loc[y_test.index.tolist()]['pred'].tolist()
			os.system("rm %s*"%(addon_string))

		else:

			pred_y = X_test.tolist()

		y_test = y_test.tolist()
		my_pred += pred_y
		my_true += y_test
		try:
			auROC = roc_auc_score(y_test,pred_y)
			auPRC = average_precision_score(y_test,pred_y)
			print ("model %s auPRC: %s. auROC: %s"%(k,auPRC,auROC))
			auPRC_list.append(auPRC)
			auROC_list.append(auROC)
		except:
			pass

	df = pd.DataFrame()
	df['true']=my_true
	df['pred']=my_pred
	df['label'] = k
	return df,auROC_list,auPRC_list

def plot_auROC_multi(df,color_dict):
	sns.set_style("white")
	plt.figure()

	for s,d in df.groupby('label'):
		plot_df = pd.DataFrame()
		x_predict,y_predict,_ = roc_curve(d['true'],d['pred'])
		auc = roc_auc_score(d['true'],d['pred'])
		# print (auc)
		plot_df['x'] = x_predict
		plot_df['y'] = y_predict
		sns.lineplot(data=plot_df,x="x",y="y",ci=0,label="%s AUC:%.2f"%(s,auc),color=color_dict[s])
	plt.plot([0, 1], [0, 1], 'k--')
	plt.xlim(0,1)
	plt.ylim(0,1)
	plt.xlabel('False positive rate')
	plt.ylabel('True positive rate')
	plt.title('ROC curve')
	plt.legend(loc='best',title="")
	# plt.savefig("auROC.png")
	plt.savefig("auROC.pdf", bbox_inches='tight')
	plt.close()

def plot_auPRC_multi(df,color_dict):
	sns.set_style("white")
	plt.figure()

	for s,d in df.groupby('label'):
		plot_df = pd.DataFrame()
		y_predict,x_predict,_ = precision_recall_curve(d['true'],d['pred'])
		auc = average_precision_score(d['true'],d['pred'])
		# print (auc)
		plot_df['x'] = x_predict
		plot_df['y'] = y_predict
		# sns.lineplot(data=plot_df,x="x",y="y",ci=0,label="%s AUC:%.2f"%(s,auc),color=color_dict[s])
		plt.step(x_predict, y_predict, where='post',label="%s AUC:%.2f"%(s,auc),color=color_dict[s])
	# plt.plot([0, 1], [0, 1], 'k--')
	plt.xlim(0,1)
	plt.ylim(0,1)
	plt.xlabel('Recall')
	plt.ylabel('Precision')
	plt.title('Precision-Recall curve')
	plt.legend(loc='best',title="")
	# plt.savefig("auPRC.png")
	plt.savefig("auPRC.pdf", bbox_inches='tight')
	plt.close()

def define_high_low(x):
	if x == 0:
		return 0
	elif x>=50:
		return 1
	else:
		print ("value out of range. Error!")
		return -1

def boxplot_paired_t_test(df,color_dict,ylabel,output):
	sns.set_style("whitegrid")
	myMin = df.min().min()
	myMax = df.max().max()
	plot_df = pd.melt(df)
	plt.figure()
	ax=sns.boxplot(x="variable",y='value',data=plot_df,palette =color_dict,linewidth=3,order=['RF-motif-ATAC','LS-GKM','CADD','DeepSEA'])
	for patch in ax.artists:
		r, g, bb, _ = patch.get_facecolor()
		patch.set_facecolor((r, g, bb, .3))

	sns.swarmplot(x="variable", y="value", data=plot_df,palette =color_dict,order=['RF-motif-ATAC','LS-GKM','CADD','DeepSEA'])
	unit=0.01

	y=myMax+unit*1
	h=unit*1.5
	plt.plot([0, 0, 1, 1], [y, y+h, y+h, y], lw=1.5, c="black")
	plt.text(0.5, y+h+unit, "Paired T-test: %.2E" % scipy.stats.ttest_rel(df['RF-motif-ATAC'],df['LS-GKM']).pvalue, ha='center', va='bottom', color="black")

	# y=myMax+unit*6
	# h=unit*1.5
	# plt.plot([0, 0, 2, 2], [y, y+h, y+h, y], lw=1.5, c="black")
	# plt.text(0.5, y+h+unit, "Paired T-test: %.2E" % scipy.stats.ttest_rel(df['RF-motif-ATAC'],df['CADD']).pvalue, ha='center', va='bottom', color="black")
	print ("RF vs. CADD Paired T-test: %.2E" % scipy.stats.ttest_rel(df['RF-motif-ATAC'],df['CADD']).pvalue)
	print ("RF vs. DeepSEA Paired T-test: %.2E" % scipy.stats.ttest_rel(df['RF-motif-ATAC'],df['DeepSEA']).pvalue)
	# y=myMax+unit*11
	# h=unit*1.5
	# plt.plot([0, 0, 3, 3], [y, y+h, y+h, y], lw=1.5, c="black")
	# plt.text(0.5, y+h+unit, "Paired T-test: %.2E" % scipy.stats.ttest_rel(df['RF-motif-ATAC'],df['DeepSEA']).pvalue, ha='center', va='bottom', color="black")

	plt.ylim(myMin-unit*5,myMax+0.1)
	plt.ylabel(ylabel)
	plt.savefig("%s.pdf"%(output), bbox_inches='tight')

## define colors for line plot and boxplot
color_dict = {}
color_dict['RF-motif-ATAC']="#fa2111"
color_dict['DeepSEA']="#213fff"
color_dict['CADD']="#00bd3c"
color_dict['LS-GKM']="#813685"

## define X Y
df = pd.read_csv("ML_data.csv",index_col=0)
target = "HbFBase"
df = df.drop(['GERP'],axis=1)
df['label'] = [define_high_low(x) for x in df[target]]
Y = df['label']
print ("Negative Size:",Y[Y==0].shape)
print ("Positive Size:",Y[Y==1].shape)
X = df.drop([target,'label'],axis=1)
RF_features = []
for c in X.columns:
	if "motif" in c:
		RF_features.append(c)
X_RF = X[RF_features]
X_DeepSEA = X['DeepSEA']
X_CADD = X['CADD']
X_gkm = X['seq']

## RF model
model,params=sklearn_RF()

auROC_dict={}
auPRC_dict={}
for k in color_dict:
	auROC_dict[k]=[]
	auPRC_dict[k]=[]
df_list = []
X_dict = {}
X_dict['RF-motif-ATAC']=X_RF
X_dict['DeepSEA']=X_DeepSEA
X_dict['CADD']=X_CADD
X_dict['LS-GKM']=X_gkm
## RF top features
plot_top_features(dp(model),X_RF,Y,"motif-ATAC")

from joblib import Parallel, delayed
## CV evaluation
def parallel_unit(model,X_dict,Y):
	random_state = np.random.randint(0,999999)
	out_dict = {}
	for k in X_dict:
		ddf,roc_list,prc_list = simple_CV_evaluation(model,X_dict[k],Y,k,random_state)
		out_dict[k]=[ddf,roc_list,prc_list]
	return out_dict
out_list = Parallel(n_jobs=-1)(delayed(parallel_unit)(model,X_dict,Y) for i in range(300))

# for i in range(500):
	# random_state = np.random.randint(0,999999)
	# for k in X_dict:
		# ddf,roc_list,prc_list = simple_CV_evaluation(model,X_dict[k],Y,k,random_state)
		# auROC_dict[k]+=roc_list
		# auPRC_dict[k]+=prc_list
		# df_list.append(ddf)
for item in out_list:
	for k in X_dict:
		auROC_dict[k]+=item[k][1]
		auPRC_dict[k]+=item[k][2]
		df_list.append(item[k][0])
plot_df = pd.concat(df_list)
plot_auROC_multi(plot_df,color_dict)
plot_auPRC_multi(plot_df,color_dict)

boxplot_paired_t_test(pd.DataFrame(auROC_dict),color_dict,"Area under ROC curve","auROC_boxplot")
boxplot_paired_t_test(pd.DataFrame(auPRC_dict),color_dict,"Area under Precision-Recall curve","auPRC_boxplot")
# auROC
# RF vs. CADD Paired T-test: 4.92E-84
# RF vs. DeepSEA Paired T-test: 2.40E-92
# auPRC
# RF vs. CADD Paired T-test: 4.52E-75
# RF vs. DeepSEA Paired T-test: 1.19E-58

reason: import numpy,import scipy,from scipy


text:
# Author: Mikita Sazanovich

from collections import namedtuple
from typing import List
import os
import pathlib
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

Experiment = namedtuple('Experiment', ['name', 'color', 'run_csvs'])

TOTAL_STEPS = 300000
BIN_SIZE = 30000
assert TOTAL_STEPS % BIN_SIZE == 0
TOTAL_BINS = TOTAL_STEPS // BIN_SIZE + 1

def plot_with_variance(values, color, label):
    xs = np.arange(TOTAL_BINS) * BIN_SIZE
    ymeans = np.zeros(TOTAL_BINS)
    ystderrs = np.zeros(TOTAL_BINS)
    for i in range(TOTAL_BINS):
        xvalues = values[i]
        if len(xvalues) == 0:
            continue
        ymean = np.mean(xvalues)
        ystd = np.std(xvalues)
        ystderr = ystd / np.sqrt(len(xvalues))
        ymeans[i], ystderrs[i] = ymean, ystderr
    plt.plot(xs, ymeans, color, label=label)
    plt.fill_between(xs, ymeans - ystderrs, ymeans + ystderrs, color=color, alpha=.2)

def init_the_plot():
    plt.figure(figsize=(10, 8))

def draw_the_plot(y_limit):
    fontsize = 16
    plt.xlabel('Training Steps', fontsize=fontsize)
    plt.xticks(fontsize=fontsize)
    plt.ylabel('Episode Discounted Reward', fontsize=fontsize)
    plt.yticks(fontsize=fontsize)
    plt.axis([None, None, 0, y_limit])
    plt.legend(loc='upper left', prop={'size': fontsize})
    plt.show()

def draw_mean_var_graphs(csv_dir: str, experiments: List[Experiment]) -> None:
    init_the_plot()
    for experiment in experiments:
        values = []
        for i in range(TOTAL_BINS):
            values.append([])
        for filename in experiment.run_csvs:
            csv_file = os.path.join(csv_dir, filename)
            with open(csv_file) as csv_file:
                df = pd.read_csv(csv_file)
            for _, row in df.iterrows():
                step = row['Step']
                if step >= TOTAL_STEPS:
                    continue
                value = row['Value']
                bin_id = int(step / BIN_SIZE) + 1
                if bin_id >= TOTAL_BINS:
                    continue
                values[bin_id].append(value)
        plot_with_variance(values, color=experiment.color, label=experiment.name)
    draw_the_plot(y_limit=60)

def draw_smoothed_graphs(rewards_dir, filenames, labels):
    init_the_plot()
    smoothing_weight = 0.95
    rewards_dir = pathlib.Path(rewards_dir)
    for filename, label in zip(filenames, labels):
        csv_path = os.path.join(rewards_dir, filename)
        steps = []
        values = []
        last = None
        with open(str(csv_path), 'r') as csv_file:
            df = pd.read_csv(csv_file)
            for _, row in df.iterrows():
                step = row['Step']
                if step > TOTAL_STEPS:
                    continue
                value = row['Value']
                if last is None:
                    last = value
                smoothed_value = last * smoothing_weight + (1 - smoothing_weight) * value
                last = smoothed_value
                steps.append(step)
                values.append(smoothed_value)
        plt.plot(steps, values, label=label)
    draw_the_plot(y_limit=80)

def verify_results_existence(results_dir, id_filenames_dict):
    for filenames in id_filenames_dict.values():
        for filename in filenames:
            result_filename = os.path.join(results_dir, filename)
            assert os.path.exists(result_filename)

def main():
    results_dir = '/Users/niksaz/4-RnD/csv-results'
    id_filenames_dict = {
        '0-demo': [
            'run_20190504-no-demo-seed-50_summaries-tag-rewards.csv',
            'run_20190506-no-demo-seed-98_summaries-tag-rewards.csv',
            'run_20190506-no-demo-seed-79_summaries-tag-rewards.csv',
            'run_20190508-no-demo-seed-92_summaries-tag-rewards.csv',
        ],
        '1-demo': [
            'run_20190415-1-random-full-sigma-02_summaries-tag-rewards.csv',
            'run_20190416-1-random-full-sigma-02_summaries-tag-rewards.csv',
            'run_20190504-1-demo-1725_summaries-tag-rewards.csv',
            'run_20190505-1-demo-1738_summaries-tag-rewards.csv',
        ],
        '50-demos': [
            'run_20190419-50-demos-merged-sigma-02_summaries-tag-rewards.csv',
            'run_20190421-50-demos-merged-sigma-02_summaries-tag-rewards.csv',
            'run_20190429-50-demos-merged-seed-90_summaries-tag-rewards.csv',
            'run_20190503-50-demos-seed-32_summaries-tag-rewards.csv',
        ],
        '150-demos': [
            'run_20190423-150-demos-merged-sigma-02_summaries-tag-rewards.csv',
            'run_20190425-150-demos-merged-sigma-02_summaries-tag-rewards.csv',
            'run_20190425-150-demos-merged-sigma-02_summaries-tag-rewards.csv'
        ],
        '254-demos': [
            'run_20190413-all-merged-095-sigma-02_summaries-tag-rewards.csv',
            'run_20190417-all-merged-095-sigma-02_summaries-tag-rewards.csv',
            'run_20190418-all-merged-095-sigma-02_summaries-tag-rewards.csv',
            'run_20190507-all-merged-095-sigma-seed-32_summaries-tag-rewards.csv',
        ],
        'clustering': [
            'run_20190427-cluster-all-1700_summaries-tag-rewards.csv',
            'run_20190428-cluster-all-1700_summaries-tag-rewards.csv',
        ],
        'stochastic': [
            'run_20190419-50-demos-merged-sigma-02_summaries-tag-rewards.csv',
            'run_20190420-50-demos-merged-sigma-02_summaries-tag-rewards.csv',
            'run_20190421-50-demos-merged-sigma-02_summaries-tag-rewards.csv'
        ],
        'random': [
            'run-20191009-random-seed-7_summaries-tag-rewards.csv',
            'run-20191010-random-seed-46_summaries-tag-rewards.csv',
            'run-20191011-random-seed-10_summaries-tag-rewards.csv',
            'run-20191015-random-seed-78_summaries-tag-rewards.csv',
        ],
    }
    verify_results_existence(results_dir, id_filenames_dict)

    # Compiling the stochasticity visual proof graph
    draw_smoothed_graphs(results_dir, id_filenames_dict['stochastic'], [
        'Run 1',
        'Run 2',
        'Run 3',
    ])
    # Compiling the demo/no-demo graph
    draw_mean_var_graphs(results_dir, [
        Experiment('No demonstrations', 'black', id_filenames_dict['0-demo']),
        Experiment('Single demonstration', 'magenta', id_filenames_dict['1-demo']),
    ])
    # Compiling the comparison of random sampling with our approach
    draw_mean_var_graphs(results_dir, [
        Experiment('Our approach', '#F5B041', id_filenames_dict['50-demos']),
        Experiment('Random sampling', '#708090', id_filenames_dict['random'])
    ])
    # Compiling the different volumes graph
    draw_mean_var_graphs(results_dir, [
        Experiment('With 1 demonstration', 'magenta', id_filenames_dict['1-demo']),
        Experiment('With 50 demonstrations', '#F5B041', id_filenames_dict['50-demos']),
        # Experiment('With 150 demonstrations', 'cyan', id_filenames_dict['150-demos']),
        Experiment('With 254 demonstrations', 'blue', id_filenames_dict['254-demos']),
        # Experiment('With 254 demonstrations clustered', 'red', id_filenames_dict['clustering']),
    ])

if __name__ == '__main__':
    main()


text:
import numpy as np
import torch
import torch.nn.functional as F

M_DAYS = {
    0: 0,
    1: 31,
    2: 59,
    3: 90,
    4: 120,
    5: 151,
    6: 181,
    7: 212,
    8: 243,
    9: 273,
    10: 304,
    11: 334,
}

def calc_interval(date_lst):
    month = [int(x.split("-")[0]) for x in date_lst]
    day = [int(x.split("-")[1]) for x in date_lst]
    day = [M_DAYS[m - 1] + d for m, d in zip(month, day)]
    day = np.array([d - day[0] for d in day])
    return np.abs(day)

def load_all_factors(name_lst, keys, center, ref="filterlesions", minus=False):
    out = []
    mask = []
    for k in keys:
        if k == "multilesions" or k == "filterlesions" or k == "deltalesions":
            data = np.stack([np.load(x.replace(ref, k)) for x in name_lst], axis=-1).astype(np.float32)
        else:
            data = np.stack([np.load(x.replace(ref, k)[:-7] + ".npy") for x in name_lst], axis=-1).astype(np.float32)
        if k == "lungsegs":
            mask.append(data)
            continue
        elif k == "lungdists" or k == "vesseldists" or k == "lesiondists":
            if minus:
                data = 1.0 - data / 255
            else:
                data /= 255
        # data = (data - data.min()) / (data.max() - data.min()).clip(0.001, None)
        if "lesions" not in k and k != "vessels":
            data = (data - data.mean()) / data.std().clip(0.001, None) + np.random.randn(*(data.shape)) * 0.1
        out.append(data)
    out = [mask[0] * x for x in out]
    out = np.concatenate(out, axis=-1)
    return out, mask[0][..., center].astype(np.uint8)
    # return out + np.random.randn(*(out.shape)) * 0.1, mask[0][..., center].astype(np.uint8)

def reshape_flatten(input_arr, kernel_size=1):
    ndim = input_arr.ndim - 1
    if len(kernel_size) == 1:
        kernel_size = (kernel_size,) * ndim
    elif len(kernel_size) != ndim:
        raise ValueError(f"the kernel size is not consist with input array")

    vol_shape = input_arr.shape[:-1]
    if np.all(kernel_size == 1):
        return np.squeeze(input_arr.reshape(np.prod(vol_shape), -1))
    else:
        padding = [x // 2 for x in kernel_size]
        to_tensor = torch.from_numpy(input_arr.transpose((-1,) + tuple(range(ndim))))[None].contiguous()
        tensor_view = F.unfold(to_tensor, kernel_size, padding=padding).transpose(1, 2)
        tensor_view = tensor_view.view(np.prod(vol_shape), -1)
        return tensor_view.numpy()

def random_select_samples(x_pool, y_pool, mask=None, sampler_p=0.5):
    numel = len(y_pool) if mask is None else mask.sum()
    pos_N = (y_pool > 0).sum()
    neg_N = numel - pos_N
    #
    # pos_indices = np.where(y_pool > 0)
    # if mask is None or (y_pool == 0 & mask).sum() < pos_N:
    #     neg_indices = np.where(y_pool == 0)
    # else:
    #     neg_indices = np.where(y_pool == 0 & mask)
    # indices = np.arange(len(neg_indices[0]))
    # indices = np.random.choice(indices, min(len(indices), 4 * pos_N))
    # neg_indices = tuple([x[indices] for x in neg_indices])
    # x_samples = np.concatenate((x_pool[pos_indices], x_pool[neg_indices]), axis=0)
    # y_samples = np.concatenate((y_pool[pos_indices], y_pool[neg_indices]), axis=0)
    # return x_samples, y_samples, np.ones_like(y_samples)

    indices = np.where(mask > 0)
    x_samples = x_pool[indices]
    y_samples = y_pool[indices]
    mask = mask[indices]
    neg_weight = sampler_p * pos_N / neg_N
    sample_weight = (mask != 0) * sampler_p + (mask == 0) * neg_weight
    return x_samples, y_samples, sample_weight

    # neg_weight = sampler_p * pos_N / neg_N
    # sample_weight = (y_pool != 0) * sampler_p + (y_pool == 0) * neg_weight
    # if mask is not None:
    #     sample_weight *= mask > 0.5
    # return x_pool.reshape(y_pool.size, -1), y_pool.reshape(-1), sample_weight.reshape(-1)


text:
__author__ = 'zhengwang'

import numpy as np
import cv2
import socket

class VideoStreamingTest(object):
    def __init__(self, host, port):

        self.server_socket = socket.socket()
        self.server_socket.bind((host, port))
        self.server_socket.listen(0)
        self.connection, self.client_address = self.server_socket.accept()
        self.connection = self.connection.makefile('rb')
        self.host_name = socket.gethostname()
        self.host_ip = socket.gethostbyname(self.host_name)
        self.streaming()

    def streaming(self):

        try:
            print("Host: ", self.host_name + ' ' + self.host_ip)
            print("Connection from: ", self.client_address)
            print("Streaming...")
            print("Press 'q' to exit")

            # need bytes here
            stream_bytes = b' '
            while True:
                stream_bytes += self.connection.read(1024)
                first = stream_bytes.find(b'\xff\xd8')
                last = stream_bytes.find(b'\xff\xd9')
                if first != -1 and last != -1:
                    jpg = stream_bytes[first:last + 2]
                    stream_bytes = stream_bytes[last + 2:]
                    image = cv2.imdecode(np.frombuffer(jpg, dtype=np.uint8), cv2.IMREAD_COLOR)
                    cv2.imshow('image', image)

                    if cv2.waitKey(1) & 0xFF == ord('q'):
                        break
        finally:
            self.connection.close()
            self.server_socket.close()

if __name__ == '__main__':
    # host, port
    h, p = "192.168.1.100", 8000
    VideoStreamingTest(h, p)


text:
#!/usr/bin/env python

# ENCODE DCC compare signal to roadmap wrapper
# Author: Daniel Kim, Jin Lee (leepc12@gmail.com)

import warnings
from matplotlib import pyplot as plt
import sys
import os
import argparse
from encode_lib_common import (
    strip_ext_bigwig, ls_l, log, mkdir_p)
import numpy as np
import pandas as pd
import scipy.stats
import matplotlib as mpl
mpl.use('Agg')

warnings.filterwarnings("ignore")

def parse_arguments():
    parser = argparse.ArgumentParser(prog='ENCODE compare signal to roadmap')
    parser.add_argument('--bigwig', type=str,
                        help='BIGWIG file (from task macs2).')
    parser.add_argument('--dnase', type=str, help='DNase file.')
    parser.add_argument('--reg2map', type=str, help='Reg2map file.')
    parser.add_argument('--reg2map-bed', type=str, help='Reg2map bed file.')
    parser.add_argument('--roadmap-meta', type=str,
                        help='Roadmap metadata file.')
    parser.add_argument('--out-dir', default='', type=str,
                        help='Output directory.')
    parser.add_argument('--log-level', default='INFO', help='Log level',
                        choices=['NOTSET', 'DEBUG', 'INFO', 'WARNING',
                                 'CRITICAL', 'ERROR', 'CRITICAL'])
    args = parser.parse_args()
    log.setLevel(args.log_level)
    log.info(sys.argv)
    return args

def compare_to_roadmap(bw_file, regions_file, reg2map_file,
                       metadata, output_prefix):
    '''
    Takes a bigwig file and signal file, gets the bwAverageOverBed,
    then compares that signal with the signal in the Roadmap
    regions
    '''

    out_file = '{0}.signal'.format(output_prefix)
    log_file = '{0}.roadmap_compare.log'.format(output_prefix)

    # First get the signal vals for the peak regions
    # remember to use a UCSC formatted bed file for regions
    bw_average_over_bed = 'bigWigAverageOverBed {0} {1} {2}'.format(
        bw_file, regions_file, out_file)
    log.info(bw_average_over_bed)
    os.system(bw_average_over_bed)

    # Read the file back in
    sample_data = pd.read_table(out_file, header=None)
    sample_mean0_col = np.array(sample_data.iloc[:, 5])

    # Then, calculate correlations with all other Roadmap samples and rank
    # the correlations
    roadmap_signals = pd.read_table(reg2map_file, compression='gzip')
    (nrow, ncol) = roadmap_signals.shape

    results = pd.DataFrame(columns=('eid', 'corr'))
    with open(log_file, 'w') as fp:
        for i in range(ncol):
            # Slice, run correlation
            roadmap_i = roadmap_signals.iloc[:, i]
            spearman_corr = scipy.stats.spearmanr(np.array(roadmap_i),
                                                  sample_mean0_col)
            results.loc[i] = [roadmap_i.name, spearman_corr[0]]
            s = '{0}\t{1}'.format(roadmap_i.name, spearman_corr)
            log.info(s)
            fp.write(s + '\n')

    # Read in metadata to make the chart more understandable
    metadata = pd.read_table(metadata)
    metadata.columns = ['eid', 'mnemonic']

    merged = pd.merge(metadata, results, on='eid')

    sorted_results = merged.sort_values('corr', ascending=True)

    # Plot results
    pos = np.array(range(ncol)) + 0.5
    fig = plt.figure(figsize=(5, int(ncol/4)))
    plt.barh(pos, sorted_results['corr'], align='center', height=1.0)
    plt.yticks(pos, sorted_results['mnemonic'].tolist(), fontsize=7)
    plt.xlabel('Spearmans correlation')
    plt.title('Signal correlation to Roadmap DNase')
    plt.axis('tight')
    ax = plt.axes()
    ax.yaxis.set_ticks_position('none')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)

    plot_img = output_prefix + '.roadmap_compare_plot.png'
    fig.savefig(plot_img, format='png', bbox_inches='tight')

    return plot_img

def main():
    # read params
    args = parse_arguments()
    BIGWIG = args.bigwig
    DNASE = args.dnase
    OUTPUT_PREFIX = os.path.join(
        args.out_dir,
        os.path.basename(strip_ext_bigwig(BIGWIG)))

    REG2MAP_BED = args.reg2map_bed if args.reg2map_bed and os.path.basename(
        args.reg2map_bed) != 'null' else DNASE
    REG2MAP = args.reg2map if args.reg2map and os.path.basename(
        args.reg2map) != 'null' else ''
    ROADMAP_META = args.roadmap_meta if args.roadmap_meta and os.path.basename(
        args.roadmap_meta) != 'null' else ''

    log.info('Initializing and making output directory...')
    mkdir_p(args.out_dir)

    compare_to_roadmap(BIGWIG, REG2MAP_BED, REG2MAP,
                       ROADMAP_META, OUTPUT_PREFIX)

    log.info('List all files in output directory...')
    ls_l(args.out_dir)

    log.info('All done.')

if __name__ == '__main__':
    main()

reason: import numpy,import scipy


text:
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, PathPatch
from matplotlib.text import TextPath
from matplotlib.transforms import Affine2D
# This import registers the 3D projection, but is otherwise unused.
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 unused import
import mpl_toolkits.mplot3d.art3d as art3d

import numpy as np
xx, yy = np.meshgrid(np.linspace(0,15,16), np.linspace(0,15,16))

import cv2
img = cv2.imread('mitosisSmall2.png')
img = cv2.resize(img,dsize=(16,16))
img=img[:,:,1]
X =  xx
Y =  yy
Z =  0*np.ones(X.shape)

data = img
fig = plt.figure(figsize=(15,5))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, color='black', linewidth=1, facecolors=plt.cm.Greys(255-data), shade=False)
ax.grid(True)
for x in np.arange(0,15):
    for y in np.arange(0,15):
        print(x,y)
        ax.text(x+0.5,y+0.5,0.1,'%.1f' % (data[int(x),int(y)]/255), fontweight='bold', fontsize=4, color=[1.0,0.3,1.0], zorder=5, horizontalalignment='center', verticalalignment='center')

#ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, facecolors=plt.cm.Greys(255-data), shade=False)
ax.view_init(elev=90., azim=84)
ax.set_zlim(0,5)
# draw grid
#for y in
#ax.imshow(img)

xx, yy = np.meshgrid(np.arange(0,5), np.arange(0,5))
filtr = np.array([[255,255, 255,255,255], [0,255, 255,255,255], [0,0,255,255,255], [0,0,0,255,255], [0,0,0,0,255]])
print(xx.shape,yy.shape,filtr.shape)
ax.set_axis_off()
#ax.plot_surface(xx-8, yy+18, 2*np.ones(xx.shape), linewidth=1, zorder=10, facecolors=plt.cm.Greys(filtr), shade=False)
xo_act = -8
yo_act = 12

for x in np.arange(0,5):
    for y in np.arange(0,5):
        print(x,y)
        ax.text(x+xo_act+0.5,y+yo_act+0.5,2.1,'%.1f' % (filtr[int(x),int(y)]/255), fontweight='bold', fontsize=4, color='blue',horizontalalignment='center', verticalalignment='center')
#plt.savefig('conv2d.svg')

# draw grid
for k in range(16):
    plt.plot(xs=[0,15],ys=[k,k], zs=[0,0],color='k', linewidth=0.2)
    plt.plot(xs=[k,k],ys=[0,15], zs=[0,0],color='k', linewidth=0.2)

for k in [0,5]:
    startx=7
    starty=5
    plt.plot(xs=[startx,startx+5],ys=[starty+k,starty+k], zs=[0.01,0.01],color=[1.0,0.3,1.0], linewidth=1.)
    plt.plot(xs=[startx+k,startx+k],ys=[starty+0,starty+5], zs=[0.01,0.01],color=[1.0,0.3,1.0], linewidth=1.)

for k in range(6):
    plt.plot(xs=[xo_act+0,xo_act+5],ys=[yo_act+k,yo_act+k], zs=[0,0],color='k', linewidth=0.2)
    plt.plot(xs=[xo_act+k,xo_act+k],ys=[yo_act+0,yo_act+5], zs=[0,0],color='k', linewidth=0.2)

plt.plot(xs=[startx+5,xo_act+5], ys=[starty, yo_act],color=[1.0,0.3,1.0], linewidth=0.4, linestyle='--')
plt.plot(xs=[startx,xo_act], ys=[starty, yo_act],color=[1.0,0.3,1.0], linewidth=0.4, linestyle='--')
plt.plot(xs=[startx+5,xo_act+5], ys=[starty+5, yo_act+5],color=[1.0,0.3,1.0], linewidth=0.4, linestyle='--')

convresult = np.mean(img[starty:starty+5, startx:startx+5] *  filtr / 255 / 255)
print('Conv result:',convresult)
xo_act2 = -25
yo_act2 = 12
ax.text(startx+xo_act2+0.5,starty+yo_act2+0.5, z=0.0, s='%.1f' % convresult,fontweight='bold', fontsize=4, color='blue',horizontalalignment='center', verticalalignment='center')

#plt.plot(xs=[xo_act2+1+startx,xo_act+5], ys=[starty+yo_act2, yo_act],color=[1.0,0.3,1.0], linewidth=0.4, linestyle='--')
plt.plot(xs=[xo_act2+startx,xo_act], ys=[starty+yo_act2, yo_act],color='b', linewidth=0.4, linestyle='--')
plt.plot(xs=[xo_act2+startx+1,xo_act+5], ys=[starty+yo_act2+1, yo_act+5],color='b', linewidth=0.4, linestyle='--')

for k in [0,5]:
    plt.plot(xs=[xo_act,xo_act+5],ys=[yo_act+k,yo_act+k], zs=[0.01,0.01],color='b', linewidth=1.)
    plt.plot(xs=[xo_act+k,xo_act+k],ys=[yo_act+0,yo_act+5], zs=[0.01,0.01],color='b', linewidth=1.)

xo_act3=xo_act2+startx
yo_act3=yo_act2+starty
for k in [0,1]:
    plt.plot(xs=[xo_act3,xo_act3+1],ys=[yo_act3+k,yo_act3+k], zs=[0.01,0.01],color='b', linewidth=1.)
    plt.plot(xs=[xo_act3+k,xo_act3+k],ys=[yo_act3+0,yo_act3+1], zs=[0.01,0.01],color='b', linewidth=1.)

for k in range(12):
    plt.plot(xs=[xo_act2+0,xo_act2+11],ys=[yo_act2+k,yo_act2+k], zs=[0,0],color='k', linewidth=0.2)
    plt.plot(xs=[xo_act2+k,xo_act2+k],ys=[yo_act2+0,yo_act2+11], zs=[0,0],color='k', linewidth=0.2)
plt.tight_layout()
plt.savefig('conv2d.svg')

max_stars_repo_name: maubreville/Thesis

max_issues_repo_name: maubreville/Thesis

max_forks_repo_name: maubreville/Thesis