text:
import warnings
from typing import List, Dict

import numpy as np
import pandas as pd

import reeds.function_libs.visualization.sampling_plots

def undersampling_occurence_potential_threshold_densityClustering(ene_traj_csvs: List[pd.DataFrame],
                                                                  max_distance_kJ: float = 300,
                                                                  sampling_fraction_treshold: float = 0.9)->List[float]:
    """
    This function is estimating the pot_tresh for all states by using DBSCAN identifying the density region containing sampling_fraction_treshold of the data.
    The mean and std of the density region will result in pot_tresh = mean+3std
    This function is rather performance expensive.
    ! WARNING REQUIRES SKLEARN -Not Used by the pipeline. - Prototype!

    Parameters
    ----------
    ene_traj_csvs : List[pd.Dataframe]
        a list of pandas dataframes containing the energy data of each state eX
    max_distance_kJ : float, optional
        the maximum distance between two points for dbscan clustering (eps) (default:300kJ).
    sampling_fraction_treshold : float, optional
        the sampling threshold for detecting undersampling (default: 0.9).

    Returns
    -------
    List[float]
        list of individual identified potential thresholds
    """
    try:
        from sklearn.cluster import DBSCAN
    except:
        raise Exception("Please provide package sklearn!")

    pot_tresh_pre_rep = []

    state_names = [i for i in ene_traj_csvs[0].columns if (i.startswith("e") and not i == "eR")]
    num_s_values = len(ene_traj_csvs)
    num_states = len(state_names)

    total_step_num = ene_traj_csvs[0].shape[0]
    threshold_fulfilled = [True for j in range(num_s_values)]

    v_pot_sampling = [[0 for y in state_names] for x in range(num_s_values)]

    for i, replica_data in enumerate(ene_traj_csvs):
        # print("")
        # print("replica " + str(i + 1) + ":")
        pot_tresh = []

        # calculate for each state the pot thresh, depending on clustering.
        for k, state in enumerate(state_names):
            vec = np.array(replica_data[state]).reshape(-1, 1)
            db = DBSCAN(eps=max_distance_kJ, min_samples=total_step_num * sampling_fraction_treshold).fit(vec)
            v_pot_sampling[i][k] = len(db.components_) / total_step_num

            if (not len(db.components_)):
                threshold_fulfilled[i] = False
            else:
                # print(k)
                threshold = np.mean(db.components_) + 3 * np.std(db.components_)
                pot_tresh.append(threshold)
        if (len(pot_tresh) == len(state_names)):
            pot_tresh_pre_rep.append(pot_tresh)

    # final averaging
    pot_thresh_per_state_and_replica = np.array([t for t in pot_tresh_pre_rep if (len(t) == num_states)]).T
    pot_thresh_per_state = np.mean(pot_thresh_per_state_and_replica, axis=1) + np.std(pot_thresh_per_state_and_replica,
                                                                                axis=1)
    return pot_thresh_per_state

def undersampling_occurence_potential_threshold_distribution_based(ene_traj_csvs: List[pd.DataFrame],
                                                                   sampling_fraction_treshold: float = 0.9,
                                                                    verbose:bool = False)->List[float]:
    """undersampling_occurence_potential_threshold_distribution_based
    This function is estimating the pot_tresh for all states by testing if around the minimal energy 90% of the data is located in a threshold of  max_distance_kJ.
    The mean and std of the density region will result in pot_tresh = min(V)+6*std(V
    This function has low demand in performance.

    Parameters
    ----------
    ene_traj_csvs : List[pd.Dataframe]
        a list of pandas dataframes containing the energy data of each state eX
    sampling_fraction_treshold : float, optional
        the sampling threshold for detecting undersampling (default: 0.9).

    Returns
    -------
    List[float]
        list of individual identified potential thresholds
    """
    pot_tresh_pre_rep = []
    state_names = [i for i in ene_traj_csvs[0].columns if (i.startswith("e") and not i == "eR")]
    num_states = len(state_names)
    total_step_num = ene_traj_csvs[0].shape[0]

    std = [6*np.std(ene_traj_csvs[-1][state]) for state in state_names]
    #mean_potential = [np.mean(ene_traj_csvs[-1][state]) for state in state_names]
    for i, replica_data in enumerate(ene_traj_csvs):
        if(verbose): print("replica " + str(i + 1) + ":")
        mean_potential = [np.min(replica_data[state]) for state in state_names]
        thresholds = [m+s for m, s in zip(mean_potential, std)]
        if(verbose): print("Checkout std: ", thresholds)
        # iterate over states, identify all undersampling replicas, calc fraction
        pot_tresh = []
        for k, state in enumerate(state_names):
            vec = np.array(replica_data[state]).reshape(-1, 1)
            treshold =  thresholds[k]

            below_thresh = vec[vec < treshold]
            below_thresh_fraction = len(below_thresh) / total_step_num
            if(verbose): print("\t belowTresh: ", state, sampling_fraction_treshold < below_thresh_fraction)
            if (sampling_fraction_treshold < below_thresh_fraction):
                threshold = np.mean(below_thresh) + 3 * np.std(below_thresh)
                pot_tresh.append(threshold)

        # only save undersampling pot_threshes
        if (len(pot_tresh) == num_states):
            pot_tresh_pre_rep.append(pot_tresh)

    # final averaging.
    pot_thresh_per_state_and_replica = np.array([t for t in pot_tresh_pre_rep if (len(t) == num_states)]).T

    print(pot_thresh_per_state_and_replica)
    if (len(pot_thresh_per_state_and_replica) > 0):
        pot_thresh_per_state = np.min(pot_thresh_per_state_and_replica, axis=1) + np.std(pot_thresh_per_state_and_replica, axis=1)
    else:
        pot_thresh_per_state = [0 for x in range(num_states)]
    return pot_thresh_per_state

def physical_occurence_potential_threshold_distribution_based(ene_traj_csv: pd.DataFrame, equilibrate_dominationState:float=0.01, _vacuum_simulation:bool=False,  verbose:bool=False)->List[float]:

    """physical_occurence_potential_threshold_distribution_based
    This function is estimating the pot_tresh for all states by testing if around the minimal energy 90% of the data is located in a threshold of  max_distance_kJ.
    The mean and std of the density region will result in pot_tresh = mean+3*std
    This function is cheap in performance.

    Parameters
    ----------
    ene_traj_csv : pd.Dataframe
        a pandas dataframe containing the energy data of each state eX
    equilibrate_dominationState : int, optional
        equilibrate the domination state for this fraction
    _vacuum_simulation=vacuum_simulation: bool, optional
        this flag switches the algorithm to use all potential energies of a state, instead of only the domination sampling.
    verbose: bool, optional
        print fun
    Returns
    -------
    List[float]
        list of individual identified potential thresholds
    """

    state_names = [i for i in ene_traj_csv.columns if (i.startswith("e") and not i == "eR")]
    total_step_num = ene_traj_csv.shape[0]
    data =ene_traj_csv

    # iterate over states, identify all undersampling replicas, calc fraction
    pot_tresh = []
    for k, state in enumerate(state_names):
        #get only the sampling, in which the target state is dominating: pre-filter noise!
        if(_vacuum_simulation):
            state_domination_sampling = data
        else:
            state_domination_sampling = data.where(data[state_names].min(axis=1) == data[state]).dropna()
        start_after_eq = int(np.round(equilibrate_dominationState*state_domination_sampling.shape[0]))
        state_domination_sampling = state_domination_sampling.iloc[start_after_eq:]

        if(verbose):
            below_thresh_fraction = state_domination_sampling.shape[0] / total_step_num
            print("State "+str(state)+" - Domination state sampling fraction ", below_thresh_fraction)

        #calculate threshold
        threshold = np.min(state_domination_sampling[state]) + 6 * np.std(state_domination_sampling[state])
        pot_tresh.append(threshold)

    return pot_tresh

def get_all_physical_occurence_potential_threshold_distribution_based(ene_trajs: List[pd.DataFrame], _vacuum_simulation:bool=False,)->List[float]:
    """
        This function is used in the state optimization approach and gives the physical potential energy threshold for occurrence sampling back, that is estimated from the optimized EDS-System.

    Parameters
    ----------
    ene_trajs: List[pd.DataFrame]
        energy trajectories from all end state optimizations
    _vacuum_simulation=vacuum_simulation: bool, optional
        this flag switches the algorithm to use all potential energies of a state, instead of only the domination sampling.
    Returns
    -------
    List[float]
        list of potential energy thresholds.
    """
    opt_pot_tresh = []
    for key, traj in enumerate(ene_trajs):
        print(key)
        pot_tresh_state = physical_occurence_potential_threshold_distribution_based(traj, _vacuum_simulation=_vacuum_simulation)

        print(pot_tresh_state)
        if (np.isnan(pot_tresh_state[key])):
            warnings.warn("A state potential threshold was NaN -> This hints on that you did not sample the state as a dominating one! Please check your simulatigons!")
        opt_pot_tresh.append(pot_tresh_state[key])
    return opt_pot_tresh

def calculate_sampling_distributions(ene_traj_csvs: List[pd.DataFrame],
                                     potential_treshold: List[float])-> Dict[int, Dict[str, Dict[int, float]]]:
    """calculate_sampling_distributions
    This function is using the dominating state sampling and occurrence state sampling definition, to calculate
    for both definitions the sampling distributions including each stat.
    Further each replica, having an occurrence sampling >= undersampling_occurence_sampling_tresh for each state is categorized as undersampling.

    Parameters
    ----------
    ene_traj_csvs: List[pd.Dataframe]
        a list of pandas dataframes containing the energy data of each state eX
    potential_treshold: List[float]
        a list of potential thresholds for undersampling
    undersampling_occurence_sampling_tresh: float, optional
        threshold for the fraction of the energies which has to be below the threshold (default 0.75)

    Returns
    -------
    Dict[int, Dict[str, Dict[int, float]]]
        the dictionary contains for all replicas, the occurrence/dominating sampling fractions of each state  and undersamplin classification
    """

    replica_sampling_dist = {}
    num_states = sum([1 for i in ene_traj_csvs[0].columns if (i.startswith("e") and not i == "eR")])
    states = ["e" + str(s) for s in range(1, 1 + num_states)]
    for replica in ene_traj_csvs:
        total_number_steps = replica.shape[0]

        # Domination sampling
        dominating_state_sampling = {x: 0 for x in range(1, 1 + num_states)}
        min_state, counts = np.unique(replica[states].idxmin(axis=1), return_counts=True)
        dominating_state_sampling.update(
            {int(key.replace("e", "")): val / total_number_steps for key, val in zip(min_state, counts)})

        # occurence samplng
        occurrence_state_sampling = {}
        for ind, state in enumerate(states):
            occurrence_sampling_frac = replica[replica[state] < potential_treshold[ind]].shape[0] / total_number_steps
            occurrence_state_sampling.update({int(state.replace("e", "")): occurrence_sampling_frac})

        # update results
        replica_sampling_dist.update({int(replica.s.replace("s", "")): {"dominating_state": dominating_state_sampling,
                                                                        "occurence_state": occurrence_state_sampling}})

    return replica_sampling_dist

def sampling_analysis(ene_traj_csvs: List[pd.DataFrame],
                      state_potential_treshold: List[float],
                      s_values: List[float],
                      out_path: str = None,
                      xmax: bool = False,
                      _visualize: bool = True,
                      verbose: bool = False, _usample_run:bool=False) -> (dict, str):
    """sampling_analysis
    This function is analysing the samplings

    Parameters
    ----------
    out_path :  str
        path out for the plotfiles
    ene_traj_csvs : List[pd.DataFrame]
        contains the energy data
    s_values :  List[float]
        list of s_values
    state_potential_treshold : List[float]
        potential energy thresholds, for considering a state as sampled
    xmax :  bool, optional
        additionally output a plot only with a smaller x_range. (0->xmax) (default False)
    _visualize: bool, opptional
        create additional plots (default True)
    verbose: bool, optional
        story time :) (default False)

    Returns
    -------
    str
        out_path
    dict
      returns a dictionary containing the information of undersampling start, undersampling_potential thresholds and the different sampling type fractions. (physical sampling occurence not included!)
    """

    # read all Vy_sx_files!
    from reeds.function_libs.visualization.utils import nice_s_vals
    if (verbose): print("\n\n Potential Threshold\n\n")

    ##glob vars
    num_states = len(state_potential_treshold)
    select_states = ["e" + str(x) for x in range(1, num_states + 1)]
    s_vals_nice = nice_s_vals(s_values)

    # show_presence of undersampling
    if (verbose): print("\n\n Sampling Timeseries\n\n")
    for ind, replica in enumerate(ene_traj_csvs):
        if (verbose): print("\t replica " + replica.s)

        dominating_state_sampling = replica[select_states].idxmin(axis=1).replace("e", "", regex=True)

        occurrence_sampling_replica = []
        for state in select_states:
            occurrence_sampling_state_replica = replica.index[
                replica[state] < state_potential_treshold[int(state.replace("e", "")) - 1]]
            occurrence_sampling_replica.append(occurrence_sampling_state_replica)

        data = {"time": replica.time, "occurrence_t": occurrence_sampling_replica, "dominating_state": dominating_state_sampling}

        if (_visualize):
            reeds.function_libs.visualization.sampling_plots.plot_t_statepres(data=data,
                                                                              title="s=" + str(s_vals_nice[ind]),
                                                                              out_path=out_path + "/sampling_timeseries_s" + str(ind + 1) + ".png")

        if (xmax):
            reeds.function_libs.visualization.sampling_plots.plot_t_statepres(data=data,
                                                                              title="s=" + str(s_vals_nice[ind]),
                                                                              out_path=out_path + "/sampling_timeseries_s" + str(ind + 1) + "_upto_" + str(
                                     xmax) + ".png",
                                                                              xlim=[0, xmax])

    # SamplingMatrix by kays
    if (verbose): print("\n\n Calculate Sampling Distributions\n\n")
    replica_sampling_distributions = calculate_sampling_distributions(ene_traj_csvs=ene_traj_csvs,
                                                                      potential_treshold=state_potential_treshold)
    if (_visualize):
        if (verbose): print("\n\n Sampling Histograms\n\n")
        for ind, x in enumerate(replica_sampling_distributions):
            reeds.function_libs.visualization.sampling_plots.plot_stateOccurence_hist(data=replica_sampling_distributions[x],
                                                                                title="s=" + str(s_vals_nice[ind]),
                                                                                out_path=out_path + "/sampling_hist_" + str(x) + ".png")
    if (_visualize and not _usample_run):
        if (verbose): print("\n\n Sampling Matrix\n\n")
        reeds.function_libs.visualization.sampling_plots.plot_stateOccurence_matrix(data=replica_sampling_distributions, out_dir=out_path, s_values=s_vals_nice,
                                                                                place_undersampling_threshold=False, title_suffix="")

    final_results = {"potentialThreshold": state_potential_treshold,
                     "samplingDistributions": replica_sampling_distributions,
                    }

    return final_results, out_path

def detect_undersampling(ene_traj_csvs: List[pd.DataFrame],
                         state_potential_treshold: List[float],
                         s_values: List[float],
                         out_path: str = None,
                         undersampling_occurence_sampling_tresh: float = 0.9,
                         xmax: bool = False,
                         _visualize: bool = True,
                         verbose: bool = False) -> (dict, str):
    """sampling_analysis
    This function is analysing the samplings

    Parameters
    ----------
    out_path :  str
        path out for the plotfiles
    ene_traj_csvs : List[pd.DataFrame]
        contains the energy data
    s_values :  List[float]
        list of s_values
    state_potential_treshold : List[float]
        potential energy thresholds, for considering a state as sampled
    xmax :  bool, optional
        additionally output a plot only with a smaller x_range. (0->xmax) (default False)
    _visualize: bool, opptional
        create additional plots (default True)
    verbose: bool, optional
        story time :) (default False)

    Returns
    -------
    str
        out_path
    dict
      returns a dictionary containing the information of undersampling start, undersampling_potential thresholds and the different sampling type fractions. (physical sampling occurence not included!)
    """

    # read all Vy_sx_files!
    from reeds.function_libs.visualization.utils import nice_s_vals
    if (verbose): print("\n\n Potential Threshold\n\n")

    ##glob vars
    sampling_stat, out_path = sampling_analysis(ene_traj_csvs=ene_traj_csvs,
                                                state_potential_treshold=state_potential_treshold,
                                                s_values=s_values,
                                                out_path=out_path,
                                                xmax= xmax,
                                                _visualize=_visualize,
                                                verbose = verbose, _usample_run=True)

    ##get undersampling id:
    found_undersampling = False
    undersampling_idx = None
    replica_sampling_distributions = sampling_stat["samplingDistributions"]
    for i in replica_sampling_distributions:
        # undersampling?
        undersampling_criterium = True if (all([occ >= undersampling_occurence_sampling_tresh for occ in replica_sampling_distributions[i]["occurence_state"].values()])) else False
        sampling_stat["samplingDistributions"][i].update({"undersampling":undersampling_criterium})

        if (undersampling_criterium):
            undersampling_idx = i
            break

    if (not found_undersampling):
        warnings.warn("Could not find undersampling!")

    if (_visualize):
        if (verbose): print("\n\n Sampling Matrix\n\n")
        s_vals_nice = nice_s_vals(s_values)
        if(_visualize):
            reeds.function_libs.visualization.sampling_plots.plot_stateOccurence_matrix(data=replica_sampling_distributions, out_dir=out_path, s_values=s_vals_nice,
                                                                                place_undersampling_threshold=True, title_suffix="")

    sampling_stat.update({"undersamplingThreshold": undersampling_idx})
    sampling_stat.update({"state_undersampling_potTresh": state_potential_treshold,
                          "undersampling_occurence_sampling_tresh": undersampling_occurence_sampling_tresh})

    return sampling_stat, out_path

max_stars_repo_name: candidechamp/reeds

max_issues_repo_name: candidechamp/reeds

max_forks_repo_name: candidechamp/reeds


text:
#!/usr/bin/python3
"""
Polyhedral set library.

This library implements convex regions of the form H x <= k, where H, x, and k
are matricies.  It also provides convenient methods to find all the verticies.
"""

__author__ = 'Austin Schuh (austin.linux@gmail.com)'

from frc971.control_loops.python import libcdd
import numpy
import string
import sys

def _PiecewiseConcat(*args):
    """Concatenates strings inside lists, elementwise.

  Given ['a', 's'] and ['d', 'f'], returns ['ad', 'sf']
  """
    return map(''.join, zip(*args))

def _SplitAndPad(s):
    """Splits a string on newlines, and pads the lines to be the same width."""
    split_string = s.split('\n')
    width = max(len(stringpiece) for stringpiece in split_string) + 1

    padded_strings = [
        stringpiece.ljust(width, ' ') for stringpiece in split_string
    ]
    return padded_strings

def _PadHeight(padded_array, min_height):
    """Adds lines of spaces to the top and bottom of an array symmetrically."""
    height = len(padded_array)
    if height < min_height:
        pad_array = [' ' * len(padded_array[0])]
        height_error = min_height - height
        return (pad_array * int(
            (height_error) / 2) + padded_array + pad_array * int(
                (height_error + 1) / 2))
    return padded_array

class HPolytope(object):
    """This object represents a H-polytope.

  Polytopes are convex regions in n-dimensional space.
  For H-polytopes, this is represented as the intersection of a set of half
  planes.  The mathematic equation that represents this is H x <= k.
  """

    def __init__(self, H, k):
        """Constructs a H-polytope from the H and k matricies.

    Args:
      H: numpy.Matrix (n by k), where n is the number of constraints, and k is
        the number of dimensions in the space.  Does not copy the matrix.
      k: numpy.Matrix (n by 1).  Does not copy the matrix.
    """
        self._H = H
        self._k = k

    @property
    def k(self):
        """Returns the k in H x <= k."""
        return self._k

    @property
    def H(self):
        """Returns the H in H x <= k."""
        return self._H

    @property
    def ndim(self):
        """Returns the dimension of the set."""
        return self._H.shape[1]

    @property
    def num_constraints(self):
        """Returns the number of constraints defining the set."""
        return self._k.shape[0]

    def IsInside(self, point):
        """Returns true if the point is inside the polytope, edges included."""
        return (self._H * point <= self._k).all()

    def Vertices(self):
        """Returns a matrix with the vertices of the set in its rows."""
        # TODO(aschuh): It would be better to write some small C/C++ function that
        # does all of this and takes in a numpy array.
        # The copy is expensive in Python and cheaper in C.

        # Create an empty matrix with the correct size.
        matrixptr = libcdd.dd_CreateMatrix(self.num_constraints, self.ndim + 1)
        matrix = matrixptr.contents

        try:
            # Copy the data into the matrix.
            for i in range(self.num_constraints):
                libcdd.dd_set_d(matrix.matrix[i][0], self._k[i, 0])
                for j in range(self.ndim):
                    libcdd.dd_set_d(matrix.matrix[i][j + 1], -self._H[i, j])

            # Set enums to the correct values.
            matrix.representation = libcdd.DD_INEQUALITY
            matrix.numbtype = libcdd.DD_REAL

            # TODO(aschuh): Set linearity if it is useful.
            # This would be useful if we had any constraints saying B - A x = 0

            # Build a Polyhedra
            polyhedraptr = libcdd.dd_DDMatrix2Poly(matrixptr)

            if not polyhedraptr:
                return None

            try:
                # Return None on error.
                # The error values are enums, so they aren't exposed.

                # Magic happens here.  Computes the vertices
                vertex_matrixptr = libcdd.dd_CopyGenerators(polyhedraptr)
                vertex_matrix = vertex_matrixptr.contents

                try:
                    # Count the number of vertices and rays in the result.
                    num_vertices = 0
                    num_rays = 0
                    for i in range(vertex_matrix.rowsize):
                        if libcdd.dd_get_d(vertex_matrix.matrix[i][0]) == 0:
                            num_rays += 1
                        else:
                            num_vertices += 1

                    # Build zeroed matricies for the new vertices and rays.
                    vertices = numpy.matrix(
                        numpy.zeros((num_vertices, vertex_matrix.colsize - 1)))
                    rays = numpy.matrix(
                        numpy.zeros((num_rays, vertex_matrix.colsize - 1)))

                    ray_index = 0
                    vertex_index = 0

                    # Copy the data out of the matrix.
                    for index in range(vertex_matrix.rowsize):
                        if libcdd.dd_get_d(
                                vertex_matrix.matrix[index][0]) == 0.0:
                            for j in range(vertex_matrix.colsize - 1):
                                rays[ray_index, j] = libcdd.dd_get_d(
                                    vertex_matrix.matrix[index][j + 1])
                            ray_index += 1
                        else:
                            for j in range(vertex_matrix.colsize - 1):
                                vertices[vertex_index, j] = libcdd.dd_get_d(
                                    vertex_matrix.matrix[index][j + 1])
                            vertex_index += 1
                finally:
                    # Free everything.
                    libcdd.dd_FreeMatrix(vertex_matrixptr)

            finally:
                libcdd.dd_FreePolyhedra(polyhedraptr)

        finally:
            libcdd.dd_FreeMatrix(matrixptr)

        # Rays are unsupported right now.  This may change in the future.
        assert (rays.shape[0] == 0)

        return vertices

    def __str__(self):
        """Returns a formatted version of the polytope.

    The dump will look something like the following, which prints out the matrix
    comparison.

    [[ 1  0]            [[12]
     [-1  0]  [[x0]  <=  [12]
     [ 0  1]   [x1]]     [12]
     [ 0 -1]]            [12]]
    """
        height = max(self.ndim, self.num_constraints)

        # Split the print up into 4 parts and concatenate them all.
        H_strings = _PadHeight(_SplitAndPad(str(self.H)), height)
        v_strings = _PadHeight(_SplitAndPad(str(self.k)), height)
        x_strings = _PadHeight(self._MakeXStrings(), height)
        cmp_strings = self._MakeCmpStrings(height)

        return '\n'.join(
            _PiecewiseConcat(H_strings, x_strings, cmp_strings, v_strings))

    def _MakeXStrings(self):
        """Builds an array of strings with constraint names in it for printing."""
        x_strings = []
        if self.ndim == 1:
            x_strings = ["[[x0]] "]
        else:
            for index in range(self.ndim):
                if index == 0:
                    x = "[[x%d]  " % index
                elif index == self.ndim - 1:
                    x = " [x%d]] " % index
                else:
                    x = " [x%d]  " % index
                x_strings.append(x)
        return x_strings

    def _MakeCmpStrings(self, height):
        """Builds an array of strings with the comparison in it for printing."""
        cmp_strings = []
        for index in range(height):
            if index == int((height - 1) / 2):
                cmp_strings.append("<= ")
            else:
                cmp_strings.append("   ")
        return cmp_strings

ext: py

reason: import numpy


text:
import torch
import numpy as np
import pickle
import logging
import platalea.basic as basic
import platalea.encoders as encoders
import platalea.dataset as dataset
import platalea.config
import json
import os

import configargparse

parser = configargparse.get_argument_parser('platalea')

config_args, unknown_args = parser.parse_known_args()

_device = platalea.config.device()

def prepare():
    logging.getLogger().setLevel('INFO')
    logging.info("Loading pytorch models")
    net_rand = basic.SpeechImage(basic.DEFAULT_CONFIG).to(_device)
    net_rand.eval()
    net_train = basic.SpeechImage(basic.DEFAULT_CONFIG)
    net_train.load_state_dict(torch.load("net.20.pt").state_dict())
    net_train.to(_device)
    net_train.eval()
    nets = [('trained', net_train), ('random', net_rand)]
    with torch.no_grad():
        save_data(nets, directory='.', batch_size=32)

def massage_transformer_data():
    logging.getLogger().setLevel('INFO')
    factors = pickle.load(open("data/in/trans/downsampling_factors.pkl", "rb"))
    layers = factors.keys()
    logging.info("Loading input")
    data = pickle.load(open("data/in/trans/global_input.pkl", "rb"))
    logging.info("Fixing IDs")
    data['audio_id'] = np.array([ i + '.wav' for i in data['audio_id']])
    logging.info("Adding IPA")
    alignment = load_alignment("data/out/trans/fa.json")
    data['ipa'] = np.array([ align2ipa(alignment[i]) for i in data['audio_id'] ])
    logging.info("Saving input")
    pickle.dump(data, open("data/out/trans/global_input.pkl", "wb"), protocol=4)
    for mode in ['trained', 'random']:
        for layer in layers:
            logging.info("Loading data for {} {}".format(mode, layer))
            data = pickle.load(open("data/in/trans/global_{}.{}.pkl".format(mode, layer), "rb"))
            data = {layer: data}
            logging.info("Saving data for {} {}".format(mode, layer))
            pickle.dump(data, open("data/out/trans/global_{}_{}.pkl".format(mode, layer), "wb"), protocol=4)

def vgs_factors():
    return {'conv': {'pad': 0, 'ksize': 6, 'stride': 2 },
                   'rnn0': None,
                   'rnn1': None,
                   'rnn2': None,
                   'rnn3': None }

def massage_vgs_asr_data():
    logging.getLogger().setLevel('INFO')
    filter_global_data("data/in/vgs-asr/")
    factors = vgs_factors()
    layers = factors.keys()
    logging.info("Loading input")
    data = pickle.load(open("data/in/vgs-asr/global_input.pkl", "rb"))
    logging.info("Adding IPA")
    alignment = load_alignment("data/out/vgs-asr/fa.json")
    data['ipa'] = np.array([ align2ipa(alignment[i]) for i in data['audio_id'] ])
    logging.info("Saving input")
    pickle.dump(data, open("data/out/vgs-asr/global_input.pkl", "wb"), protocol=4)
    for mode in ['trained', 'random']:
        logging.info("Loading data for {}".format(mode))
        data = pickle.load(open("data/in/vgs-asr/global_{}.pkl".format(mode), "rb"))
        for layer in layers:
            ldata = {layer: data[layer]}
            logging.info("Saving data for {} {}".format(mode, layer))
            pickle.dump(data, open("data/out/vgs-asr/global_{}_{}.pkl".format(mode, layer), "wb"), protocol=4)

def save_data(nets, directory, batch_size=32):
    save_global_data(nets, directory=directory, batch_size=batch_size) # FIXME adapt this per directory too
    save_local_data(directory=directory)

def save_global_data(nets,  directory='.', batch_size=32):
    """Generate data for training a phoneme decoding model."""
    logging.info("Loading alignments")
    data = load_alignment("{}/fa.json".format(directory))
    logging.info("Loading audio features")
    val = dataset.Flickr8KData(root=config_args.data_root, split='val')
    #
    alignments = [ data[sent['audio_id']] for sent in val ]
    # Only consider cases where alignement does not fail
    alignments = [item for item in alignments if good_alignment(item) ]
    sentids = set(item['audio_id'] for item in alignments)
    audio = [ sent['audio'] for sent in val if sent['audio_id'] in sentids ]
    audio_np = [ a.numpy() for a in audio]

    ## Global data

    global_input = dict(audio_id = np.array([datum['audio_id'] for datum in alignments]),
                       ipa =      np.array([align2ipa(datum)  for datum in alignments]),
                       text =     np.array([datum['transcript'] for datum in alignments]),
                       audio = np.array(audio_np))
    pickle.dump(global_input, open("global_input.pkl", "wb"), protocol=4)

    for mode, net in nets:
        global_act = collect_activations(net, audio, batch_size=batch_size)
        for layer in global_act:
            logging.info("Saving global data in {}/global_{}_{}.pkl".format(directory, mode, layer))
            pickle.dump({layer: global_act[layer]}, open("{}/global_{}_{}.pkl".format(directory, mode, layer), "wb"), protocol=4)

def filter_global_data(directory):
    """Remove sentences with OOV items from data."""
    logging.getLogger().setLevel('INFO')
    logging.info("Loading raw data")
    global_input =  pickle.load(open("{}/global_input_raw.pkl".format(directory), "rb"))
    logging.info("Loading alignments")
    adata = load_alignment("{}/fa.json".format(directory))

    logging.info("Filtering out failed alignments and OOV")
    alignments = [ adata.get(i, adata.get(i+'.wav')) for i in global_input['audio_id'] ]
    # Only consider cases where alignement does not fail
    # Only consider cases with no OOV items
    alignments = [item for item in alignments if good_alignment(item) ]
    sentids = set(item['audio_id'] for item in alignments)
    ## Global data

    include = np.array([ i in sentids for i in global_input['audio_id'] ])

    #filtered = { key: np.array(value)[include] for key, value in global_input.items() }
    # Hack to fix broken format
    filtered = {}
    for key, value in global_input.items():
        if key == "audio":
            value = [ v.numpy() for v in value ]
        filtered[key] = np.array(value)[include]

    logging.info("Saving filtered data")
    pickle.dump(filtered, open("{}/global_input.pkl".format(directory), "wb"), protocol=4)

    for name in ['trained', 'random']:
        global_act = pickle.load(open("{}/global_{}_raw.pkl".format(directory, name), "rb"))
        filtered = { key: np.array(value)[include] for key, value in global_act.items() }
        logging.info("Saving filtered global data in global_{}.pkl".format(name))
        pickle.dump(filtered, open("{}/global_{}.pkl".format(directory, name), "wb"), protocol=4)

def good_alignment(item):
    for word in item['words']:
        if word['case'] != 'success' or word['alignedWord'] == '<unk>':
            return False
    return True

def make_indexer(factors, layer):
     def inout(pad, ksize, stride, L):
         return ((L + 2*pad - 1*(ksize-1) - 1) // stride + 1)
     def index(ms):
         t = ms//10
         for l in factors:
             if factors[l] is None:
                 pass
             else:
                 pad = factors[l]['pad']
                 ksize = factors[l]['ksize']
                 stride = factors[l]['stride']
                 t = inout(pad, ksize, stride, t)
             if l == layer:
                 break
         return t
     return index

def save_local_data(directory):
    logging.getLogger().setLevel('INFO')
    logging.info("Loading alignments")
    global_input =  pickle.load(open("{}/global_input.pkl".format(directory), "rb"))
    adata = load_alignment("{}/fa.json".format(directory))
    alignments = [ adata.get(i, adata.get(i+'.wav')) for i in global_input['audio_id'] ]
    ## Local data
    local_data = {}
    logging.info("Computing local data for MFCC")
    y, X = phoneme_activations(global_input['audio'], alignments, index=lambda ms: ms//10, framewise=True)
    local_input = check_nan(features=X, labels=y)
    pickle.dump(local_input, open("{}/local_input.pkl".format(directory), "wb"), protocol=4)
    try:
        factors = pickle.load(open("{}/downsampling_factors.pkl".format(directory), "rb"))
    except FileNotFoundError:
        # Default VGS settings
        factors = vgs_factors()
    for mode in ['trained', 'random']:
        for layer in factors.keys():
            if layer == "conv1":
                pass # This data is too big
            else:
                global_act = pickle.load(open("{}/global_{}_{}.pkl".format(directory, mode, layer), "rb"))
                local_act = {}
                index = make_indexer(factors, layer)
                logging.info("Computing local data for {}, {}".format(mode, layer))
                y, X = phoneme_activations(global_act[layer], alignments, index=index, framewise=True)
                local_act[layer] = check_nan(features=X, labels=y)
                logging.info("Saving local data in local_{}_{}.pkl".format(mode, layer))
                pickle.dump(local_act, open("{}/local_{}_{}.pkl".format(directory, mode, layer), "wb"), protocol=4)

def load_alignment(path):
    data = {}
    for line in open(path):
        item = json.loads(line)
        item['audio_id'] = os.path.basename(item['audiopath'])
        data[item['audio_id']] = item
    return data

def phoneme_activations(activations, alignments, index=lambda ms: ms//10, framewise=True):
    """Return array of phoneme labels and array of corresponding mean-pooled activation states."""
    labels = []
    states = []
    for activation, alignment in zip(activations, alignments):
        # extract phoneme labels and activations for current utterance
        if framewise:
            fr = list(frames(alignment, activation, index=index))
        else:
            fr = list(slices(alignment, activation, index=index))
        if len(fr) > 0:
            y, X = zip(*fr)
            y = np.array(y)
            X = np.stack(X)
            labels.append(y)
            states.append(X)
    return np.concatenate(labels), np.concatenate(states)

def align2ipa(datum):
    """Extract IPA transcription from alignment information for a sentence."""
    from platalea.analysis.ipa import arpa2ipa
    result = []
    for word in datum['words']:
        for phoneme in word['phones']:
            result.append(arpa2ipa(phoneme['phone'].split('_')[0], '_'))
    return ''.join(result)

def slices(utt, rep, index, aggregate=lambda x: x.mean(axis=0)):
    """Return sequence of slices associated with phoneme labels, given an
       alignment object `utt`, a representation array `rep`, and
       indexing function `index`, and an aggregating function\
       `aggregate`.
    """
    for phoneme in phones(utt):
        phone, start, end = phoneme
        assert index(start)<index(end)+1, "Something funny: {} {} {} {}".format(start, end, index(start), index(end))
        yield (phone, aggregate(rep[index(start):index(end)+1]))

def frames(utt, rep, index):
     """Return pair sequence of (phoneme label, frame), given an
        alignment object `utt`, a representation array `rep`, and
       indexing function `index`.
     """
     for phoneme in phones(utt):
         phone, start, end = phoneme
         assert index(start)<index(end)+1, "Something funny: {} {} {} {}".format(start, end, index(start), index(end))
         for j in range(index(start), index(end)+1):
             if j < rep.shape[0]:
                 yield (phone, rep[j])
             else:
                 logging.warning("Index out of bounds: {} {}".format(j, rep.shape))

def phones(utt):
    """Return sequence of phoneme labels associated with start and end
     time corresponding to the alignment JSON object `utt`.

    """
    for word in utt['words']:
        pos = word['start']
        for phone in word['phones']:
            start = pos
            end = pos + phone['duration']
            pos = end
            label = phone['phone'].split('_')[0]
            if label != 'oov':
                yield (label, int(start*1000), int(end*1000))

def check_nan(labels, features):
    # Get rid of NaNs
    ix = np.isnan(features.sum(axis=1))
    logging.info("Found {} NaNs".format(sum(ix)))
    X = features[~ix]
    y = labels[~ix]
    return dict(features=X, labels=y)

def collect_activations(net, audio, batch_size=32):
    data = torch.utils.data.DataLoader(dataset=audio,
                                       batch_size=batch_size,
                                       shuffle=False,
                                       num_workers=0,
                                       collate_fn=dataset.batch_audio)
    out = {}
    for au, l in data:
        act = net.SpeechEncoder.introspect(au.to(_device), l.to(_device))
        for k in act:
            if k not in out:
                out[k] = []
            out[k]  += [ item.detach().cpu().numpy() for item in act[k] ]
    return { k: np.array(v) for k,v in out.items() }

def spec(d):
    if type(d) != type(dict()):
        return type(d)
    else:
        return { key:spec(val) for key, val in d.items() }

max_stars_repo_name: gchrupala/platalea

max_issues_repo_name: gchrupala/platalea

max_forks_repo_name: gchrupala/platalea


text:
# -*- coding: utf-8 -*-
"""
Plot of Rayleigh curves in several diagrams
@author: j.gressier
"""

import numpy                 as np
import matplotlib.pyplot     as plt
import aero.Rayleigh         as ray

npoints = 100
gam     = 1.4

fig=plt.figure(1, figsize=(10,8))
fig.suptitle('Mach number evolution when heating, $\gamma = %.1f$'%gam, fontsize=12, y=0.93)

for M0 in [.1, .2, .3, .4, .5, .6, .7, .8, .9, 1.2, 1.5, 2., 3., 4., 6., 10., 20.]:
#for M0 in [.1, .5]:
    Timax = ray.maxTiratio_Mach(M0, gam)
    Ti    = np.log10(np.logspace(1., Timax, npoints+1))
    Mach  = ray.SubMach_TiTicri(Ti/Timax, gam) if M0 < 1 else ray.SupMach_TiTicri(Ti/Timax, gam)
    plt.plot(Ti, Mach, '-', color='#ff0000')

plt.axis([1., 10, .1, 10.])
plt.ylabel('Mach number', fontsize=10)
plt.xlabel('$T_i/T_{i0}$', fontsize=10)
plt.loglog()
#plt.ticklabel_format(axis='both', style='plain')
#plt.minorticks_on()
plt.grid(which='major', linestyle=':', alpha=0.5)
plt.grid(which='minor', linestyle=':', alpha=0.5)
fig.savefig('Rayleigh-heating.pdf', bbox_inches='tight')
plt.show()

max_stars_repo_name: PierreMignerot/aerokit

max_issues_repo_name: PierreMignerot/aerokit

max_forks_repo_name: PierreMignerot/aerokit


text:
#!/usr/bin/env python

import rospy
import cv2
import torch
import numpy as np
from network import DensenetUnetHybrid
import image_utils

from std_msgs.msg import Header
from sensor_msgs.msg import Image, CameraInfo
from cv_bridge import CvBridge, CvBridgeError

class Predictor:
    def __init__(self, republish_color=False):
        self.image_callback_called = False
        self.republish_color = republish_color ### only so that image and depth are synchronized. it saves using an approximate filter subscriber on other nodes
        self.pI = PredictInterface()
        rospy.init_node('predictor', anonymous=True)
        rospy.Subscriber("image_in", Image, self.image_callback) ## to image
        if self.republish_color:
            self.image_pub = rospy.Publisher('image/image_raw', Image, queue_size=1)
            self.camera_info_image_pub = rospy.Publisher('image/camera_info', CameraInfo, queue_size=1)

        self.depth_pub = rospy.Publisher('depth/image_raw', Image, queue_size=1)
        self.camera_info_depth_pub = rospy.Publisher('depth/camera_info', CameraInfo, queue_size=1)
        self.bridge = CvBridge()
        rospy.loginfo('Predictor instantiated.')
        rospy.spin()

    def image_callback(self, msg):
        try:
            the_now_time = rospy.Time.now()
            if not self.image_callback_called:
                rospy.loginfo('Received at least one image.')
                self.image_callback_called = True
            cv_image = self.bridge.imgmsg_to_cv2(msg, "bgr8")
            outimage = self.pI.predict(cv_image)
            #i need to stamp or they will not synchronize!
            outDepthMsg = self.bridge.cv2_to_imgmsg(outimage, "passthrough")
            outDepthMsg.header.stamp = the_now_time
            outDepthMsg.header.frame_id = "camera_rgb_optical_frame" ## i kinda need to read this from the original camerrainfo

            self.depth_pub.publish(outDepthMsg)
            stdDepthCameraInfoMsg = CameraInfo()
            stdDepthCameraInfoMsg.width = 304
            stdDepthCameraInfoMsg.height = 228
            stdDepthCameraInfoMsg.distortion_model = "plumb_bob"
            stdDepthCameraInfoMsg.D = [0.0, 0.0, 0.0, 0.0, 0.0]
            stdDepthCameraInfoMsg.K = [575.8157348632812, 0.0, 314.5, 0.0, 575.8157348632812, 235.5, 0.0, 0.0, 1.0]
            stdDepthCameraInfoMsg.R = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
            stdDepthCameraInfoMsg.P = [575.8157348632812, 0.0, 314.5, 0.0, 0.0, 575.8157348632812, 235.5, 0.0, 0.0, 0.0, 1.0, 0.0]
            stdDepthCameraInfoMsg.header = Header()
            stdDepthCameraInfoMsg.header.stamp = the_now_time
            stdDepthCameraInfoMsg.header.frame_id = "camera_rgb_optical_frame" ## i kinda need to read this from the original camerrainfo
            self.camera_info_depth_pub.publish(stdDepthCameraInfoMsg)
            if self.republish_color:
                outImageMsg = self.bridge.cv2_to_imgmsg(cv_image, "bgr8")
                outImageMsg.header.stamp = the_now_time
                outImageMsg.header.frame_id = "camera_rgb_optical_frame"
                self.image_pub.publish(outImageMsg)
                stdImageCameraInfoMsg = CameraInfo()
                stdImageCameraInfoMsg.width = 640
                stdImageCameraInfoMsg.height = 480
                stdImageCameraInfoMsg.distortion_model = "plumb_bob"
                stdImageCameraInfoMsg.D = [0.0, 0.0, 0.0, 0.0, 0.0]
                stdImageCameraInfoMsg.K = [575.8157348632812, 0.0, 314.5, 0.0, 575.8157348632812, 235.5, 0.0, 0.0, 1.0]
                stdImageCameraInfoMsg.R = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
                stdImageCameraInfoMsg.P = [575.8157348632812, 0.0, 314.5, 0.0, 0.0, 575.8157348632812, 235.5, 0.0, 0.0, 0.0, 1.0, 0.0]
                stdImageCameraInfoMsg.header = Header()
                stdImageCameraInfoMsg.header.stamp = the_now_time
                stdImageCameraInfoMsg.header.frame_id = "camera_rgb_optical_frame" ## i kinda need to read this from the original camerrainfo
                self.camera_info_image_pub.publish(stdImageCameraInfoMsg)

        except CvBridgeError as e:
          print(e)

class PredictInterface:
    def __init__(self):
        # switch to CUDA device
        self.device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')
        rospy.loginfo('Use GPU: {}'.format(str(self.device) != 'cpu'))

        # load model
        rospy.loginfo('Loading model...')
        self.model = DensenetUnetHybrid.load_pretrained(device=self.device)
        self.model.eval()

        self.predict_called = False

    def predict(self, img):
        img = image_utils.scale_image(img)
        img = image_utils.center_crop(img)
        inp = image_utils.img_transform(img)
        inp = inp[None, :, :, :].to(self.device)

        # inference
        if not self.predict_called:
            rospy.loginfo('Running the image through the network...')
            self.predict_called = True
        output = self.model(inp)

        # transform and output the results
        output = output.cpu()[0].data.numpy()
        pred = np.transpose(output, (1, 2, 0))
        return pred[:, :, 0]

def predict_img(img_path):
    """Inference a single image."""

    # load image
    img = cv2.imread(img_path)[..., ::-1]

    pI = PredictInterface()

    output = pI.predict(img)
    image_utils.show_img_and_pred(img, output)

if __name__ == '__main__':
    try:
        p = Predictor(republish_color=True)

    except rospy.ROSInterruptException:
        pass


text:
"""
*-------------------------------------------------------------------------------------
*                           #### MACHINERY PROBLEM ####
**------------------------------------------------------------------------------------
*                                                           Juan Javaloyes     2018
*                                                           javaloyes.juan@gmail.com

* Source: Grado Ingeniería Química. Universidad de Alicante.
          Asignatura: Simulaicón, optimización y diseño de procesos químicos

"""

"""
 PROBLEM STATMENT:
    A company manufacture four types of machinery. The factory is divided in three sections.
	The first section has available 960 h/week, the second 1110 h/week and the third 400 h/week.
	Each machinery unit requires the time given in "time_x_section" at each section.
	Determine the number of units of machinery for each type that should be manufacture per week
	to maximize profit.
"""

from pyomo.environ import *
m = ConcreteModel()

M = m.M = Set(initialize =  ['m1', 'm2', 'm3', 'm4'], ordered = True)
S = m.S = Set(initialize =  ['s1', 's2', 's3'], ordered = True)

profit   = {'m1':  12, 'm2':    8, 'm3':  12, 'm4':17}
max_time = {'s1': 960, 's2': 1110, 's3': 400}

time_x_section = { ('m1','s1'): 6 , ('m1','s2'): 3 , ('m1','s3'): 2,
                   ('m2','s1'): 4 , ('m2','s2'): 3 , ('m2','s3'): 1,
				   ('m3','s1'): 4 , ('m3','s2'): 6 , ('m3','s3'): 2,
				   ('m4','s1'): 8 , ('m4','s2'): 9 , ('m4','s3'): 1}

x = m.x = Var( M, within = PositiveIntegers )

m.value = Objective( expr = sum( profit[i] * m.x[i] for i in M),
                     sense = maximize )

def constraint_rule(m, j):
	return sum( time_x_section[i,j] * x[i] for i in M) <= max_time[j]
m.constraint = Constraint(S, rule = constraint_rule)

SolverFactory('glpk').solve(m, tee = True)

m.pprint()

max_stars_repo_path: # 02 -Case Studies/# 00 - Machinery problem/machinery_problem.py

max_stars_repo_name: CAChemE/pyomo-tutorial

max_issues_repo_name: CAChemE/pyomo-tutorial

max_forks_repo_name: CAChemE/pyomo-tutorial

reason: from pyomo


text:
# author: Ovidiu Mura
# date: May 22, 2019

# Reference: http://chris35wills.github.io/simple_fft_filter/
def filter_amps_below(signal, amp_below):
    # signal: signal to be filtered
    # amp: the threshold amplitude; the amps below this amp will be considered
    a = []
    for i in range(1,238):
        if(i<237):
            b=0
            if(len(a)!=0):
                b = a[-1]
            if(signal[i]<b+amp_below):
                a.insert(i, signal[i])
    a[0] = 0
    return a

# https://plot.ly/python/fft-filters/
from scipy import signal
import matplotlib.pyplot as plt
import numpy as np

import wave
from scipy.io import wavfile
from scipy import signal
from scipy.fftpack import fft,fftfreq, ifft, rfft, rfftfreq, irfft

# It implements the Goertzel, the Low pass, High pass, and Band pass filters.
class Filter:

    def __init__(self):
        self.signal = None
        self.goertzel_freqs = {}
        self.norm = None
        self.coeffs = None

    # https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lfilter.html
    # https://tomroelandts.com/articles/how-to-create-a-simple-low-pass-filter
    def lowpass(self):
        samplerate, x0 = wavfile.read("SineWave_440Hz.wav")
        #xn = np.array(x0)
        xn = fft(x0)
        t = np.arange(len(xn))/float(samplerate)
        cut_off = 500.0
        Wn = (float(cut_off) / (float(samplerate) / 2.0))

        b, a = signal.butter(3, Wn, 'low')

        y = signal.filtfilt(b, a, xn)

        yy = ifft(y)
        xx = fftfreq(len(xn), 1/samplerate)

        plt.title("Sine Wave - Low pass filter")
        plt.plot(xx[:10000].real, x0[:10000].real, 'red')
        plt.plot(xx[:10000].real, yy[:10000].real, 'green')
        plt.legend(('sine wave', 'low pass filtered sine wave'), loc='best')
        plt.xlabel("no of samples")
        plt.ylabel("amplitudes")
        plt.grid(True)
        plt.show()

    # https://www.programcreek.com/python/example/59508/scipy.signal.butter
    def highpass(self):
        cut_off = 500.0
        samplerate, xn = wavfile.read("SineWave_440Hz.wav")

        f1 = fft(xn)

        Wn = (float(cut_off) / (float(samplerate) / 2.0))

        b, a = signal.butter(3, Wn, 'high')

        z = signal.filtfilt(b, a, f1)

        f2 = ifft(z)
        xx = fftfreq(len(f2), 1/samplerate)
        plt.title("Sine Wave - High pass filter")
        plt.plot(xx[:10000].real, xn[:10000].real, 'r')
        plt.plot(xx[:10000].real, f2[:10000].real, 'green')
        plt.legend(('sine wave', 'high pass filtered sine wave'), loc='best')
        plt.xlabel("no of samples")
        plt.ylabel("amplitudes")
        plt.grid(True)
        plt.show()

    # https://github.com/neurotechuoft/Wall-EEG/blob/master/Code/OpenBCIPy/src/butter_tingz.py
    def bandpass(self):
        samplerate, xn = wavfile.read("SineWave_440Hz.wav")
        nyq = 0.5 * 250
        lowcut = 10
        highcut = 30
        low = lowcut / nyq
        high = highcut / nyq

        x = fft(xn)

        b, a = signal.butter(3, [low, high], btype='band')

        filtered_data = signal.lfilter(b, a, x)

        y = ifft(filtered_data)
        xx = fftfreq(len(xn), 1/samplerate)
        plt.title("Sine Wave - Band pass filter")
        plt.plot(xx[:30000].real, xn[:30000].real, 'green')
        plt.plot(xx[:30000].real, y[:30000].real, 'r')
        plt.legend(('sine wave', 'band pass filtered sine wave'), loc='best')
        plt.xlabel("no of samples")
        plt.ylabel("amplitudes")
        plt.grid(True)
        plt.show()

    # https://en.wikipedia.org/wiki/Goertzel_algorithm , Equation (6)
    def goertzel(self, x, f):
        """x is an array of samples, f is the target frequency.
        Returns the output magnitude."""
        import math
        w0 = float((2*math.pi*f)/44100)
        n = len(x)
        self.norm = np.exp(1j * w0 * n)
        self.coeffs = np.exp(np.array([-1j * w0 * k for k in range(n)]))
        y = np.abs(self.norm * np.dot(self.coeffs, x))
        self.goertzel_freqs[f] = y
        return y

# It reads the sign frequency from a given file using Goertzel filter.
def get_sine_freq(file_name='SineWave_440Hz.wav', p=True):
    f = Filter()

    samplerate, data = wavfile.read(file_name)

    freqs = {}

    for i in range(150):
        m = f.goertzel(data, 420+i)
        freqs[420+i] = m

    m = max(freqs.values())
    idx = list(freqs.values()).index(m)
    freq = list(freqs.keys())[idx]

    if(p):
        plt.title("Sine Wave - Goertzel filter")
        plt.plot(freqs.keys(), freqs.values(), color='green')
        plt.ylabel("amplitudes")
        plt.xlabel("frequency (Hz)")
        plt.show()
    return freq

max_stars_repo_name: ovidiu-mura/emergency-sound-detection

max_issues_repo_name: ovidiu-mura/emergency-sound-detection

max_forks_repo_name: ovidiu-mura/emergency-sound-detection

reason: import numpy,from scipy


text:
import numpy as np
import os
import csv

import pickle
import sys

import matplotlib
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
from matplotlib import cm

from mpl_toolkits.mplot3d import Axes3D

import scipy.interpolate

from pyevtk.hl import pointsToVTK

import pyvtk

DATA_FORMAT_VERSION = 2

var_params = []

#path = "output/restart_002"
# record_times = [0, 50, 100]  # in Myrs
# record_depths = [31e3, 95e3, 45e3]

PLOT = 21
PLOT_TIME = 0
PLOT_DEPTH_BOTTOM = 5
PLOT_DEPTH = 0
PLOT_TIME1 = 0
PLOT_TIME2 = 1
PLOT_DEPTH_BELOW_THRUST = 0
ADD_INITIAL_DEPTH = True
# 1 = non-dim near-bottom temperature, at time PLOT_TIME
# 2 = non-dim change in temperature at PLOT_DEPTH_BELOW_THRUST, from PLOT_TIME1 to PLOT_TIME2

#path = "output/steadystate_003"
#record_times = [1800, 1900]
#record_depths = [35e3, 75e3, 125e3]

path = "output/001-ini"
record_times = np.array([10,180,200])
record_depths = np.array([10e3,99e3,20e3])

depth_adjust = record_depths * 0.0

ndatapoints = 0

rec_values = []
rec_times = []
rec_depths = []
rec_params = []
rec_ranks = []

for d in os.walk(path):
    meta_read = False
    params_read = False

    start_time = None
    irec = 0

    if ("params" in d[2]):
        params = []
        rparamfile = open(d[0] + "/" + "params")
        for row in rparamfile:
            params.append(eval(row))
        rparamfile.close()
        if PLOT == 20 or PLOT == 21:
            if ADD_INITIAL_DEPTH:
                depth_adjust = record_depths * 0.0
                depth_adjust = depth_adjust + params[0]
            else:
                depth_adjust = record_depths * 0.0
        else:
            depth_adjust = record_depths * 0.0
    else:
        # skipped
        continue

    if ("meta" in d[2]):
        rmetafile = open(d[0] + "/" + "meta")
        csvmetar = csv.reader(rmetafile, delimiter=",", quotechar='"')
        irow = 0
        alltimesread = False
        for row in csvmetar:
            irow = irow + 1
            if irow == 1:
                if int(row[0]) != DATA_FORMAT_VERSION:
                    raise Exception("Incompatible data format")
                else:
                    continue
            elif irow == 2:
                mpi_rank = int(row[0])
                meta_read = True
                continue
            elif irow == 3:
                # just the header
                pass

            if irow == 4:
                # we've got data!
                start_time = float(row[1])
                rec_values.append([])
                rec_times.append([])
                rec_depths.append([])
                rec_params.append([])

                for p in params:
                    rec_params[-1].append(p)
                rec_ranks.append(mpi_rank)
                ndatapoints = ndatapoints + 1

            if irow > 3:
                if (float(row[1])-start_time) >= record_times[irec]:
                    irec = irec + 1
                    rec_values[-1].append([])
                    rec_depths[-1].append([])
                    time = float(row[1]) - start_time
                    rec_times[-1].append(time)
                    rnodefile = open(d[0] + "/" + "nodes." + str(row[0]), "rb")
                    csvnoder = csv.reader(rnodefile, delimiter=",", quotechar='"')
                    inoderow = 0
                    idepthrec = 0
                    skip_header = True
                    DEPTHCOL = 2
                    TEMPCOL = 3
                    alldepthsread = False
                    lastdepth = -np.inf
                    for noderow in csvnoder:
                        if skip_header:
                            skip_header = False
                            continue
                        if inoderow == 0:
                            for i in range(len(record_depths)):
                                if float(noderow[DEPTHCOL]) > record_depths[i] + depth_adjust[i]:
                                    # first depth recorder is larger than a requested recording depth
                                    idepthrec = i
                                    rec_values[-1][-1].append(np.nan)
                                    rec_depths[-1][-1].append(np.nan)
                                    lastdepth = float(noderow[DEPTHCOL])
                                    inoderow = inoderow + 1
                                    continue
                        #if idepthrec >= len(record_depths):
                        #    # we have already recorded more than requested amount of depths
                        #    rec_depths[-1][-1].append(-999)
                        #    rec_values[-1][-1].append(-999)
                        #    idepthrec = idepthrec + 1
                        #    lastdepth = float(noderow[DEPTHCOL])
                        elif float(noderow[DEPTHCOL]) > record_depths[idepthrec] + depth_adjust[idepthrec]:
                            idepthrec = idepthrec + 1
                            rec_values[-1][-1].append(float(noderow[TEMPCOL]))
                            rec_depths[-1][-1].append(float(noderow[DEPTHCOL]))
                            lastdepth = float(noderow[DEPTHCOL])
                        if idepthrec >= len(record_depths):
                            alldepthsread = True
                            break
                        inoderow = inoderow + 1
                    rnodefile.close()
                    while len(rec_depths[-1][-1]) < len(record_depths):
                        # fill unfound depths with NaNs
                        rec_depths[-1][-1].append(np.nan)
                        rec_values[-1][-1].append(np.nan)

                    if irec >= len(record_times):
                        alltimesread = True
                        break
        # all times checked
        # fill unfound times with NaNs
        while len(rec_values[-1]) < len(record_times):
            rec_values[-1].append([])
            rec_depths[-1].append([])
            rec_times[-1].append(np.nan)
            while len(rec_depths[-1][-1]) < len(record_depths):
                rec_depths[-1][-1].append(np.nan)
                rec_values[-1][-1].append(np.nan)

        rmetafile.close()

######
## SECONDARY HANDLING
######

print ndatapoints

np_rec_values = np.array(rec_values)
print "---"
print np_rec_values.shape

np_rec_ranks = np.array(rec_ranks)
np_rec_times = np.array(rec_times)
np_rec_depths = np.array(rec_depths)

np.set_printoptions(threshold=np.nan)

#Tscale = np_rec_values[:,0,1]
#dT = (np_rec_values[:,1,1] - np_rec_values[:,0,1])/Tscale

if PLOT == 1:
    plottitle = "bottom temp change from time zero to t=" + str(record_times[PLOT_TIME])
    Tscale = np_rec_values[:,0,PLOT_DEPTH_BOTTOM]
    dT = np_rec_values[:,PLOT_TIME,PLOT_DEPTH_BOTTOM]/Tscale
    p1 = np.array(rec_params)[:,0]
    p2 = np.array(rec_params)[:,1]
elif PLOT == 2:
    plottitle = "temp change from t=" + str(record_times[PLOT_TIME1]) + " to t=" + str(record_times[PLOT_TIME2]) + " at z = " + str(record_depths[PLOT_DEPTH_BELOW_THRUST])
    Tscale = np_rec_values[:,0,PLOT_DEPTH_BOTTOM]
    dT = np_rec_values[:,PLOT_TIME2,PLOT_DEPTH_BELOW_THRUST] - np_rec_values[:,PLOT_TIME1,PLOT_DEPTH_BELOW_THRUST]
    dT = dT / Tscale
    p1 = np.array(rec_params)[:,0]
    p2 = np.array(rec_params)[:,1]
elif PLOT == 20:
    plottitle = "what are we plotting here? " + str(record_times[PLOT_TIME])
    Tscale = 1350.0
    print np_rec_values.shape
    print np_rec_times.shape
    pickedT = np_rec_values[:,PLOT_TIME,PLOT_DEPTH]
    ip0 = []
    ip1 = []
    ip2 = []
    iy0 = []
    for i in range(0,len(rec_params)):
        if rec_params[i][2] > -1:
            ip0.append(rec_params[i][0])   # thickness of the overthrust sheet
            ip1.append(rec_params[i][1][1]-rec_params[i][1][0]) # original extent of the lithosphere
            ip2.append(rec_params[i][2]) # erosion speed
            iy0.append(pickedT[i]) # value itself
    ip0 = np.array(ip0)
    ip1 = np.array(ip1)
    ip2 = np.array(ip2)
    iy0 = np.array(iy0)

    p0 = ip0
    p0lab = "overthrust sheet thickness"
    p1 = ip1
    p1lab = "orig. lithosphere thickness"
    p2 = ip2
    p2lab = "erosion speed"
    y0 = iy0
    y0lab = "temperature"

    print ip0
    print ip1
    print ip2
    print iy0

elif PLOT == 21:
    plottitle = "what are we plotting here? " + str(record_times[PLOT_TIME])
    Tscale = 1350.0
    print np_rec_values.shape
    print np_rec_times.shape
    pickedT = np_rec_values[:,PLOT_TIME,PLOT_DEPTH]
    ip0 = []
    ip1 = []
    ip2 = []
    iy0 = []
    for i in range(0,len(rec_params)):
        if rec_params[i][2] > -1:
            ip0.append(rec_params[i][0])   # thickness of the overthrust sheet
            ip1.append(rec_params[i][1]) # bottom temp
            ip2.append(rec_params[i][2]) # erosion speed
            iy0.append(pickedT[i]) # value itself
    ip0 = np.array(ip0)
    ip1 = np.array(ip1)
    ip2 = np.array(ip2)
    iy0 = np.array(iy0)

    p0 = ip0
    p0lab = "overthrust sheet thickness"
    p1 = ip1
    p1lab = "lithosphere bottom temp"
    p2 = ip2
    p2lab = "erosion speed"
    y0 = iy0
    y0lab = "temperature"

    print ip0
    print ip1
    print ip2
    print iy0

x0i = np.linspace(min(p0), max(p0))
x1i = np.linspace(min(p1), max(p1))
x2i = np.linspace(min(p2), max(p2))
x0, x1, x2 = np.meshgrid(x0i, x1i, x2i)

p1 = np.ma.array(p1, mask=np.isnan(y0))
p2 = np.ma.array(p2, mask=np.isnan(y0))
p0 = np.ma.array(p0, mask=np.isnan(y0))
y0 = np.ma.array(y0, mask=np.isnan(y0))
print "----"
print len(p0)
print len(p1)
print len(p2)
print len(y0)
print np.ma.vstack((p0,p1,p2)).T.shape
orig_x = np.ma.vstack((p0,p1,p2)).T
print "---"
##print len(x0), len(x0[0]), len(x0[0][0])
queryx = np.vstack((x0.flatten(), x1.flatten(), x2.flatten())).T
##print queryx.shape
ipolate = scipy.interpolate.LinearNDInterpolator(points=orig_x, values=y0)
zi = ipolate(queryx)
##print len(queryx)
##print len(zi)
#fig = plt.figure()
#ax = fig.add_subplot(111, projection='3d')
#ax.scatter(queryx[:,0], queryx[:,1], queryx[:,2])
#plt.show()
#zi = scipy.interpolate.griddata((p1,p2), y0, (xi[None,:],yi[:,None]), method='cubic')
NN = len(queryx[:,0])
picks = np.arange(0,NN,1)
pointsToVTK("./points", queryx[picks,0], queryx[picks,1], queryx[picks,2], data = {"temp" : zi[picks]})
print queryx.shape
print zi.shape
#sys.exit()

do3d = False

for x in np.unique(queryx[picks,2]):
    idx = np.where(queryx[picks,2] == x)
    fig = plt.figure()
    if do3d:
        ax = fig.add_subplot(111,projection='3d')
        ax.scatter((queryx[picks,0])[idx], (queryx[picks,1])[idx], (queryx[picks,2])[idx], c=(zi[picks])[idx], cmap=plt.hot(), s=5, edgecolors='none')
    else:
        ax = fig.add_subplot(111)
        sc = ax.scatter((queryx[picks,0])[idx], (queryx[picks,1])[idx], c=(zi[picks])[idx], cmap=plt.hot(), edgecolors='none')
        plt.colorbar(sc)
        #cs = ax.contour((queryx[picks,0])[idx], (queryx[picks,1])[idx], (zi[picks])[idx])
        #ax.clabel(cs)
        ax.set_xlabel(p0lab)
        ax.set_ylabel(p1lab)
    plt.show()

#contour(xi, yi, zi)

#plt.figure()
#CS = plt.contour(X, Y, Z)
#plt.show()

max_stars_repo_name: larskaislaniemi/PyLitho

max_issues_repo_name: larskaislaniemi/PyLitho

max_forks_repo_name: larskaislaniemi/PyLitho

reason: import numpy,import scipy


text:
import os
import errno
import numpy as np
import itertools
import six

from scipy import sparse

import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.patches import Wedge
from sklearn.metrics import confusion_matrix

from textwrap import wrap

def savefig_and_close(fig, figname, path='', bbox_extra_artists=None):
    filename = os.path.join(path, figname)
    fig.savefig(filename, bbox_extra_artists=bbox_extra_artists,
                bbox_inches='tight')
    fig.clear()
    plt.close(fig)

def newfig(name, **kwargs):
    fig = plt.figure(name, **kwargs)
    fig.clf()
    return fig

def savefig(fig, path='figures', prefix='weak_labels_', extension='svg'):
    fig.tight_layout()
    name = fig.get_label()
    filename = "{}{}.{}".format(prefix, name, extension)
    try:
        os.makedirs(path)
    except OSError as exception:
        if exception.errno != errno.EEXIST:
            raise
    fig.savefig(os.path.join(path, filename))

def plot_data(x, y, loc='best', save=True, title='data', cmap='Paired'):
    if sparse.issparse(x):
        x = x.toarray()

    fig = newfig('data', figsize=(5, 3))
    ax = fig.add_subplot(111)

    classes = np.unique(y)
    n_c = float(len(classes))

    cmap = cm.get_cmap(cmap)
    for i, y_i in enumerate(classes):

        ax.scatter(x[(y == y_i).flatten(), 0], x[(y == y_i).flatten(), 1],
                   c=cmap(i/n_c), s=30, edgecolors=None, alpha=.8, lw=0.1,
                   label='{}'.format(y_i))
    ax.set_xlabel('$x_0$')
    ax.set_ylabel('$x_1$')
    ax.set_title(title)
    ax.axis('equal')
    ax.legend(loc=loc)
    ax.grid(True)
    if save:
        savefig(fig)
    return fig

class MyFloat(float):
    def _remove_leading_zero(self, value, string):
        if 1 > value > -1:
            string = string.replace('0', '', 1)
        return string

    def __str__(self):
        string = super(MyFloat, self).__str__()
        return self._remove_leading_zero(self, string)

    def __format__(self, format_string):
        string = super(MyFloat, self).__format__(format_string)
        return self._remove_leading_zero(self, string)

# TODO use this or other heatmap to visualize confusion matrix
def plot_df_heatmap(df, normalize=None, title='Heat-map',
                    cmap=plt.cm.Blues, colorbar=False, ylabel=None,
                    xlabel=None):
    """
    This function prints and plots the confusion matrix.
    Normalization can be applied by setting `normalize=True`.

    normalize : 'rows', 'cols' (default=None)
    """
    rows = df.index.values
    columns = df.columns.values
    M = df.values

    if normalize == 'rows':
        M = M.astype('float') / M.sum(axis=1)[:, np.newaxis]
    if normalize == 'cols':
        M = M.astype('float') / M.sum(axis=0)[np.newaxis, :]

    xlabel = df.columns.name
    ylabel = df.index.name

    return plot_heatmap(M, columns=columns, rows=rows, cmap=cmap,
                        colorbar=colorbar, title=title, ylabel=ylabel,
                        xlabel=xlabel)

def plot_confusion_matrix(M, columns=None, rows=None, cmap=plt.cm.Blues,
                          colorbar=False, fig=None, title='Heat-map',
                          ylabel='True label', xlabel='Predicted label',
                          **kwargs):
    return plot_heatmap(M=M, columns=columns, rows=rows, cmap=cmap,
                        colorbar=colorbar, fig=fig, title=title, ylabel=ylabel,
                        xlabel=xlabel, **kwargs)

def plot_heatmap(M, columns=None, rows=None, cmap=plt.cm.Blues, colorbar=False,
                 fig=None, ax=None, title='Heat-map', ylabel=None, xlabel=None):
    if columns is None:
        columns = [str(i) for i in range(M.shape[1])]
    if rows is None:
        rows = [str(i) for i in range(M.shape[0])]

    h_size = 5 # len(columns)*.5 + 2
    v_size = 4 # len(rows)*.5 + 2

    if fig is None:
        fig = plt.figure(figsize=(h_size, v_size))
    if ax is None:
        ax = fig.add_subplot(111)
    im = ax.imshow(M, interpolation='nearest', cmap=cmap)
    if colorbar:
        fig.colorbar(im)
    if h_size < 4:
        title = "\n".join(wrap(title, 30))
    ax.set_title(title)
    column_tick_marks = np.arange(len(columns))
    ax.set_xticks(column_tick_marks)
    ax.set_xticklabels(columns, rotation=45, ha='right')
    row_tick_marks = np.arange(len(rows))
    ax.set_yticks(row_tick_marks)
    ax.set_yticklabels(rows)
    if xlabel is not None:
        ax.set_xlabel(xlabel)
    if ylabel is not None:
        ax.set_ylabel(ylabel)

    thresh = np.nanmin(M) + ((np.nanmax(M)-np.nanmin(M)) / 2.)
    are_ints = M.dtype in ['int', 'int32', 'int64']
    for i, j in itertools.product(range(M.shape[0]), range(M.shape[1])):
        # fontsize is adjusted for different number of digits
        if are_ints:
            num_text = str(M[i, j])
        else:
            num_text = '{:0.2f}'.format(MyFloat(M[i, j]))

        if np.isfinite(M[i, j]):
            ax.text(j, i, num_text, horizontalalignment="center",
                    verticalalignment="center", color="white" if M[i, j] >
                    thresh else "black", fontsize=16-len(num_text))

    fig.tight_layout()
    return fig

def dual_half_circle(center, radius, angle=0, ax=None, colors=('w','k'),
                     **kwargs):
    """
    Add two half circles to the axes *ax* (or the current axes) with the
    specified facecolors *colors* rotated at *angle* (in degrees).
    """
    if ax is None:
        ax = plt.gca()
    theta1, theta2 = angle, angle + 180
    w1 = Wedge(center, radius, theta1, theta2, fc=colors[0], **kwargs)
    w2 = Wedge(center, radius, theta2, theta1, fc=colors[1], **kwargs)
    for wedge in [w1, w2]:
        #ax.add_artist(wedge)
        ax.add_patch(wedge)
    return [w1, w2]

def test_dual_half_circle_main():
    fig, ax = plt.subplots()
    dual_half_circle((0.5, 0.5), radius=0.3, angle=90, ax=ax)
    ax.axis('equal')
    plt.show()

def plot_multilabel_scatter(X, Y, cmap=cm.get_cmap('tab20'), edgecolor='k',
                            linewidth=0.4, title=None, fig=None, ax=None,
                            radius_scaler=20.0, **kwargs):
    X_std = X.std(axis=0)
    if X.shape[1] > 2:
        biggest_variance = np.argsort(X_std)[-2:]
        X_std = X_std[biggest_variance]
        X = X[:,biggest_variance]

    X_min = X.min(axis=0)
    X_max = X.max(axis=0)
    n_classes = Y.shape[1]

    radius = ((X_max - X_min)/radius_scaler)[:2].min()
    #radius = (X.max() - X.min())/radius_scaler

    if fig is None:
        fig = plt.figure(figsize=(4, 3))
    if ax is None:
        ax = fig.add_subplot(111)

    for x, y in zip(X, Y):
        theta2s = np.cumsum(np.true_divide(y, y.sum())*360.0)
        theta1 = 0
        if np.isfinite(theta2s[0]):
            for i, theta2 in enumerate(theta2s):
                if theta1 != theta2:
                    w = Wedge(x[:2], radius, theta1, theta2, ec=edgecolor, lw=linewidth,
                              fc=cmap(i), **kwargs)
                    ax.add_patch(w)
                    theta1 = theta2
        else:
            # Not belong to any class
            print('Do not belong to any class')
            w = Wedge(x[:2], radius, 0, 360, ec='black', lw=linewidth,
                      fc='white', **kwargs)
            ax.add_patch(w)

    ax.set_xlim([X_min[0]-X_std[0], X_max[0]+X_std[0]])
    ax.set_ylim([X_min[1]-X_std[1], X_max[1]+X_std[1]])
    ax.axis('equal')
    if title is not None:
        ax.set_title(title)
    return fig

def test_multilabel_plot():
    X = np.array([[0,0], [0,1], [1,0], [1,1]])
    Y = np.array([[0, 0, 1],
                  [0, 1, 1],
                  [1, 0, 0],
                  [1, 1, 1]])
    plot_multilabel_scatter(X, Y)
    plt.show()

def plot_errorbar(data, x=None, fmt='--o', title='errorbar', elinewidth=1.0,
                  perrorevery=0.2, legend=None, fig=None, **kwargs):
    """

    paramters
        data: np.array or list of np.array
            If it is a list, each np.array is considered as an errorbar line

        errorevery: float
            Percentage of errorbars with respect to the number of samples
    """
    return_fig = fig is None

    if fig is None:
        fig = newfig(title, figsize=(5, 3))

    ax = fig.add_subplot(111)
    ax.set_title(title)

    if type(data) is np.ndarray:
        data = (data,)

    for i, matrix in enumerate(data):
        errorevery = int(1 / perrorevery)
        if errorevery < 1:
            errorevery = 1

        if x is None:
            x = range(matrix.shape[1])

        means = matrix.mean(axis=0)
        stds = matrix.std(axis=0)

        ax.errorbar(x=x, y=means, yerr=stds, elinewidth=(len(data)-i)*elinewidth,
                    errorevery=errorevery, capthick=(len(data)-i)*elinewidth,
                    capsize=4, **kwargs)
    if legend is not None:
        ax.legend(legend)

    if return_fig:
        return fig
    return ax

def render_mpl_table(data, col_width=3.0, row_height=0.625, font_size=14,
                     header_color='#40466e', row_colors=['#f1f1f2', 'w'],
                     edge_color='w', bbox=[0, 0, 1, 1], header_columns=0,
                     ax=None, **kwargs):
    """
    source: https://stackoverflow.com/questions/19726663/how-to-save-the-pandas-dataframe-series-data-as-a-figure
    """
    if ax is None:
        size = (np.array(data.shape[::-1]) + np.array([0, 1])) * np.array([col_width, row_height])
        fig, ax = plt.subplots(figsize=size)
        ax.axis('off')

    mpl_table = ax.table(cellText=data.values, bbox=bbox, colLabels=data.columns, **kwargs)

    mpl_table.auto_set_font_size(False)
    mpl_table.set_fontsize(font_size)

    for k, cell in six.iteritems(mpl_table._cells):
        cell.set_edgecolor(edge_color)
        if k[0] == 0 or k[1] < header_columns:
            cell.set_text_props(weight='bold', color='w')
            cell.set_facecolor(header_color)
        else:
            cell.set_facecolor(row_colors[k[0] % len(row_colors)])
    return ax

def export_df(df, filename, path='', extension='pdf',
              stylesheets=['style.css']):
    """
    df : pandas.DataFrame
    """
    filename = os.path.join(path, filename)
    ax = render_mpl_table(df)
    fig = ax.figure
    fig.tight_layout()
    fig.savefig("{}.{}".format(filename, extension))

def plot_history(history, model=None, X_test=None, y_test=None, best_measure='val_loss', fig=None):
    selected_epoch = np.argmin(history.history[best_measure])

    def show_best(ax, history, field, selected_epoch=None):
        if 'acc' in field:
            best_epoch = np.argmax(history.history[field])
        else:
            best_epoch = np.argmin(history.history[field])
        _ = ax.scatter(best_epoch, history.history[field][best_epoch],
                       color='white', edgecolor='red', marker='*',
                       s=200, zorder=3, label='Best val.')
        if selected_epoch is not None:
            _ = ax.scatter(selected_epoch, history.history[field][selected_epoch],
                           color='gold', edgecolor='black', marker='*',
                           s=150, zorder=3, label='Selected')

    metrics = [key for key in history.history.keys() if not key.startswith('val_')]

    n_fig = len(metrics) + (model is not None)

    if fig is None:
        fig = plt.figure(figsize=(n_fig*3, 3))

    for i, metric in enumerate(metrics):
        ax = fig.add_subplot(1, n_fig, i + 1)
        ax.set_title(metric)
        _ = ax.plot(history.history[metric], label='Train.')
        metric = 'val_' + metric
        if metric in history.history.keys():
            _ = ax.plot(history.history[metric], label='Val.')
            show_best(ax, history, metric, selected_epoch=selected_epoch)
        ax.legend()

    if model is not None:
        # FIXME predict destroys the history
        clf_proba_wt_test = model.predict(X_test)
        # FIXME I am currently restoring it from history
        model.history = history
        clf_pred_wt_test = np.argmax(clf_proba_wt_test, axis=1)
        cm = confusion_matrix(y_test, clf_pred_wt_test)

        ax = fig.add_subplot(1, n_fig, n_fig)
        acc = (y_test == clf_pred_wt_test).mean()
        _ = plot_confusion_matrix(cm, ax=ax, title='Test acc. {:.3}'.format(acc))
    return fig

max_stars_repo_name: perellonieto/RecyclingLabels

max_issues_repo_name: perellonieto/RecyclingLabels

max_forks_repo_name: perellonieto/RecyclingLabels


text:
'''
Copyright 2020 CITA

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

20201019_desnity.py

Shows Historgam of CT scan density

Author: Sebastian Gatz

'''

import numpy as np
from vedo import *
from vedo.pyplot import histogram
import os
import glob
import cv2
from natsort import natsorted, ns

def main():

    all_loaded_images = []

    x_scale = 1
    y_scale = 1
    z_scale = 1
    res = 10 #read every nth pixel in x,y,z

    #set path
    #all_images = "C:/Users/sgat/OneDrive - KADK/03_CITA/59_DataViz/03_data/all/"
    #all_images = "C:/Users/sgat/OneDrive - KADK/03_CITA/59_DataViz/03_data/zebicon-7a12a3/31853_2020-09-18_Wood_1_Top/"
    #all_images = "C:/Users/sgat/OneDrive - KADK/Sebastian_after_Sync/06_bones_to_flesh/browsing male/(VKH) CT Images (494 X 281)/"
    #all_images = "C:/Users/sgat/OneDrive - KADK/03_CITA/59_DataViz/03_data/31853_2020-09-18_Cellulose_2_Top/"
    #all_images = "C:/Users/sgat/OneDrive - KADK/03_CITA/59_DataViz/03_data/31853_2020-09-18_Cellulose_1_Top/"
    #all_images = "C:/Users/sgat/OneDrive - KADK/03_CITA/59_DataViz/03_data/zebicon-ccc195/31853_2020-10_02_Fungar_1_Top/"
    all_images = "C:/Users/sgat/OneDrive - KADK/03_CITA/59_DataViz/03_data/zebicon-ccc195/31853_2020-10_02_Fungar_2_Top/"
    #all_images = "C:/Users/sgat/OneDrive - KADK/03_CITA/59_DataViz/03_data/zebicon-ccc195/cropped_jpeg/"

    os.chdir(all_images)

    #load all images
    types = ('*.gif', '*.png', '*.jpg', "*bmp", "*tif")
    files = []
    for f in types:
        files.extend(glob.glob(f))
    print (len(files))

    files = natsorted(files, key=lambda y: y.lower())

    #create image stack
    img_stack = 0
    for i in range(0,len(files),res):
        if i % 100 == 0:
            print ("img loaded: ", i)
        img = cv2.imread(files[i],0)

        if res > 1:
            img = cv2.resize(img, (0,0), fx=1/res, fy=1/res)

        if i == 0:
            img.fill(0)
            img_stack = img
        else:
            img_stack = np.dstack((img_stack,img))

    img = cv2.imread(files[0],0)
    img = cv2.resize(img, (0,0), fx=1/res, fy=1/res)
    img.fill(0)
    img_stack = np.dstack((img_stack,img))
    print ("image stack created. shape:", img_stack.shape)

    #prepare data for historgam
    data = img_stack.copy()
    data = data.flatten()
    data = data [data > 5] #only data highe than 5

    #create historgam
    hist = histogram(data,
                     bins=25,
                     errors=False,
                     c = "white",
                     bc = "white",
                     ytitle= "Voxel Counts",
                     xtitle= "Density"
                    )

    #create volume
    vol = Volume(img_stack, spacing=(x_scale,y_scale,z_scale), mapper="smart", mode=1, alpha= [0.0, 0.9, 1.0], c= "coolwarm")

    #render and show
    vp = Plotter(N=2, axes=False, bg="black", size="fullscreen", sharecam=False) #(1000,1000))
    vp.show(vol, "Voxel Render", axes=0, at=0)
    vp.show(hist,"Density Historgam", axes=0, at=1)
    vp.show(interactive=True)

if __name__ == "__main__":
    main()

max_stars_repo_name: CITA-cph/deep-sight

max_issues_repo_name: CITA-cph/deep-sight

max_forks_repo_name: CITA-cph/deep-sight


text:
# library of convenience plotting funcs that are used for making various plots for all optical photostimulation/imaging experiments

# imports
import os
from typing import Union

import numpy as np
import functools
import random
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns

# global plotting params
params = {'legend.fontsize': 'x-large',
          'axes.labelsize': 'x-large',
          'axes.titlesize': 'x-large',
          'xtick.labelsize': 'x-large',
          'ytick.labelsize': 'x-large'}
plt.rcParams.update(params)
sns.set()
sns.set_style('white')

# %% UTILITY FUNCS

# wrapper for piping plots in and out of figures
def plotting_decorator(figsize=(3, 3), nrows=1, ncols=1, apply_image_frame_options=False, apply_heatmap_options=False):
    def plotting_decorator(plotting_func):
        """
        Wrapper to help simplify creating plots from matplotlib.pyplot

        :param plotting_func: plotting convenience function to be wrapped
        :return: fig and ax objects if requested, and/or shows the figure as specified

        Examples:
        ---------
        1) in this example the fig and ax will be taken directly from the kwargs inside the inner wrapper
        >>> @plotting_decorator
        >>> def example_decorated_plot(title='', **kwargs):
        >>>     fig, ax = kwargs['fig'], kwargs['ax']  # need to add atleast this line to each plotting function's definition
        >>>     print(f'|-kwargs inside example_decorated_plot definition: {kwargs}')
        >>>     ax.plot(np.random.rand(10))
        >>>     ax.set_title(title)

        >>> example_decorated_plot()

        """

        @functools.wraps(plotting_func)
        def inner(*args, **kwargs):

            return_fig_obj = False

            # set number of rows, cols and figsize
            if 'nrows' in [*kwargs]:
                nrows_ = kwargs['nrows']
            else:
                nrows_ = nrows

            if 'ncols' in [*kwargs]:
                ncols_ = kwargs['ncols']
            else:
                ncols_ = ncols

            if 'figsize' in [*kwargs]:
                figsize_ = kwargs['figsize']
            else:
                figsize_ = figsize

            # create or retrieve the fig, ax objects --> end up in kwargs to use into the plotting func call below
            if 'fig' in [*kwargs] and 'ax' in [*kwargs]:
                if kwargs['fig'] is None or kwargs['ax'] is None:
                    # print('\-creating fig, ax [1]')
                    kwargs['fig'], kwargs['axs'] = plt.subplots(nrows=nrows_, ncols=ncols_, figsize=figsize_)
                else: pass
            elif ncols_ > 1 or nrows_ > 1:
                kwargs['fig'], kwargs['axs'] = plt.subplots(nrows=nrows_, ncols=ncols_, figsize=figsize_)
            else:
                kwargs['fig'], kwargs['ax'] = plt.subplots(figsize=figsize_)

            print(f'\- executing plotting_func')
            res = plotting_func(**kwargs)  # these kwargs are the original kwargs defined at the respective plotting_func call + any additional kwargs defined in inner()

            kwargs['ax'].set_title(kwargs['title'], wrap=True) if 'title' in [*kwargs] else None

            kwargs['fig'].tight_layout(pad=1.8)

            if 'show' in [*kwargs]:
                if kwargs['show'] is True:
                    kwargs['fig'].show()
                    return res
                else:
                    if res is not None:
                        return (kwargs['fig'], kwargs['ax'], res)
                    else:
                        return (kwargs['fig'], kwargs['ax'])
            else:
                kwargs['fig'].show()
                return res

        return inner
    return plotting_decorator

def save_figure(fig, save_path_full: str = None):
    print(f'\nsaving figure to: {save_path_full}')
    os.makedirs(save_path_full)
    fig.savefig(save_path_full)

# custom colorbar for heatmaps
from matplotlib.colors import LinearSegmentedColormap
def _make_colormap(seq):
    """Return a LinearSegmentedColormap
    seq: a sequence of floats and RGB-tuples. The floats should be increasing
    and in the interval (0,1).
    """
    seq = [(None,) * 3, 0.0] + list(seq) + [1.0, (None,) * 3]
    cdict = {'red': [], 'green': [], 'blue': []}
    for i, item in enumerate(seq):
        if isinstance(item, float):
            r1, g1, b1 = seq[i - 1]
            r2, g2, b2 = seq[i + 1]
            cdict['red'].append([item, r1, r2])
            cdict['green'].append([item, g1, g2])
            cdict['blue'].append([item, b1, b2])
    return LinearSegmentedColormap('CustomMap', cdict)

# generate an array of random colors
def _get_random_color(pastel_factor=0.5):
    return [(x + pastel_factor) / (1.0 + pastel_factor) for x in [random.uniform(0, 1.0) for i in [1, 2, 3]]]

def _color_distance(c1, c2):
    return sum([abs(x[0] - x[1]) for x in zip(c1, c2)])

def _generate_new_color(existing_colors, pastel_factor=0.5):
    max_distance = None
    best_color = None
    for i in range(0, 100):
        color = _get_random_color(pastel_factor=pastel_factor)
        if not existing_colors:
            return color
        best_distance = min([_color_distance(color, c) for c in existing_colors])
        if not max_distance or best_distance > max_distance:
            max_distance = best_distance
            best_color = color
    return best_color

def make_random_color_array(n_colors):
    """
    Generates a list of random colors for an input number of colors required.

    :param n_colors: # of colors required
    :return: list of colors in RGB
    """
    colors = []
    for i in range(0, n_colors):
        colors.append(_generate_new_color(colors, pastel_factor=0.2))
    return colors

image_frame_ops = {
        'axes.spines.left': False,
        'axes.spines.bottom': False,
        'axes.spines.top': False,
        'axes.spines.right': False,
        'legend.frameon': False,
        'figure.subplot.wspace': .01,
        'figure.subplot.hspace': .01,
        'figure.figsize': (18, 13),
        'ytick.major.left': False,
        'xtick.major.bottom': False}

# Figure Style settings for notebook.
def image_frame_options(fig, ax):

    # mpl.pyplot.rcdefaults()

    # mpl.rcParams.update({
    #     'axes.spines.left': False,
    #     'axes.spines.bottom': False,
    #     'axes.spines.top': False,
    #     'axes.spines.right': False,
    #     'legend.frameon': False,
    #     'figure.subplot.wspace': .01,
    #     'figure.subplot.hspace': .01,
    #     'figure.figsize': (18, 13),
    #     'ytick.major.left': False,
    #     'xtick.major.bottom': False
    #
    # })

    #
    # mpl.rcParams.update({
    #     'axes.spines.left': False,
    #     'axes.spines.bottom': False,
    #     'axes.spines.top': False,
    #     'axes.spines.right': False,
    #     'legend.frameon': False,
    #     'figure.subplot.wspace': .01,
    #     'figure.subplot.hspace': .01,
    #     'ytick.major.left': False,
    #     'xtick.major.bottom': False
    # })

    ax.spines.left = False
    ax.spines.bottom = False
    ax.spines.top = False
    ax.spines.right = False
    fig.subplots_adjust(hspace=0.01)
    fig.subplots_adjust(wspace=0.01)
    ax.set_xticks([])
    ax.set_yticks([])
    # ax.xaxis.set_major_locator([])
    # ax.yaxis.set_major_locator([])
    # ax.xaxis.set_minor_locator([])
    # ax.yaxis.set_minor_locator([])

def dataplot_frame_options():
    import matplotlib as mpl

    mpl.rcParams.update({
        'axes.spines.top': False,
        'axes.spines.right': False,
        'legend.fontsize': 'x-large',
        'axes.labelsize': 'x-large',
        'axes.titlesize': 'x-large',
        'xtick.labelsize': 'x-large',
        'ytick.labelsize': 'x-large',
        'legend.frameon': False,
        'figure.subplot.wspace': .01,
        'figure.subplot.hspace': .01,
    })
    sns.set()
    sns.set_style('white')

def dataplot_ax_options(ax, data_length: int, **kwargs):
    """
    :param
        **kwargs:

    """
    if ax:
        ax.margins(0.1)

        # set x_axis label
        # change x-axis to time (secs) if time is requested
        if 'x_axis' in [*kwargs]:
            x_axis = kwargs['x_axis']
            if ('time' in x_axis or 'Time' in x_axis) and 'trialobj' in [*kwargs]:
                trialobj = kwargs['trialobj']

                # change x axis ticks to every 30 seconds
                labels = list(range(0, int(data_length // trialobj.imparams.fps), 30))
                ax.set_xticks(ticks=[(label * trialobj.imparams.fps) for label in labels])

                ax.set_xticklabels(labels)
                ax.set_xlabel('Time (secs)')
            else:
                ax.set_xlabel(x_axis)

        # set y_axis label
        ax.set_ylabel(kwargs['y_axis']) if 'y_axis' in [*kwargs] else None

        # set x and y axis limits
        if 'xlims' in [*kwargs]: ax.set_xlim(kwargs['xlims'])
        if 'ylims' in [*kwargs]: ax.set_ylim(kwargs['ylims'])

    else:
        pass

def heatmap_options():
    # import matplotlib as mpl

    jet = mpl.cm.get_cmap('jet')
    jet.set_bad(color='k')

# %% GENERAL PLOTTING FUNCS
### plot the location of provided coordinates
@plotting_decorator(figsize=(5, 5))
def plot_coordinates(coords: list,  frame_x: int, frame_y: int, background: np.ndarray = None, fig=None, ax=None, **kwargs):
    """
    plot coordinate locations

    :param targets_coords: ls containing (x,y) coordinates of targets to plot
    :param background: np.array on which to plot coordinates, default is black background (optional)
    :param kwargs:
    """

    if background is None:
        background = np.zeros((frame_x, frame_y), dtype='uint16')
        ax.imshow(background, cmap='gray')
    else:
        ax.imshow(background, cmap='gray')

    if 'edgecolors' in [*kwargs]:
        edgecolors = kwargs['edgecolors']
    else:
        edgecolors = 'yellowgreen'
    for (x, y) in coords:
        ax.scatter(x=x, y=y, edgecolors=edgecolors, facecolors='none', linewidths=2.0)

    ax.margins(0)
    fig.tight_layout()

    if 'title' in [*kwargs]:
        if kwargs['title'] is not None:
            ax.set_title(kwargs['title'])
        else:
            pass

max_stars_repo_name: prajayshah/funcsforprajay

max_issues_repo_name: prajayshah/funcsforprajay

max_forks_repo_name: prajayshah/funcsforprajay