text:
# -*- coding: utf-8 -*-
"""
Created on Tue Sep 29 22:59:45 2015

@author: Kevin
"""

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import os
import json
import itertools
try:
    import pykml.parser
    pykml_installed = True
except ImportError:
    pykml_installed = False

class bus:
    """
    CTA Bus Class
    =============
    Class with functions and data related to Chicago Transit Authority (CTA)
    data pertaining to bus rides broken down by day and route.

    Do this first
    -------------
    This documentation assumes you have imported and loaded up the
    objects by calling

    >>> import cta
    >>> bus = cta.bus()

    The Data
    ---------
    This class provides the data as a pandas dataframe. It is accessed by

    >>> bus.data

    With the exception of "daytype", the columns are the names of routes.
    The daytype column is a classification of the date into
        'W' : normal weekday (not holiday)
        'A' : Saturday
        'U' : Sunday or holiday
    It appears as if a holiday is considered as one of New Year's Day,
    Memorial Day, Independence Day, Labor Day, Thanksgiving, and Christmas Day

    The rest of the columns are the number of riders on the corresponding
    day and route. A value of NaN indicates that the bus was not running.
    Route 1001 is reserved for "shuttle buses used for construction or other
    unforeseen events".

    The data was downloaded as a CSV from
    https://data.cityofchicago.org/
    specifically, the CTA-Ridership-Bus-Routes-Daily-Totals-by-Route dataset.

    Of note, the data-viewer available at the supplied link seems to indicate
    that no routes have 0 ridership for any day (although quite a few have
    just 1 rider reported), but the CSV indicates otherwise.

    There are 9 routes in total that have 0 ridership reported for any day:

    >>> iszero = bus.data == 0
    >>> for c in iszero.columns:
    >>>     if any(iszero[c]):
    >>>     print(c)
    1001
    106
    168
    169
    290
    290S
    X3
    X4
    X98

    There is also a README associated with the data-set. This can be found at
    https://data.cityofchicago.org/api/assets/DBD13076-48DD-4AF8-AD91-4130456E96E9?

    Methods
    -------
    plot_rides
        Plots ridership values for specified routes as a time series
    plot_fft
        Plots the Fast Fourier Transform for specified routes
    plot_shapes
        Plots the geographic shape of the specified routes
    routes
        Returns a list of all routes
    shapes
        Returns a dict of Nx2x2 arrays of line segments describing
        the shape of each route as (lat, long) pairs.
    stops
        Returns a dict of lists of (lat, long) pairs describing
        the location of the bus stops
    condense_shape
        Returns a condensed version of an Nx2x2 array in the format
        of the values returned by shapes. The condensed version
        is a list of Nx2 arrays describing a path as (lat, long) pairs.
    routes_within_dist
        Returns a dict describing which routes are within a given
        distance from each other.

    Examples
    --------
    Plot some riderhip values.

    >>> bus.plot_rides([9, 'j14', 2, '15'])
    >>> bus.plot_rides([9, 'j14', 2, '15'], resamp='W')
    >>> bus.plot_rides([9, 'j14', 2, '15'], resamp='W',
    ...                 fillzero=True)
    >>> bus.plot_rides([9, 'j14', 2, '15'], resamp='W',
    ...                 fillzero=True, stacked=True)

    Plot the FFT for some buses.

    >>> bus.plot_fft([48, 49, 50])

    Plot the shapes of the routes. The lines are made more transparent
    inversely proportional to mean daily ridership.

    >>> bus.plot_shapes(bus.routes())

    Plot only the shapes of the routes within 50 meters of route 55.

    >>> bus.plot_shapes(bus.routes_within_dist(50, 55)['55'])
    """

    def __init__(self, filename=None):
        """
        __init__
        ========

        Loads the data from the CSV file. If filename is not given then
        ../data/bus_route_daily_totals.csv is used.
        """
        if filename is None:
            filename = os.path.join(os.path.dirname(__file__),
                                    '../data/bus_route_daily_totals.csv')
        self.data = self.__read_bus_data(filename)

    def plot_rides(self, routes, resamp=None, fillzero=False, stacked=False, ax=None):
        """
        Plots ridership for all specified routes on the same axis.

        Parameters
        ----------
        routes   : list of routes, can be integers or strings, i.e. [2, 6, 'x28']
        resamp   : String representing how to downsample for plotting. Use 'W' for
                   weekly, 'M' for monthly, 'Q' for quarterly, or 'AS' for yearly.
                   DEFAULT: resamp = None # no resampling
                   http://pandas.pydata.org/pandas-docs/stable/timeseries.html#offset-aliases
        fillzero : Set to True to fill in missing ridership data with 0. Data is
                   considered missing for any day the bus was not running.
                   This can include weekends for commuter buses, or times when
                   the bus route did not exist.
                   DEFAULT: fillzero = False
        stacked  : True to produce a stacked area plot showing total amounts as
                   well as individual route contributions.
        ax       : axes object to plot to.
                   DEFAULT: ax = None

        Returns
        -------
        ax : axes object the object was plotted on

        Notes
        -----
        Matplotlib stacked area charts do not play well with NaN values.
        If one of the specified routes is NaN at a given coordinate, then
        all routes listed *after* that route in the function call will not
        be plotted. Thus, if you have set stacked=True then it is
        recommended that you also set fillzero=True.
        """
        # create new figure if no axis provided
        if ax is None:
            fig = plt.figure()
            ax = fig.add_subplot(111)

        # be flexible with "scalar" inputs
        if not isinstance(routes, (list, np.ndarray)):
            routes = [routes]

        # resample and fill missing with 0 as necessary
        data = self.data.copy()
        if fillzero:
            data.fillna(value=0, inplace=True)
        if resamp is not None:
            data = data.resample(resamp) # mean resampling is default

        # Convert routes to strings if they are not already
        routes = map(lambda x: str(x).upper(), routes)

        # make sure route is valid
        for i in xrange(len(routes)):
            route = routes[i]
            if route not in data.columns:
                print('Invalid route ' + route + ', check bus.routes() for complete list')
                del(routes[routes.index(route)])

        if not routes: # no routes to plot
            return ax

        # do this awful mess to deal with stacked plot legends not showing up
        if stacked:
            poly_collections = ax.stackplot(data.index, data[routes].T, linewidth=0)
            legend_proxies = []
            for i, poly in enumerate(poly_collections):
                # poly.set_gid(legend_labels[i])
                legend_proxies.append(plt.Rectangle((0, 0), 1, 1, fc=poly.get_facecolor()[0]))
            ax.legend(legend_proxies, routes, ncol=max([(len(routes)/30), 1]),
                      bbox_to_anchor=[1, .5], loc='center left')
        else:
            data[routes].plot(ax=ax)

        plt.show()

        return ax

    def plot_fft(self, routes, ax=None):
        """
        Plot the [discrete/fast] fourier transform of bus riderships per day
        for all the specified routes on the optionally specified axis. The FFT
        is normalized to sum to 1 so that different routes can be compared
        more meaningfully.

        Parameters
        ----------
        routes : list of routes, can be integers or strings, i.e. [2, 6, 'x28']
                 routes can also be a scalar, i.e. just "x28" or 2
        ax     : axes object to plot to.
                 DEFAULT: ax = None

        Returns
        -------
        ax    : axes object the object was plotted on
        """
        if ax is None:
            fig = plt.figure()
            ax = fig.add_subplot(111)
        hold = ax._hold
        ax.hold(True)

        if not isinstance(routes, (list, np.ndarray)):
            routes = [routes]

        # fill missing with 0 as necessary
        data = self.data.fillna(value=0, inplace=False)

        # Convert routes to strings if they are not already
        routes = map(lambda x: str(x).upper(), routes)

        legend_labels = []
        for route in routes:

            time_x = data[route]
            y = np.absolute(np.fft.rfft(time_x))
            y = y / np.sum(y)
            x = np.fft.rfftfreq(time_x.size, 1)

            legend_labels.append(route)
            plt.plot(x, y, axes=ax, gid=route)

        plt.legend(legend_labels)
        ax.hold(hold)

        plt.show()

        return ax

    def plot_shapes(self, routes=None, transparency=True):
        '''
        Plots the shapes of the specified routes using lat/long
        pairs. The routes can have their transparency set based
        on how many people have ridden the bus since the data
        starts in 2001.

        Parameters
        ----------
        routes       : List of routes to plot. The default is to plot
                       all possible routes. Note that some routes are
                       in the daily ridership data but NOT in the
                       route shape data. This can happen because the
                       route is older than the shape data, or because
                       it is too small.
        transparency : Boolean. Set to True if you want the routes
                       to have variable transparency as described
                       above. Set to False to make all routes opaque.

        Returns
        -------
        ax : matplotlib axes object containing the plot.

        Notes
        -----
        The axis is set to have a square aspect ratio in order to
        preserve the shape of the routes. Use ax.set_aspect('auto')
        to set it back to the default value.
        '''
        plt.figure()
        ax = plt.gca()

        shapes = self.shapes()

        if routes is None:
            routes = self.routes()

        if type(routes) is not list:
            routes = [routes]

        max_ridership = self.data.fillna(0).mean().max()

        for idx, route in enumerate(routes):
            r = str(route).upper()
            if r not in shapes.keys():
                print('route ' + r + ' was not found in the shapefile.')
                continue

            if transparency:
                alpha = self.data[r].fillna(0).mean()
                alpha = alpha / max_ridership
            else:
                alpha = 1

            line_segs = shapes[r]
            cond_shape = self.condense_shape(line_segs)
            for shape in cond_shape:
                ax.plot(shape[:,1], shape[:,0], alpha=alpha, color=(.8,.9,1), linewidth=1)

        ax.set_aspect('equal', 'datalim')
        ax.set_axis_bgcolor((0,0,0))

        plt.show()

        return ax

    def routes(self):
        """
        Returns a list of routes as strings.
        """
        return self.data.columns.tolist()[1:] # [1:] gets rid of "daytype" col

    @staticmethod
    def condense_shape(shape):
        '''
        Takes in a list of line segments and condenses them by joining
        consecutive line segments using a greedy algorithm.

        Parameters
        ----------
        shape : Nx2x2 numpy array corresponding to line segments. For example,
                one of the values in the dictionary returned by shapes.

        Returns
        -------
        cond_shape : List of Mx2 arrays where line segments have been
                     joined together if they share a coordinate.

        Notes
        -----
        This method employs a simple greedy algorithm to build up the lists.
        There is no guarantee that the number of elements in the returned
        list is optimally low.
        '''
        cond_shape = []
        curr_path = shape[0].tolist()
        shape = shape[1:]

        # Simple function to test whether the shape and the path
        # share a coordinate at the specified indices
        def shares_coord(shp, pth, shp_idx, pth_idx):
            share_coords = map(lambda x: all(x),
                               shp[:,shp_idx,:] == pth[pth_idx])
            try:
                return share_coords.index(True)
            except ValueError:
                return None

        while shape.shape[0]:
            did_extend_shape = False
            for shape_index, path_index in zip([0,0,1,1],[0,-1,0,-1]):
                shared = shares_coord(shape, curr_path, shape_index, path_index)
                if shared is None:
                    continue

                # shape and curr_path have the same coordinate at
                # shape[shared,shape_index,:] and curr_path[path_index]
                # We need to append the other end of the line segment
                # in shape[shared,:,:], which is indexed by
                # shape[shared, not shape_index, :]
                shape_index = not shape_index

                if path_index == 0: # need to prepend new coordinate
                    curr_path = [shape[shared,shape_index].tolist()] + curr_path
                else: # need to append new coordinate
                    curr_path = curr_path + [shape[shared,shape_index].tolist()]

                # remove coordinate from shape
                shape = np.vstack((shape[0:shared], shape[shared+1:]))

                did_extend_shape = True
                break

            if not did_extend_shape:
                # we were not able to extend the current path
                # store the path and start a new one.
                cond_shape.append(curr_path)
                curr_path = shape[0].tolist()
                shape = shape[1:]

        return map(lambda x: np.array(x), cond_shape)

    def shapes(self, filename=None):
        '''
        Returns a set of shapes (i.e. list of (lat,lon) coordinate pairs)
        describing the bus routes. The expected input is a JSON file with
        a similar structure as ../data/CTABusRoutes.json. The code used
        to retrieve the information from the initial KML file is left in
        the function in case it needs to be re-run, but should not be
        used in general.

        Parameters
        ----------
        filename : path to the .kml file that should be parsed.
                   DEFAULT: ../data/CTABusRoutes.json.

        Returns
        -------
        route_coords : A dictionary of shapes for bus routes. The keys
                       are the bus route names and the values are Nx2x2
                       numpy arrays. If X is one of these arrays, then
                       X[i] describes a line segment such that X[i,0,:]
                       and X[i,1,:] are the two ends of the line segment.
                       The line segments are "dictionary sorted", i.e.
                       sorted along their first coordinate and then their
                       second coordinate.
        '''
        if filename is None:
            filename = os.path.join(os.path.dirname(__file__),
                                    '../data/CTABusRoutes.json')
        try:
            # Try to load the JSON file, prefered over the KML
            # file because this does not depend on a package
            # outside of the built-in python packages.
            with open(filename, 'r') as f:
                return {k: np.array(v) for k, v in json.load(f).iteritems()}
        except IOError:
            if not pykml_installed:
                print('Warning! JSON file not found and pykml is not installed!')
                return None

            # If JSON file is not there, parse the KML file.
            filename = filename[:-4] + 'kml'
            try:
                with open(filename, 'r') as f:
                    root = pykml.parser.fromstring(f.read())
            except IOError:
                print('Warning! JSON file and KML files not found')
                return None

            route_coords = {}
            for route in root.Document.Folder.iterchildren():
                if 'Placemark' in route.tag:
                    name = str(route.name)
                    line_strings = route.MultiGeometry.getchildren()
                    coords = [l.coordinates.text for l in line_strings]

                    # Convert string to numpy array
                    coords = map(lambda x: self.__parse_coords(x), coords)

                    # Change from lists of points to lists of line segments
                    # i.e., pairs of consecutive points
                    for i, arr in enumerate(coords):
                        coords[i] = [np.array(arr[j:j+2,:]) for j in range(len(arr)-1)]

                    # Flatten the list of line segments
                    coords = np.array([item for sublist in coords for item in sublist])

                    # do a dictionary sort on the coordinates
                    # facilitates removal of duplicates later.
                    for i in range(len(coords)):
                        if coords[i][0,0] < coords[i][1,0]:
                            continue
                        if coords[i][0,0] > coords[i][1,0] or coords[i][0,1] > coords[i][1,1]:
                            coords[i] = np.flipud(coords[i])

                    route_coords[name] = coords

        return route_coords

    def stops(self, filename=None):
        '''
        Returns a dict of shapes (i.e. list of (lat,lon) coordinate pairs)
        describing bus stop locations. The expected input is a JSON file
        with a similar structure as ../data/CTABusStops.json. The code
        used to retrieve the information from the initial KML file
        is left in the function in case it needs to be re-run, but
        should not be used in general.

        Parameters
        ----------
        filename : path to the .kml file that should be parsed.
                   DEFAULT: ../data/CTABusRoutes.json

        Returns
        -------
        route_stops : A dictionary of stops for bus routes. The keys
                      are the bus route names and the values are
                      Nx2 numpy arrays of [lat,long] pairs.
        '''
        if filename is None:
            filename = os.path.join(os.path.dirname(__file__),
                                    '../data/CTABusStops.json')

        try:
            # Try to load the JSON file, this should be the only
            # file source controlled in git
            # JSON file was 1/50th the size of the KML file
            with open(filename, 'r') as f:
                return {k: np.array(v) for k, v in json.load(f).iteritems()}
        except IOError:
            if not pykml_installed:
                print('Warning! JSON file not found and pykml is not installed!')
                return None

            # If JSON file is not there, parse the KML file.
            filename = filename[:-4] + 'kml'

            try:
                with open(filename, 'r') as f:
                    root = pykml.parser.fromstring(f.read())
            except IOError:
                print('Warning! JSON file and KML files not found')
                return None

            route_stops = {}
            for stop in root.Document.Folder.iterchildren():
                if 'Placemark' in stop.tag:
                    descrip = stop.description.text
                    if not 'ROUTES' in descrip:
                        print('Missing ROUTES in "' + descrip + '"')
                        continue
                    descrip = descrip[descrip.index('ROUTES')+6:]
                    descrip = descrip[:descrip.index('</tr>')]

                    routes = []
                    curr_num = ''
                    for c in descrip:
                        if c.isdigit() or c.isupper():
                            curr_num = curr_num + c
                        else:
                            if len(curr_num):
                                routes.append(curr_num)
                                curr_num = ''

                    coords = self.__parse_coords(stop.Point.coordinates.text)
                    coords = coords.tolist() # don't want these as numpy array

                    for route in routes:
                        if route in route_stops.keys():
                            route_stops[route].append(coords)
                        else:
                            route_stops[route] = [coords]
            return {k: np.array(v) for k, v in route_stops.iteritems()}

    def routes_within_dist(self, d, routes=None, accurate=False):
        """
        Returns all routes that are within the distance d.

        The distance between two routes is defined as the
        minimum distance between any two pairs of stops.

        Parameters
        ----------
        d        : The minimum distance allowed between stops, in meters
        routes   : If not None, then only routes within distance d
                   of the specified routes will be returned. If left
                   as None, then all routes will be calculated.
        accurate : Passed to lat_long_dist(). Setting accurate to True
                   increases run-time by a factor of ~4.
                   When using d=50, the majority of routes are not
                   affected (137 of 182) by the value of accurate.
                   accurate=False will always err on the side of
                   inclusion.

        Returns
        -------
        within_dist : within_dist is a dict such that if r2 is in
                      within_dist[r1], then r2 is within distance
                      d of r1.

        Notes
        -----
        The naive implementation of this algorithm would be to simply
        find the minimum distance between all pairs of stops for two
        routes, and then test if this is less than or equal to d.
        However, this would result in a great number of pairs to check.
        The search space is reduced by first finding a bounding box
        for each route (i.e. the smallest rectangle such that all
        stops are in the interior or boundary of the rectangle). We
        can quickly find the distance between two rectangles, and
        only if the distance between rectangles is <= d for a pair
        of routes do we dig deeper by comparing pairs of stops.

        TODO
        ----
        Investigate if using the reverse triangle inequality via a
        pre-computed table of distance between pairs within routes
        would result in any speedups.
        """
        if routes is None:
            routes = self.routes()
            default_routes = True
        else:
            default_routes = False
        if not type(routes) is list:
            routes = [routes]
        routes = map(lambda x: str(x).upper(), routes)

        stops = self.stops()

        all_routes = self.routes()
        all_routes = [x for x in all_routes if x in stops.keys()]

        bounding_boxes = dict()
        for route in all_routes:
            min_lat = np.min(stops[route][:,0])
            max_lat = np.max(stops[route][:,0])
            min_lon = np.min(stops[route][:,1])
            max_lon = np.max(stops[route][:,1])
            bounding_boxes[route] = mercator_projection(np.array(
                [[min_lat, min_lon],[max_lat, max_lon]]))

        # fudge factor, adjust for projection errors when creating the
        # bounding box, only relevant when accurate is True
        if accurate:
            fudge_factor = 1.1
        else:
            fudge_factor = 1

        within_dist = {k: [] for k in routes}
        already_checked = {k: [] for k in stops.keys()}
        for route1 in routes:
            # Reduces the search space as follows:
            # 1. Create bounding boxes
            # 2. Test for distance between boxes of < d
            # 3. If they pass that test, dig deeper
            if route1 not in stops.keys():
                if not default_routes:
                    print('Warning! route "' + route1 + '" has no stop locations')
                continue
            for route2 in all_routes:
                if route1 in already_checked[route2]:
                    # If it is within the distance, we will have already added it
                    continue
                already_checked[route2].append(route1)
                box_dist = box_distance(bounding_boxes[route1], bounding_boxes[route2])

                if box_dist <= fudge_factor * d:
                    curr_dist = np.inf
                    for pair in itertools.product(stops[route1].tolist(),
                                                  stops[route2].tolist()):
                        curr_dist = min(curr_dist,
                                        lat_long_dist(pair[0], pair[1],
                                                      accurate=accurate))
                        if curr_dist <= d:
                            within_dist[route1].append(route2)
                            if route2 in routes:
                                within_dist[route2].append(route1)
                            break
        return within_dist

    @staticmethod
    def __parse_coords(coords_text):
        coords = coords_text.split(' ')[1:]
        coords = map(lambda c: c.split(',')[0:2], coords)
        coords = map(lambda xy_pair: [float(c) for c in xy_pair[::-1]], coords)
        return np.array(coords)

    @staticmethod
    def __read_bus_data(filename):
        data = pd.read_csv(filename)
        data['date'] = date_lookup(data['date'])
        daytypes = data.drop_duplicates(subset='date')[['date', 'daytype']]
        daytypes = daytypes.set_index('date')
        data = data.pivot(index='date', columns='route', values='rides')
        data.insert(0,0,daytypes)
        data.rename(columns={0:'daytype'}, inplace=True)
        return data

    def __str__(self):
        route_strs = ', '.join(self.routes())
        return 'cta bus daily ridership data for routes [' + route_strs + '].'

class train:
    """
    Class with functions and data related to Chicago Transit Authority (CTA)
    data pertaining to train rides broken down by day and station.

    Warning! This class has not been looked at in quite a while. The
    focus has been on the bus class.
    """

    def __init__(self, rides_filename=None, station_filename=None):
        """
        Loads the data from CSV files.
        """
        if rides_filename is None:
            rides_filename = os.path.join(os.path.dirname(__file__),
                                          '../data/L_station_daily_entry_totals.csv')
        self.data = self.__read_train_data(rides_filename)

        if station_filename is None:
            station_filename = os.path.join(os.path.dirname(__file__),
                                            '../data/L_Stops.csv')
        self.stop_data = self.__get_L_stop_names(station_filename)

    def plot_stops(self, stops, resamp=None, fillzero=False, stacked=False, ax=None):
        """
        Plots all specified routes on the same axis.

        Parameters
        ----------
        stops    : list of stops, can be either integers corresponding to the
                   station ID (denoted MAP_ID in self.stop_data) or a string
                   corresponding to the station name. The station name will
                   attempt to be matched to the STOP_DESCRIPTIVE_NAME, STATION_NAME,
                   and STOP_NAME in that order (again, these labels are found
                   in self.stop_data).
                   See parse_stop for more information
        resamp   : String representing how to downsample for plotting. Use 'W' for
                   weekly, 'M' for monthly, 'Q' for quarterly, or 'AS' for yearly.
                   See the link below for more options.
                   http://pandas.pydata.org/pandas-docs/stable/timeseries.html#offset-aliases
                   DEFAULT: resamp = None # no resampling
        fillzero : Set to True to fill in missing ridership data with 0. Data is
                   considered missing for any day the train station was not open.
                   This can include weekends for commuter stations, or times when
                   the station did not exist.
                   DEFAULT: fillzero = False
        stacked  : True to produce a stacked area plot showing total amounts as
                   well as individual route contributions.
        ax       : axes object to plot to.
                   DEFAULT: ax = None

        Returns
        -------
        ax : axes object the object was plotted on
        """
        if ax is None:
            fig = plt.figure()
            ax = fig.add_subplot(111)

        if not isinstance(stops, (list, np.ndarray)):
            stops = [stops]

        # resample and fill missing with 0 as necessary
        data = self.data.copy()
        if fillzero:
            data.fillna(value=0, inplace=True)
        if resamp is not None:
            data = data.resample(resamp) # mean resampling is default

        # Convert routes to strings if they are not already
        stops = map(lambda x: str(x).upper(), stops)

        # make sure stop is valid
        for i in xrange(len(stops)):
            stop = stops[i]
            stop_col_name = self.parse_stop(stop)
            if stop_col_name:
                stops[i] = stop_col_name
            else:
                print('Invalid stop ' + stop + ', check *.stops() for complete list')
                del(stops[stops.index(stop)])

        if not stops:
            return ax

        # do this awful mess to deal with stacked plot legends not showing up
        if stacked:
            poly_collections = ax.stackplot(data.index, data[stops].T, linewidth=0)
            legend_proxies = []
            for i, poly in enumerate(poly_collections):
                # poly.set_gid(legend_labels[i])
                legend_proxies.append(plt.Rectangle((0, 0), 1, 1, fc=poly.get_facecolor()[0]))
            ax.legend(legend_proxies, stops, ncol=max([(len(stops)/30), 1]),
                      bbox_to_anchor=[1, .5], loc='center left')
        else:
            data[stops].plot(ax=ax)

        return ax

    def plot_fft(self, stops, ax=None):
        """
        Plot the [discrete/fast] fourier transform of train riderships per day
        for all the specified stops on the optionally specified axes object. The
        FFT is normalized so that it sums to 1. This allows different stops to be
        compared more meaningfully.

        Parameters
        ----------
        stops : list of stops, can be integers or strings, i.e. [2, 6, 'x28']
                stops can also be a scalar, i.e. just "x28" or 2
        ax    : axes object to plot to.
                DEFAULT: ax = None

        Returns
        -------
        ax    : axes object the object was plotted on
        lines : list of line objects plotted for each route.
        """
        if ax is None:
            fig = plt.figure()
            ax = fig.add_subplot(111)
        hold = ax._hold
        ax.hold(True)

        if not isinstance(stops, (list, np.ndarray)):
            stops = [stops]

        # fill missing with 0 as necessary
        data = self.data.copy()
        data.fillna(value=0, inplace=True)

        legend_labels = []
        for s in stops:
            stop = self.parse_stop(s)
            if not stop:
                print('Invalid stop ' + str(s) + ', check *.stops() for complete list')
                continue
            time_x = data[stop]
            y = np.absolute(np.fft.rfft(time_x))
            y = y / np.sum(y)
            x = np.fft.rfftfreq(time_x.size, 1)

            legend_labels.append(stop)
            plt.plot(x, y, axes=ax, gid=stop)

        plt.legend(legend_labels)
        ax.hold(hold)

        return ax

    def get_stop_by_color(self, colors):
        """
        Get a list of stations that have all the specified
        colors running through it.

        Parameters
        ----------
        colors : List of strings, i.e. 'RED', 'PURPLE'.

        Returns
        -------
        stops : List of stops that have every specified color
        """
        if type(colors) is not list:
            colors = [colors]
        colors = map(lambda x: x.upper(), colors)
        try:
            grouped = self.stop_data[['MAP_ID', 'STATION_NAME'] + colors].groupby('MAP_ID')
            grouped = grouped.max() # maximum is basically an "or" of the group
            filt = grouped[colors.pop()]
            for c in colors:
                filt = filt & grouped[c]
            data = grouped[filt].reset_index(level=0)
        except KeyError:
            print('Unexpected color, try *.stop_data.columns[7:-1] for complete list')
            return None
        return data['MAP_ID'].map(lambda x: self.parse_stop(x)).tolist()

    def parse_stop(self, stop):
        """
        Parses a stop ID or name into a standard format.

        Parameters
        ----------
        stop : list of stops, can be either integers corresponding to the
               station ID (denoted MAP_ID in self.stop_data) or
               a string corresponding to the station name. The station name will
               attempt to be matched to the STOP_DESCRIPTIVE_NAME, STATION_NAME, and
               STOP_NAME in that order (again, these labels are found
               in self.stop_data).

        Returns
        -------
        String in the format "MAP_ID: STATION_NAME" matching the column name in self.data
        """

        data = self.stop_data.copy()
        cols = self.data.columns[1:]
        col_num = map(lambda x: int(x[0:5]), cols)

        try:
            if type(stop) is str:
                if len(stop) > 5:
                    stop_num = int(stop[0:5])
                else:
                    stop_num = int(stop)
            else:
                stop_num = stop
            map_id_equal = data['MAP_ID'] == stop_num
            if any(map_id_equal):
                occurence = data[map_id_equal].iloc[0]
                return cols[col_num.index(occurence['MAP_ID'])]
        except (TypeError, ValueError):
            pass

        data['lower_case_column'] = data['STATION_DESCRIPTIVE_NAME'].map(lambda x: x.lower())
        try:
            stop_descriptive_name_equal = data['lower_case_column'] == stop.lower()
            if any(stop_descriptive_name_equal):
                occurence = data[stop_descriptive_name_equal].iloc[0]
                return cols[col_num.index(occurence['MAP_ID'])]
        except (TypeError, AttributeError, ValueError):
            pass

        data['lower_case_column'] = data['STATION_NAME'].map(lambda x: x.lower())
        try:
            station_name_equal = data['lower_case_column'] == stop.lower()
            if any(station_name_equal):
                occurence = data[station_name_equal].iloc[0]
                return cols[col_num.index(occurence['MAP_ID'])]
        except (TypeError, AttributeError, ValueError):
            pass

        data['lower_case_column'] = data['STOP_NAME'].map(lambda x: x.lower())
        try:
            stop_name_equal = data['lower_case_column'] == stop.lower()
            if any(stop_name_equal):
                occurence = data[stop_name_equal].iloc[0]
                return cols[col_num.index(occurence['MAP_ID'])]
        except (TypeError, AttributeError, ValueError):
            pass

        return False

    def stops(self):
        """
        Returns a list of all of the available stops.
        """
        return self.data.columns.tolist()[1:]

    @staticmethod
    def __read_train_data(filename):
        data = pd.read_csv(filename)
        data['date'] = date_lookup(data['date'])

        # Combine the ID and station-name fields into one column
        ids = data['station_id'].map(lambda x: str(x))
        names = data['station_name']
        id_names = ids + ': ' + names
        data.drop(['station_id', 'station_name'], inplace=True, axis=1)
        data['station_id_name'] = id_names

        # Save the daytypes
        daytypes = data.drop_duplicates(subset='date')[['date', 'daytype']]
        daytypes = daytypes.set_index('date')

        # Pivot, insert the daytypes back in
        data = data.pivot_table(index='date', columns='station_id_name', values='rides')
        data.insert(0,0,daytypes)
        data.rename(columns={0:'daytype'}, inplace=True)
        return data

    @staticmethod
    def __get_L_stop_names(filename):
        data = pd.read_csv(filename)

        # Convert this from string "(x,y)" to actual tuple
        data['Location'] = data['Location'].map(lambda x:
            (float(x[1:-1].split(',')[0]), float(x[1:-1].split(',')[1])))

        # Rename those awful default names
        data.rename(columns={'G':'GREEN', 'BRN':'BROWN', 'P':'PURPLE',
                             'Pexp':'PINK_EXPRESS', 'Y':'YELLOW',
                             'Pnk':'PINK', 'O':'ORANGE'},
                    inplace=True)

        return data

# Utility Functions

def date_lookup(s):
    """
    This is an extremely fast approach to datetime parsing.
    For large data, the same dates are often repeated. Rather than
    re-parse these, we store all unique dates, parse them, and
    use a lookup to convert all dates.

    Thanks to fixxxer, found at
    http://stackoverflow.com/questions/29882573
    """
    dates = {date:pd.to_datetime(date) for date in s.unique()}
    return s.apply(lambda v: dates[v])

def lat_long_dist(x, y, accurate=False, metric='euclidean'):
    """
    Computes the distance between two (lat,long) coordinates.

    Parameters
    ----------
    x        : lat/long pair number 1, in decimal degrees
    y        : lat/long pair number 2, in decimal degrees
    accurate : True to use a highly accurate calculation
               of the distance between the two points.
               If accurate is False, then the mercator
               projection will be used. You can expect
               a relative error of about 0.5%.
    metric   : The distance metric to use. Can be either
               'euclidean' or 'manhattan', the latter is
               only valid if accurate=False.

    Returns
    -------
    d : the distance between x and y in meters

    Sources
    -------
    [1] http://www.movable-type.co.uk/scripts/latlong.html
    [2] https://en.wikipedia.org/wiki/Earth_radius

    ToDo
    ----
    Compute the manhattan distance in the accurate version.
    Should probably complete in this order:
        1. Find the rectangle assuming Earth is a sphere
        2. Find the rectangle assuming Earth is ellipsoid
    With 2 being much harder than 1 (probably).
    """

    if accurate:
        if not metric == 'euclidean':
            raise ValueError('Invalid metric argument with accurate=True')
        phi1  = x[0] * np.pi / 180.0
        phi2  = y[0] * np.pi / 180.0
        delta_phi    = phi2 - phi1
        delta_lambda = (y[1] - x[1]) * np.pi / 180.0

        eq_rad  = 6378137.0 # radius at equator in meters
        pol_rad = 6356752.3 # radius at poles in meters
        R = np.sqrt(
        ((eq_rad**2 * np.cos(phi1))**2 + (pol_rad**2 * np.sin(phi1))**2) /
        ((eq_rad * np.cos(phi1))**2 + (pol_rad * np.sin(phi1))**2))
        # R is the radius of the eart at latitude phi1

        a = np.sin(delta_phi)**2.0 + \
            np.cos(phi1) * np.cos(phi2) * np.sin(delta_lambda) ** 2.0

        c = np.arctan2(a ** .5, (1 - a)**.5)

        d = R * c
    else:
        delta_lat = (y[0] - x[0]) * 110574.0
        delta_lon = (y[1] - x[1]) * 111320.0 * np.cos(x[0])

        if metric == 'euclidean':
            d = np.sqrt(delta_lat**2 + delta_lon**2)
        elif metric == 'manhattan':
            d = abs(delta_lat) + abs(delta_lon)
        else:
            raise ValueError('Unknown distance metric')

    return d

def mercator_projection(lat_longs, phi=0.730191653):
    """
    Compute the mercator projection of the Nx2 matrix lat_longs

    Parameters
    ----------
    lat_longs : Nx2 matrix of (lat,long) pairs
    phi       : The central latitude for the projection,
                in radians.
                DEFAULT: 0.730191653, the latitude of Chicago

    Returns
    -------
    projected : Nx2 matrix of cartesian coordinates in
                the mercator projection. Units are meters.
    """
    lat_longs = np.array(lat_longs)

    lat_longs[:,0] = lat_longs[:,0] * 110574.0
    lat_longs[:,1] = lat_longs[:,1] * 111320.0 * np.cos(phi)

    return lat_longs

def box_distance(box1, box2):
    """
    Computes the minimum distance between two boxes, i.e.
    rectangles whose sides are parallel to the x-y axes.

    Parameters
    ----------
    box1 : 2x2 array of coordinates of the first rectangle
           Should be the lower left and upper right coord
    box2 : 2x2 array of coordinates of the second rectangle
           Should be the lower left and upper right coord

    Returns
    -------
    dist : The minimum distance between the boxes
    """
    x1 = box1[0,0]
    y1 = box1[0,1]
    x1b = box1[1,0]
    y1b = box1[1,1]

    x2 = box2[0,0]
    y2 = box2[0,1]
    x2b = box2[1,0]
    y2b = box2[1,1]

    # See http://stackoverflow.com/questions/4978323
    left = x2b < x1
    right = x1b < x2
    bottom = y2b < y1
    top = y1b < y2
    if top and left:
        return np.sqrt((x1 - x2b) ** 2 + (y1b - y2) ** 2)
    elif left and bottom:
        return np.sqrt((x1 - x2b) ** 2 + (y1 - y2b) ** 2)
    elif bottom and right:
        return np.sqrt((x1b - x2) ** 2 + (y1 - y2b) ** 2)
    elif right and top:
        return np.sqrt((x1b - x2) ** 2 + (y1b - y2) ** 2)
    elif left:
        return x1 - x2b
    elif right:
        return x2 - x1b
    elif bottom:
        return y1 - y2b
    elif top:
        return y2 - y1b
    else: # rectangles overlap
        return 0

max_stars_repo_name: kbrose/chicago-data

max_issues_repo_name: kbrose/dataViz

max_forks_repo_name: kbrose/chicago-data


text:
#! /usr/bin/env python
# coding: utf-8
"""
This is Ros node for pwm control rc car
"""
from enum import Enum
from std_msgs.msg import Bool
from sensor_msgs.msg import Imu
from geometry_msgs.msg import Twist, TwistStamped
from rc_car_msgs.msg import CarParams, CarPwmContol
import math
import numpy as np
import rospy
import time

from geometry_msgs.msg import TwistStamped
import tf
import matplotlib.pyplot as plt  # creation plots
import copy

# ---------------------------------------------------------------------------------------------------------------------
# ---------- Глобальные переменные

# ---------------------------------------------------------------------------------------------------------------------
# ---------- Функции подписки на топики

norm_velocity = 0.0
norm_acceleration = 0.0
current_pwm = CarPwmContol()

srez = -100
# ---------------------------------------------------------------------------------------------------------------------
# ---------- Вспомогательные функции
def velocity_clb(data):
    """
    Get current velocity from FCU
    :param data: velocity from NED
    """
    global norm_velocity
    norm_velocity = np.linalg.norm(np.array([data.twist.linear.x,
                                             data.twist.linear.y]))

def imu_clb(data):
    """
    Get current velocity from FCU
    :param data: velocity from NED
    """
    global norm_acceleration
    norm_acceleration = np.linalg.norm(np.array([data.linear_acceleration.x,
                                                 data.linear_acceleration.y]))

def pmw_motor_clb(data):
    """
    Get current velocity from FCU
    :param data: velocity from NED
    """
    global current_pwm
    current_pwm = data

# ---------------------------------------------------------------------------------------------------------------------
# ---------- Главная функция
def main():
    """
    Основной цикл узла ROS.

    @return: код завершения программы
    """
    # Инициализируем узел ROS
    rospy.init_node("rc_reg_plotter")

    # Подписываемся на топики
    rospy.Subscriber("/pwm_output", CarPwmContol, pmw_motor_clb)
    rospy.Subscriber("/mavros/imu/data", Imu, imu_clb)
    rospy.Subscriber("mavros/local_position/velocity", TwistStamped, velocity_clb)

    # Инициализация "слушателя" навигационных преобразований
    tf_listener = tf.TransformListener()

    # Инициализируем спискиpwm_output для построения графиков
    velocity_list = list()
    accel_list = list()
    pwm_list = list()
    time_plot = list()

    old_ros_time = time.time()
    # Основной цикл
    rate = rospy.Rate(10)
    t_ = 0.0
    while not rospy.is_shutdown():
        # Сичтываем навигационные данные

        # Рассчитываем время
        dt = time.time() - old_ros_time
        old_ros_time = time.time()

        t_ += dt
        velocity_list.append(norm_velocity)
        accel_list.append(norm_acceleration)
        pwm_list.append(current_pwm.MotorPWM)
        time_plot.append(t_)
	if srez > 0:
		velocity_list = velocity_list[srez::]
		accel_list = accel_list[srez::]
		pwm_list = pwm_list[srez::]
		time_plot = time_plot[srez::]

        # print ('t %s' % t_)

        # Рисуем графики
        ax1 = plt.subplot2grid((2,1), (0,0))
        ax1.plot(time_plot, velocity_list, label='velocity')
        ax1.plot(time_plot, accel_list, label='accel')
        ax1.set_ylabel('x, m/s')
        ax1.set_xlabel('t, s')
        ax1.legend()
        ax1.grid()

        ax2 = plt.subplot2grid((2,1), (1,0))
        ax2.plot(time_plot, pwm_list, label='pwm')
        ax2.set_title('PWM')
        ax2.set_ylabel('y, ms')
        ax2.set_xlabel('t, s')
        ax2.legend()
        ax2.grid()

        # ax3 = plt.subplot2grid((3, 3), (0, 2))
        # ax3.plot(time_plot, drone_z, label='drone z')
        # ax3.plot(time_plot, goal_z, label='goal z')
        # ax3.set_title('Coords: drone Z and goal Z')
        # ax3.set_ylabel('z, m')
        # ax3.set_xlabel('t, s')
        # ax3.legend()
        # ax3.grid()

        plt.pause(1.0 / 10.0)

        # rate.sleep()
    return 0

if __name__ == '__main__':
    try:
        main()
    except rospy.ROSInterruptException:
        pass

ext: py

reason: import numpy


text:
import numpy as np
import cvxpy as cp
import torch
import torch.nn as nn
import torch.autograd as autograd
import torch.optim as optim
import torch.nn.functional as F
from cvxpylayers.torch import CvxpyLayer

def solve(price, amb, T, pn, a1, a3, max_theta, min_theta, max_power, theta_0, tensor=False):
    p = cp.Variable(T)
    y = cp.Variable(T)
    theta = cp.Variable(T)
    obj = price@p + pn * cp.sum_squares(y)
    cons = [y >= 0, p >= 0]
    for i in range(T):
        if i == 0:
            cons.append(theta[i] == a1*theta_0 + (1-a1)*amb[i] - a3*p[i] )
        else:
            cons.append(theta[i] == a1*theta[i-1] + (1-a1)*amb[i] - a3*p[i] )
        cons.append(p[i] <= max_power)
        cons.append(theta[i] <= max_theta + y[i])
        cons.append(theta[i] >= min_theta - y[i])
    problem = cp.Problem(cp.Minimize(obj), cons)
    problem.solve()
    if tensor:
        return torch.tensor(p.value),torch.tensor(theta.value),torch.tensor(y.value)
    return p.value, theta.value, y.value

def model(price_data, amb_data, theta_0, T = 24*5):
    i = 0

    pn_p = cp.Parameter(1,nonneg=True)
    a1_p = cp.Parameter(1)
    a3_p = cp.Parameter(1)
    max_theta_p = cp.Parameter(1)
    min_theta_p = cp.Parameter(1)
    max_power_p = cp.Parameter(1)

    # problem formulation
    price = price_data[i*T:(i+1)*T]
    amb = amb_data[i*T:(i+1)*T]

    p_v = cp.Variable(T)
    y_v = cp.Variable(T)
    theta_v = cp.Variable(T)
    obj = price@p_v + pn_p * cp.sum_squares(y_v)
    cons = [y_v >= 0, p_v >= 0]
    for i in range(T):
        if i == 0:
            cons.append(theta_v[i] == a1_p*theta_0 + (1-a1_p)*amb[i] - a3_p*p_v[i])
        else:
            cons.append(theta_v[i] == a1_p*theta_v[i-1] + (1-a1_p)*amb[i] - a3_p*p_v[i])
        cons.append(p_v[i] <= max_power_p)
        cons.append(theta_v[i] <= max_theta_p + y_v[i])
        cons.append(theta_v[i] >= min_theta_p - y_v[i])
    problem = cp.Problem(cp.Minimize(obj), cons)
    assert problem.is_dpp()
    layer = CvxpyLayer(problem, [pn_p,a1_p,a3_p,max_theta_p,min_theta_p,max_power_p], [p_v,theta_v,y_v])
    return layer

reason: import numpy,import cvxpy,from cvxpy


text:
# Copyright (c) 2021 ralabo.jp
# This software is released under the MIT License.
# see https://opensource.org/licenses/mit-license.php
# ====================================================
import os
import sys
import re
import json
import numpy as np
import pandas as pd
from copy import deepcopy
from collections import defaultdict

from scipy.stats import hmean
from natsort import natsorted, natsort_keygen
from sklearn import metrics
from glob import glob

from utils.score import gamma_ppf
from common import com, TargetDomain
from flow.flow_provider import get_flow_instance
from flow.transfer_learn import switch_training_param

class IterFiles(object):
    def __init__(self):
        pass

    def result_dirs(self,csv_files):
        """
        result_dir が同じファイルを集める

        return result_dir, selected_files
        """
        selected_files = []
        prev_result_dir = None
        for x in csv_files:
            result_dir = x[0] # get dir
            if result_dir != prev_result_dir:
                if selected_files:
                    yield prev_result_dir, selected_files
                    selected_files = []

            selected_files += [x]
            prev_result_dir = result_dir

        if selected_files:
            yield prev_result_dir, selected_files

    def matchine_types(self,csv_files):
        """
        return machine_type, selected_files
        """
        selected_files = []
        prev_machine_type = None
        for x in csv_files:
            basename = x[1] # get base
            machine_type = basename.split('_')[2] # machine_type
            if machine_type != prev_machine_type:
                if selected_files:
                    yield prev_machine_type, selected_files
                    selected_files = []

            selected_files += [x]
            prev_machine_type = machine_type

        if selected_files:
            yield prev_machine_type, selected_files

    def target_dirs(self,csv_files):
        """
        return target_dir, selected_files
        target_dir: eg. 'source_test', 'target_test'
        """
        selected_files = []
        prev_dir_name = None
        for x in csv_files:
            basename = x[1] # get base
            target_dir = '_'.join(basename.split('_')[5:7]) # target_dir
            if target_dir != prev_dir_name:
                if selected_files:
                    yield prev_dir_name, selected_files
                    selected_files = []
            selected_files += [x]
            prev_dir_name = target_dir

        if selected_files:
            yield prev_dir_name, selected_files

    def section_names(self,csv_files):
        """
        return section_name, selected_files
        section_name: eg. 'section_00', 'section_01', ...
        """
        selected_files = []
        prev_section_name = None
        for x in csv_files:
            basename = x[1] # get base
            section_name = '_'.join(basename.split('_')[3:5]) # section_name
            if section_name != prev_section_name:
                if selected_files:
                    yield prev_section_name, selected_files
                    selected_files = []
            selected_files += [x]
            prev_section_name = section_name

        if selected_files:
            yield prev_section_name, selected_files

def get_thresh(machine_type, section_index, target, decision_thr, flow):
    domain = flow.domain(machine_type, section_index, target)
    score_distr_param_path = com.score_distr_param_path(domain, target=target)
    with open(score_distr_param_path) as f:
        gamma_params = json.load(f)

    thresh = gamma_ppf(decision_thr, gamma_params)
    assert not np.isnan(thresh), 'thresh is nan.'\
        ' gamma param=(q={q} a={a} loc={loc} scale={scale} thresh={thresh} file={file})'\
        .format(q=decision_thr, **gamma_params, file=score_distr_param_path)
    return thresh, gamma_params

def make_result_files(combinations):
    """
    _result_*.csv を収集し result_files を作成する

    yamlで指定された result_directory:
        result_directory: <somepath>/{config}/{feature}/{net}/{fit}/{decision}
    にある
    anomaly_score_path()関数で定義されるファイル:
        anormaly_score_{machine_type}_section_{section_name}_{target_dir}.csv
    を収集し、スコアを計算する
    スコアは dataframe として返す

    return: result_path_params
        a list of tuples, (result_path, param, transfer_learning)
    """
    com.save_param() # 作業前にcom.paramを保存する

    # target独立(転移学習ではない)の結果を得る
    csv_files = make_result_file_common(combinations, transfer_learning=False)

    # target依存(転移学習)の結果を得る
    csv_files += make_result_file_common(combinations, transfer_learning=True)

    com.restore_param() # 作業後に保存したcom.paramを復帰

    assert csv_files,'no _result_*.csv files.'
    return csv_files

def make_result_file_common(combinations, transfer_learning):
    """
    transfer_learning:
        trueなら転移学習の結果ファイルをリストアップする
        falseなら通常学習の結果ファイルをリストアップする
    return:
        a list of tuples, [(result_path, param.copy(), transfer_learning),...]
    """
    if transfer_learning:
        # 転移学習結果を探索するためのワイルドカードを作成する
        tl_result_directory = com.param.get('transfer_learning') and \
                              com.param['transfer_learning'].get('result_directory')
        if tl_result_directory:
            # true時、転移学習用paramに切り替える
            fit_key, eval_key = switch_training_param(transfer_learning)
        else:
            return []

    # 探索するファイルのワイルドカードを作成する
    result_dir = com.subst_model_config(com.param['result_directory'])
    result_dir = re.sub(r'{\w+}','*',result_dir)
    anomaly_score_csv = com.anomaly_score_path(machine_type='*',
                                               #section_name='*', # exclude 03, 04, 05
                                               section_name='section_0[012]',
                                               target_dir='*',
                                               result_directory=result_dir)
    # glob でファイルを探索する
    csv_files = glob(anomaly_score_csv)
    csv_files = [x.replace('\\','/') for x in csv_files] # 区切りを'/'にそろえる
    csv_files = sorted(csv_files) # 順番を並べ替える
    # ファイル名を (dir, base, ext) に分解
    csv_files = [(os.path.dirname(x),)+tuple(os.path.splitext(os.path.basename(x)))
                 for x in csv_files]
    if combinations:
        # フィルターする
        def _filter_fun(x):
            for c in combinations: # 組み合わせの中を探索する
                if com.is_match_combination(x,c): # 条件がマッチするとき
                    return True

        csv_files = list(filter(_filter_fun, csv_files ))

    # ディレクトリ名ごとにスコアを計算する
    iter = IterFiles()
    result_path_params = []
    for result_dir, files_1 in iter.result_dirs(csv_files):
        com.restore_param() # 上位で保存していたcom.paramを復帰
        # transfer_learning==true時、転移学習用paramに切り替える
        fit_key, eval_key = switch_training_param(transfer_learning)
        # result_dirからハイパーパラメータを取得する
        cond_tuples = com.get_cond_tuples_from_dirname(result_dir)
        com.update_param_from_cond_tuples(cond_tuples) # update com.param
        if '_base_param' in com.param:
            # cond_tuplesのkeyとcom.param[key]の組を辞書base_param_とする
            base_param_ = {k:com.param[k] for k,v in cond_tuples.items()}
            # 上記辞書base_param_で com.param['_base_param']を更新する
            com.param['_base_param'].update(base_param_)
            # com.param['_base_param']['cond']を更新する
            com.update_cond_params(com.param['_base_param'])

        param = com.param # update local param
        # get decision parameters
        decision_thr = param['decision']['decision_thr']
        max_fpr = param['decision']['max_fpr']
        # get flow for model_for
        flow = get_flow_instance(param)
        # machine_types > target_dirs > section_names  ごとに集計する
        # initialize lines in csv for AUC and pAUC
        csv_lines = []
        performance_over_all = [] if com.mode else None
        for machine_type, files_2 in iter.matchine_types(files_1):
            # determine decision threshold
            # load anomaly score distribution for training
            if com.mode:
                # results for each machine type
                csv_lines.append([machine_type])
                csv_lines.append(["section", "target", "AUC", "pAUC", "precision",
                                  "recall", "F1 score"])
                performance = []
            for target_dir, files_3 in iter.target_dirs(files_2):
                for section_name,files_4 in iter.section_names(files_3):
                    assert len(files_4) == 1 , 'not unique,  check source code'

                    section_index = section_name.split('_')[1] # eg. '00', '01', ...
                    target = target_dir.split('_')[0] # eg. 'source' or 'target'

                    # get threshold
                    thresh, gamma_params = get_thresh(machine_type, section_index, target,
                                                      decision_thr, flow)

                    # read anomaly_score_csv
                    dir, base, ext = files_4[0]
                    anomaly_score_csv = dir + '/' + base + ext # compose a full path
                    df = pd.read_csv(anomaly_score_csv,header=None)
                    df.columns = ['wav','a_score']
                    if com.mode:
                        # get y_true from wav_filename
                        # section_00_target_test_normal_0000.wav, anormaly_score
                        # 正解はファイル名('wav'欄)中の'anomaly'か'normal'かで取得する
                        y_true = [1 if x.split('_')[4]=='anomaly' else 0 for x in df.loc[:,'wav']]
                        # 異常スコア('a_score'欄)を取得する
                        y_pred = df.loc[:,'a_score'] # anomaly score

                        # append AUC and pAUC to lists
                        auc = metrics.roc_auc_score(y_true, y_pred)
                        p_auc = metrics.roc_auc_score(y_true, y_pred,
                                                      max_fpr=max_fpr)
                        tn, fp, fn, tp = metrics.confusion_matrix(
                            y_true, [1 if x > thresh
                                     else 0 for x in y_pred]).ravel()
                        prec = tp / np.maximum(tp + fp, sys.float_info.epsilon)
                        recall = tp / np.maximum(tp + fn, sys.float_info.epsilon)
                        f1 = 2.0 * prec * recall / np.maximum(prec + recall,
                                                              sys.float_info.epsilon)
                        section_index_ = section_name.split("_", 1)[1]
                        target_ = target_dir.split("_", 1)[0]
                        csv_lines.append([section_index_, # 'section'
                                          target_ + ' ' + section_index_, # 'target'
                                          auc, p_auc, prec, recall, f1])
                        performance.append([auc, p_auc, prec, recall, f1])
                        performance_over_all.append([auc, p_auc, prec, recall, f1])
            # end for dirname
            if com.mode:
                # calculate averages for AUCs and pAUCs
                amean_performance = np.mean(np.array(performance, dtype=float), axis=0)
                csv_lines.append(["arithmetic mean", ""] + list(amean_performance))
                hmean_performance = hmean(np.maximum(np.array(performance, dtype=float),
                                                     sys.float_info.epsilon), axis=0)
                csv_lines.append(["harmonic mean", ""] + list(hmean_performance))
                csv_lines.append([])
        # end for matchie_type
        if com.mode:
            csv_lines.append(["", "", "AUC", "pAUC", "precision", "recall", "F1 score"])
            # calculate averages for AUCs and pAUCs
            amean_performance = np.mean(np.array(performance_over_all, dtype=float), axis=0)
            csv_lines.append(["arithmetic mean over all machine types, sections, and domains",
                                ""] + list(amean_performance))
            hmean_performance = hmean(np.maximum(np.array(performance_over_all, dtype=float),
                                                 sys.float_info.epsilon), axis=0)
            csv_lines.append(["harmonic mean over all machine types, sections, and domains",
                                ""] + list(hmean_performance))
            csv_lines.append([])

            # output results
            result_basename = param['result_file'].format(**cond_tuples)
            result_path = result_dir + os.sep + result_basename
            com.to_csv(csv_lines, result_path)
            com.print(2,'save',os.path.basename(result_path))
            param_ = deepcopy(param)
            # check param_ consisttency
            if param_.get('_base_param'):
                assert param_['model_for'] == param_['_base_param']['model_for']
                assert str(tuple(v for k,v in param_['model_for'].items())) \
                       == param_['cond']['model_for']
                assert str(tuple(v for k,v in param_['model_for'].items())) \
                       == param_['_base_param']['cond']['model_for']
            # append new row
            result_path_params += [(result_path, param_, transfer_learning)]

    return result_path_params

def collect_results(combinations):
    result_path_params = make_result_files(combinations)

    # ファイル内容をデータフレームに読み込む
    col_names = ['section','target','AUC','pAUC','precision','recall','F1 score']
    rows=[]
    for result_path_, param_, transfer_learning_ in result_path_params:

        if param_.get('_base_param'):
            assert param_['model_for'] == param_['_base_param']['model_for']

        # 拡張子なしのベースネーム
        basename = os.path.splitext(os.path.basename(result_path_))[0]
        #print('reading',basename)

        # csv を読む
        df = pd.read_csv(result_path_,names=col_names,error_bad_lines=False)

        # machin_type と　basename 欄を加える
        df.insert(0,'machine_type','') # insert column
        df.insert(0,'basename','') # insert column
        df.insert(0,'result_path_','') # insert column
        df.insert(0,'param_','') # insert column

        # 行を加工する
        for idx,row in df.iterrows():
            section = row['section']
            if section=='section' or pd.isna(section):
                pass
            elif 'mean' in section or section.isdigit():
                if 'over all' in section:
                    machine_type = '00 over all'

                row['basename']=basename
                row['machine_type']=machine_type
                row['section']=section
                row['result_path_']=result_path_
                row['param_']=param_
                rows.append(row)
            elif section:
                machine_type = section

    df = pd.DataFrame(rows)

    return df

def normalize_results(df):
    """
    # #####################################
    # 扱いやすくするための後処理
    # #####################################
    """
    # 行をソートする
    df = df.sort_values(['machine_type','target','section','basename'],key=natsort_keygen())

    # sectionの変化する点に空行を入れる(グラフを見やすくするため)
    out_rows = []
    prev_section = None
    prev_machine_type = None
    col_names = df.columns
    for idx,row in df.iterrows():

        # section間に空行を挿入
        section = row['section']
        if pd.isna(section):
           pass
        elif section != prev_section:
            out_rows.append(pd.Series(dtype=object))
            prev_section = section

        # machine_type 間にヘッダー行を挿入
        machine_type = row['machine_type']
        if pd.isna(machine_type):
           pass
        elif machine_type != prev_machine_type:
            header = row.copy()
            header.values[:] = col_names
            out_rows.append(header)
            prev_machine_type = machine_type

        # 空白の'target'(target or souce) に'section'を簡略し転記する
        if pd.isna(row['target']):
            if pd.isna(section):
                pass
            else:
                if section.startswith('arithmetic mean over all'):
                    row['target'] = 'a. m. all '
                elif section.startswith('harmonic mean over all'):
                    row['target'] = 'h. m. all '
                elif section == 'arithmetic mean':
                    row['target'] = 'a. m. '+row['machine_type']
                elif section == 'harmonic mean':
                    row['target'] = 'h. m. '+row['machine_type']

        out_rows.append(row)

    out_df = pd.DataFrame(out_rows)

    return out_df

def insert_mean_score_each_machine_type(df, score_key='AUC'):
    """
    マシンタイプ別の平均スコアAUC列を追加する(pd.DataFrame.assign()を使用)

    """
    # dfを舐めて、追加するデータを集める

    # arithmetic mean用の辞書を用意する
    # harmonic mean 用の辞書も用意する
    mean_data = dict()
    mean_data['a_mean'] = defaultdict(list)
    mean_data['h_mean'] = defaultdict(list)

    # 条件記入用に5列分のkeyを作成する
    note_keys = ['{}{}'.format(score_key, j) for j in range(5,0,-1)]

    # 最初にnote_keysを辞書に登録する
    for key in note_keys:
        mean_data['a_mean'][key] = None
        mean_data['h_mean'][key] = None

    # dataframeを調べて各machineのarithmetic meanとharmonic meanを追加する
    for idx, row in df.iterrows():
        section = row['section']
        if not pd.isna(section) and section.startswith('arithmetic mean'):
            key = row['machine_type'] + ' ' + score_key
            mean_data['a_mean'][key].append(row[score_key]) # arithmetic mean of score
        if not pd.isna(section) and section.startswith('harmonic mean'):
            key = row['machine_type'] + ' ' + score_key
            mean_data['h_mean'] [key].append(row[score_key]) # harmonic mean of score

    # 実験条件数(行数)を求める
    num_conds = len(mean_data['a_mean'][key]) # keyは最後のkeyを使う

    # 条件コラムに空白行をセットする
    for key in note_keys:
        mean_data['a_mean'][key] = [np.nan] * num_conds # 空白行
        mean_data['h_mean'][key] = [np.nan] * num_conds # 空白行

    # 集めたデータが4行目からに来るように空行とヘッダーを挿入する
    assign_data = defaultdict(list)
    for key in mean_data['a_mean'].keys():
        data = [np.nan] # 本当のヘッダーの次の空行
        data += [key] # ヘッダー(kがヘッダー文字である)を追加
        data += mean_data['a_mean'][key] # arithmetic meanを追加
        data += [key] # 個別ヘッダー行を入れる
        data += mean_data['h_mean'][key] # harmonic mean を追加
        data += [np.nan]*(len(df)-len(data)) # 残りの行を空行で埋める
        assign_data[key] = data # 列 k を更新する

    # 列を追加する
    df = df.assign(**assign_data)

    return df

def add_anomaly_distr_param(df):
    """
    異常スコア分布パラメータを追記する
    """

    # index をユニークに更新する
    df = df.reset_index(drop=True)
    for idx,row in df.iterrows():
        machine_type = row['machine_type']
        basename = row['basename']
        basename = basename if isinstance(basename,str) and basename.startswith('_result_') else None
        section = row['section']
        section = section if isinstance(section,str) and section[0].isdigit() else None

        if df.loc[idx,'basename'] == 'basename':
            prev_head_idx = idx # 最後の見つかったヘッダー行の位置

        # 個別 machine_type, section, domain の行であるとき
        if basename and section and machine_type != '00 over all':
            param_ = row['param_']

            if param_.get('_base_param'):
                assert param_['model_for'] == param_['_base_param']['model_for']
                assert str(tuple(v for k,v in param_['model_for'].items())) \
                       == param_['cond']['model_for']
                assert str(tuple(v for k,v in param_['model_for'].items())) \
                       == param_['_base_param']['cond']['model_for']

            com.update_param(param_) # update com.param
            # 分布パラメータファイルパスを取得する(flow->domainを求め取得)
            target = row['target'].split()[0] # eg. row['target'] looks like 'source 00'
            section_index = section
            flow = get_flow_instance(com.param)
            domain = flow.domain(machine_type, section_index, target)
            score_distr_param_path = com.score_distr_param_path(domain)

            # 分布パラメータを取得する
            with open(score_distr_param_path) as f:
                gamma_param = json.load(f)

            # 分布パラメータを追記する
            for k in gamma_param.keys():
                df.loc[idx,k] = gamma_param[k]

            if pd.isna(df.loc[prev_idx,'basename']):# 前行が空行の時
                for k in gamma_param.keys():
                    if pd.isna(df.loc[prev_head_idx, k]): # 未記入
                        df.loc[prev_head_idx, k] = k
                    else: # すでに記入済み
                        break # 他も記入済みのはずなのでループを終了

        prev_idx = idx

    return df

def add_cond_tuples(df):
    """
    列の最後に新しい列を追加する
    'basename'の条件をtupleに分離して
     可能ならばサブパラメータに分解して列の最後に追加する
    """
    # index をユニークなものにする
    df = df.reset_index(drop=True)
    new_columns = defaultdict(dict)
    first = True
    for idx, row in df.copy().iterrows():# df はループ内で更新されることに注意
        machine_type = row['machine_type']

        # 個別 machine_type の平均部分であるとき
        basename = row['basename']
        if not pd.isna(basename) and basename != 'basename':
            com.update_param(row['param_']) # update com.param
            # 'result_path_'から条件を cond_tuples に抽出する
            result_dir = os.path.dirname(row['result_path_'])
            cond_tuples = com.get_cond_tuples_from_dirname(result_dir)

            # 条件名をキーとする階層辞書を作成する
            # cond_params::= dict(cond, cond_value or dict(subparam, subparam_value) )
            cond_params =  {}
            for k,v in cond_tuples.items():
                if v.startswith('('): # is v tuple ?
                    # make a dict(subparam, subparam_value)
                    sub_param = com.get_cond_sub_params(cond_tuples, k)
                    if k.startswith('tl_'):
                        # sub_paramのkeyにプレフィックス'tl_'を付ける
                        sub_param = {'tl_'+k : v for k,v in sub_param.items()}
                    cond_params[k] = sub_param
                else: # v is a value
                    cond_params[k] = v

            # 新しい列を追加する(全体のヘッダーを自動作成)
            if first: # 最初のみ
                kwds = {} # 追加する列を{kwd:value}形式で作成
                rows = [np.nan]*len(df) # 追加する行(空白)
                for k1, v1 in cond_params.items():
                    if isinstance(v1, dict):
                        for k2 in v1.keys():
                            kwds[k2] = rows
                    else:
                        kwds[k1] = rows
                df = df.assign(**kwds) # 新しい列を追加
                first = False

            # 現在行idxに値を記入する
            for k1, v1 in cond_params.items():
                if isinstance(v1, dict):
                    for k2, v2 in v1.items():
                        df.loc[idx, k2] = str(v2)
                else:
                    df.loc[idx, k1] = str(v1)

            # 前行がヘッダー行であれば前行にヘッダーを追記する
            if df.loc[prev_idx,'basename']=='basename': # prev is a header
                for k1, v1 in cond_params.items():
                    if isinstance(v1, dict):
                        for k2, v2 in v1.items():
                            df.loc[prev_idx, k2] = k2
                    else:
                        df.loc[prev_idx, k1] = k1

        prev_idx = idx

    return df

def main(combinations, save_to='csv', **options):
    """
    すべての処理をする
    """
    com.parse_command_line(**options)

    # #####################
    # result_csv を収集する
    # #####################
    df = collect_results(combinations) # result_csv を収集、dataframeを取得
    df = normalize_results(df) # 扱いやすくするため、情報を追加
    df = insert_mean_score_each_machine_type(df, score_key='AUC') # マシンタイプ別の平均スコア列を追加
    df = insert_mean_score_each_machine_type(df, score_key='pAUC') # マシンタイプ別の平均スコア列を追加
    df = add_anomaly_distr_param(df) # 異常スコア分布パラメータを追記する(machine_type=='00 over all')
    df = add_cond_tuples(df) # 条件tupleを列に追加

    del df['result_path_']
    del df['param_']

    if save_to == 'ods':
        summary_ods = com.summary_ods_path() # 保存先
        with pd.ExcelWriter(summary_ods, engine='odf') as writer:
            df.to_excel(writer, index=False)
        com.print(2,'save',summary_ods)
    elif save_to == 'xlsx':
        summary_xlsx = com.summary_xlsx_path() # 保存先
        with pd.ExcelWriter(summary_xlsx) as writer:
            df.to_excel(writer, index=False)
        com.print(2,'save',summary_xlsx)
    elif save_to == 'csv':
        summary_csv = com.summary_csv_path() # 保存先
        df.to_csv(summary_csv, index=False)
        com.print(2,'save',summary_csv)

if __name__=='__main__':
    sys.argv += re.sub("\s+"," ","""
    --dev
    --defaults ./run/ave/ave_config.yaml
    --machine_types ToyCar
    --nunits 4096 2048 1024 512
    --n_frames 5
    --epochs 20
    --kld_weight 1.0
    --lr 0.0005
    --zdim 256
    --edim 32
    """).split() # split()は先頭と末尾の空白は無視するようだ

    main(combinations=None)

    if 0: # debug only
        import os
        os.system(r"../../../../Temp/dcase2021_task2_pytorch/ave/debug/_summary_debug.csv".replace('/','\\'))

reason: import numpy,from scipy


text:
import os
import sys
import numpy as np
import torch
import torch.nn.functional as F
from torch import nn
BASE_DIR = os.path.dirname(os.path.abspath(__file__))
sys.path.append(os.path.join(BASE_DIR, 'sampling'))
from data import Tree
from sampling import furthest_point_sample

class PointNet(nn.Module):
    def __init__(self, feat_len):
        super(PointNet, self).__init__()

        self.conv1 = nn.Conv1d(3, 64, 1)
        self.conv2 = nn.Conv1d(64, 128, 1)
        self.conv3 = nn.Conv1d(128, 1024, 1)

        self.mlp = nn.Linear(1024, feat_len)

    def forward(self, x):
        x = x.permute(0, 2, 1)

        x = F.leaky_relu(self.conv1(x))
        x = F.leaky_relu(self.conv2(x))
        x = F.leaky_relu(self.conv3(x))

        x = x.max(dim=-1)[0]

        x = F.leaky_relu(self.mlp(x))

        return x

class PartEncoder(nn.Module):

    def __init__(self, feat_len):
        super(PartEncoder, self).__init__()

        self.conv1 = nn.Conv1d(3, 64, 1)
        self.conv2 = nn.Conv1d(64, 128, 1)
        self.conv3 = nn.Conv1d(128, 128, 1)
        self.conv4 = nn.Conv1d(128, feat_len, 1)

        self.fc = nn.Linear(feat_len, feat_len)

    def forward(self, pc):
        net = pc.transpose(2, 1)

        net = F.leaky_relu(self.conv1(net))
        net = F.leaky_relu(self.conv2(net))
        net = F.leaky_relu(self.conv3(net))
        net = F.leaky_relu(self.conv4(net))

        net = net.max(dim=2)[0]

        net = F.leaky_relu(self.fc(net))

        return net

class ChildrenEncoder(nn.Module):

    def __init__(self, feat_len, hidden_len, symmetric_type):
        super(ChildrenEncoder, self).__init__()

        print(f'Using Symmetric Type: {symmetric_type}')
        self.symmetric_type = symmetric_type

        self.child_op = nn.Linear(feat_len + Tree.num_sem, hidden_len)
        self.second = nn.Linear(hidden_len, feat_len)

    def forward(self, child_feats):
        child_feats = F.leaky_relu(self.child_op(child_feats))

        if self.symmetric_type == 'max':
            parent_feat = child_feats.max(dim=1)[0]
        elif self.symmetric_type == 'sum':
            parent_feat = child_feats.sum(dim=1)
        elif self.symmetric_type == 'avg':
            parent_feat = child_feats.mean(dim=1)
        else:
            raise ValueError(f'Unknown symmetric type: {self.symmetric_type}')

        parent_feat = F.leaky_relu(self.second(parent_feat))

        return parent_feat

class Network(nn.Module):

    def __init__(self, args, device):
        super(Network, self).__init__()
        self.args = args
        self.device = device

        self.part_encoder = PartEncoder(args.feat_len)
        self.children_encoder = ChildrenEncoder(args.feat_len, args.hidden_len, args.symmetric_type)

        self.pointnet_encoder = PointNet(args.feat_len)

        self.mlp1 = nn.Linear(args.feat_len, 1)
        self.mlp2 = nn.Linear(args.feat_len, 1)

        self.final_activation = args.final_activation

    def encode(self, node, leaf_feats):
        if len(node.children) == 0:
            return leaf_feats[:, node.geo_id]
        else:
            child_feats = []
            batch_size = leaf_feats.shape[0]
            for cnode in node.children:
                cur_child_feat = torch.cat([self.encode(cnode, leaf_feats), cnode.get_semantic_one_hot().repeat(batch_size, 1)], dim=1)
                child_feats.append(cur_child_feat.unsqueeze(dim=1))
            child_feats = torch.cat(child_feats, dim=1)

            return self.children_encoder(child_feats)

    def forward(self, node, pcs):
        batch_size = pcs.shape[0]
        part_cnt = pcs.shape[1]

        ### StructureNet Discriminator
        # encode all leaf geometry
        leaf_feats = self.part_encoder(pcs.view(batch_size * part_cnt, -1, 3)).view(batch_size, part_cnt, -1)

        # encode the hierarchy until root
        sn_net = self.encode(node, leaf_feats)

        ### holistic PointNet Discriminator
        # fps
        shape_pcs = pcs.view(batch_size, -1, 3)
        with torch.no_grad():
            shape_pc_id1 = torch.arange(batch_size).unsqueeze(1).repeat(1, self.args.num_point_per_shape).long().view(-1).to(self.device)
            shape_pc_id2 = furthest_point_sample(shape_pcs, self.args.num_point_per_shape).long().view(-1)
        shape_pcs = shape_pcs[shape_pc_id1, shape_pc_id2].view(batch_size, self.args.num_point_per_shape, 3)

        # pointnet encode
        pn_net = self.pointnet_encoder(shape_pcs)

        ### balance two discriminators
        # [SN Dis] produce real/fake scores
        sn_net = self.mlp1(sn_net).squeeze(1)
        if self.final_activation:
            sn_net = torch.sigmoid(sn_net)

        # [PN Dis] produce real/fake scores
        pn_net = self.mlp2(pn_net).squeeze(1)
        if self.final_activation:
            pn_net = torch.sigmoid(pn_net)

        # aggregate two scores
        net = sn_net + pn_net * self.args.pointnet_dis_score_multiplier

        return net, sn_net, pn_net


text:
"""
This module implements the Sequential Least SQuares Programming optimization
algorithm (SLSQP), originally developed by Dieter Kraft.
See http://www.netlib.org/toms/733

Functions
---------
.. autosummary::
   :toctree: generated/

    approx_jacobian
    fmin_slsqp

"""

from __future__ import division, print_function, absolute_import

__all__ = ['approx_jacobian','fmin_slsqp']

import numpy as np
from scipy.optimize._slsqp import slsqp
from numpy import zeros, array, linalg, append, asfarray, concatenate, finfo, \
                  sqrt, vstack, exp, inf, where, isfinite, atleast_1d
from .optimize import wrap_function, OptimizeResult, _check_unknown_options

__docformat__ = "restructuredtext en"

_epsilon = sqrt(finfo(float).eps)

def approx_jacobian(x,func,epsilon,*args):
    """
    Approximate the Jacobian matrix of a callable function.

    Parameters
    ----------
    x : array_like
        The state vector at which to compute the Jacobian matrix.
    func : callable f(x,*args)
        The vector-valued function.
    epsilon : float
        The perturbation used to determine the partial derivatives.
    args : sequence
        Additional arguments passed to func.

    Returns
    -------
    An array of dimensions ``(lenf, lenx)`` where ``lenf`` is the length
    of the outputs of `func`, and ``lenx`` is the number of elements in
    `x`.

    Notes
    -----
    The approximation is done using forward differences.

    """
    x0 = asfarray(x)
    f0 = atleast_1d(func(*((x0,)+args)))
    jac = zeros([len(x0),len(f0)])
    dx = zeros(len(x0))
    for i in range(len(x0)):
        dx[i] = epsilon
        jac[i] = (func(*((x0+dx,)+args)) - f0)/epsilon
        dx[i] = 0.0

    return jac.transpose()

def fmin_slsqp(func, x0, eqcons=(), f_eqcons=None, ieqcons=(), f_ieqcons=None,
               bounds=(), fprime=None, fprime_eqcons=None,
               fprime_ieqcons=None, args=(), iter=100, acc=1.0E-6,
               iprint=1, disp=None, full_output=0, epsilon=_epsilon,
               callback=None):
    """
    Minimize a function using Sequential Least SQuares Programming

    Python interface function for the SLSQP Optimization subroutine
    originally implemented by Dieter Kraft.

    Parameters
    ----------
    func : callable f(x,*args)
        Objective function.  Must return a scalar.
    x0 : 1-D ndarray of float
        Initial guess for the independent variable(s).
    eqcons : list, optional
        A list of functions of length n such that
        eqcons[j](x,*args) == 0.0 in a successfully optimized
        problem.
    f_eqcons : callable f(x,*args), optional
        Returns a 1-D array in which each element must equal 0.0 in a
        successfully optimized problem.  If f_eqcons is specified,
        eqcons is ignored.
    ieqcons : list, optional
        A list of functions of length n such that
        ieqcons[j](x,*args) >= 0.0 in a successfully optimized
        problem.
    f_ieqcons : callable f(x,*args), optional
        Returns a 1-D ndarray in which each element must be greater or
        equal to 0.0 in a successfully optimized problem.  If
        f_ieqcons is specified, ieqcons is ignored.
    bounds : list, optional
        A list of tuples specifying the lower and upper bound
        for each independent variable [(xl0, xu0),(xl1, xu1),...]
        Infinite values will be interpreted as large floating values.
    fprime : callable `f(x,*args)`, optional
        A function that evaluates the partial derivatives of func.
    fprime_eqcons : callable `f(x,*args)`, optional
        A function of the form `f(x, *args)` that returns the m by n
        array of equality constraint normals.  If not provided,
        the normals will be approximated. The array returned by
        fprime_eqcons should be sized as ( len(eqcons), len(x0) ).
    fprime_ieqcons : callable `f(x,*args)`, optional
        A function of the form `f(x, *args)` that returns the m by n
        array of inequality constraint normals.  If not provided,
        the normals will be approximated. The array returned by
        fprime_ieqcons should be sized as ( len(ieqcons), len(x0) ).
    args : sequence, optional
        Additional arguments passed to func and fprime.
    iter : int, optional
        The maximum number of iterations.
    acc : float, optional
        Requested accuracy.
    iprint : int, optional
        The verbosity of fmin_slsqp :

        * iprint <= 0 : Silent operation
        * iprint == 1 : Print summary upon completion (default)
        * iprint >= 2 : Print status of each iterate and summary
    disp : int, optional
        Over-rides the iprint interface (preferred).
    full_output : bool, optional
        If False, return only the minimizer of func (default).
        Otherwise, output final objective function and summary
        information.
    epsilon : float, optional
        The step size for finite-difference derivative estimates.
    callback : callable, optional
        Called after each iteration, as ``callback(x)``, where ``x`` is the
        current parameter vector.

    Returns
    -------
    out : ndarray of float
        The final minimizer of func.
    fx : ndarray of float, if full_output is true
        The final value of the objective function.
    its : int, if full_output is true
        The number of iterations.
    imode : int, if full_output is true
        The exit mode from the optimizer (see below).
    smode : string, if full_output is true
        Message describing the exit mode from the optimizer.

    See also
    --------
    minimize: Interface to minimization algorithms for multivariate
        functions. See the 'SLSQP' `method` in particular.

    Notes
    -----
    Exit modes are defined as follows ::

        -1 : Gradient evaluation required (g & a)
         0 : Optimization terminated successfully.
         1 : Function evaluation required (f & c)
         2 : More equality constraints than independent variables
         3 : More than 3*n iterations in LSQ subproblem
         4 : Inequality constraints incompatible
         5 : Singular matrix E in LSQ subproblem
         6 : Singular matrix C in LSQ subproblem
         7 : Rank-deficient equality constraint subproblem HFTI
         8 : Positive directional derivative for linesearch
         9 : Iteration limit exceeded

    Examples
    --------
    Examples are given :ref:`in the tutorial <tutorial-sqlsp>`.

    """
    if disp is not None:
        iprint = disp
    opts = {'maxiter': iter,
            'ftol': acc,
            'iprint': iprint,
            'disp': iprint != 0,
            'eps': epsilon,
            'callback': callback}

    # Build the constraints as a tuple of dictionaries
    cons = ()
    # 1. constraints of the 1st kind (eqcons, ieqcons); no jacobian; take
    #    the same extra arguments as the objective function.
    cons += tuple({'type': 'eq', 'fun': c, 'args': args} for c in eqcons)
    cons += tuple({'type': 'ineq', 'fun': c, 'args': args} for c in ieqcons)
    # 2. constraints of the 2nd kind (f_eqcons, f_ieqcons) and their jacobian
    #    (fprime_eqcons, fprime_ieqcons); also take the same extra arguments
    #    as the objective function.
    if f_eqcons:
        cons += ({'type': 'eq', 'fun': f_eqcons, 'jac': fprime_eqcons,
                  'args': args}, )
    if f_ieqcons:
        cons += ({'type': 'ineq', 'fun': f_ieqcons, 'jac': fprime_ieqcons,
                  'args': args}, )

    res = _minimize_slsqp(func, x0, args, jac=fprime, bounds=bounds,
                          constraints=cons, **opts)
    if full_output:
        return res['x'], res['fun'], res['nit'], res['status'], res['message']
    else:
        return res['x']

def _minimize_slsqp(func, x0, args=(), jac=None, bounds=None,
                    constraints=(),
                    maxiter=100, ftol=1.0E-6, iprint=1, disp=False,
                    eps=_epsilon, callback=None,
                    **unknown_options):
    """
    Minimize a scalar function of one or more variables using Sequential
    Least SQuares Programming (SLSQP).

    Options
    -------
    ftol : float
        Precision goal for the value of f in the stopping criterion.
    eps : float
        Step size used for numerical approximation of the jacobian.
    disp : bool
        Set to True to print convergence messages. If False,
        `verbosity` is ignored and set to 0.
    maxiter : int
        Maximum number of iterations.

    """
    _check_unknown_options(unknown_options)
    fprime = jac
    iter = maxiter
    acc = ftol
    epsilon = eps

    if not disp:
        iprint = 0

    # Constraints are triaged per type into a dictionnary of tuples
    if isinstance(constraints, dict):
        constraints = (constraints, )

    cons = {'eq': (), 'ineq': ()}
    for ic, con in enumerate(constraints):
        # check type
        try:
            ctype = con['type'].lower()
        except KeyError:
            raise KeyError('Constraint %d has no type defined.' % ic)
        except TypeError:
            raise TypeError('Constraints must be defined using a '
                            'dictionary.')
        except AttributeError:
            raise TypeError("Constraint's type must be a string.")
        else:
            if ctype not in ['eq', 'ineq']:
                raise ValueError("Unknown constraint type '%s'." % con['type'])

        # check function
        if 'fun' not in con:
            raise ValueError('Constraint %d has no function defined.' % ic)

        # check jacobian
        cjac = con.get('jac')
        if cjac is None:
            # approximate jacobian function.  The factory function is needed
            # to keep a reference to `fun`, see gh-4240.
            def cjac_factory(fun):
                def cjac(x, *args):
                    return approx_jacobian(x, fun, epsilon, *args)
                return cjac
            cjac = cjac_factory(con['fun'])

        # update constraints' dictionary
        cons[ctype] += ({'fun': con['fun'],
                         'jac': cjac,
                         'args': con.get('args', ())}, )

    exit_modes = {-1: "Gradient evaluation required (g & a)",
                    0: "Optimization terminated successfully.",
                    1: "Function evaluation required (f & c)",
                    2: "More equality constraints than independent variables",
                    3: "More than 3*n iterations in LSQ subproblem",
                    4: "Inequality constraints incompatible",
                    5: "Singular matrix E in LSQ subproblem",
                    6: "Singular matrix C in LSQ subproblem",
                    7: "Rank-deficient equality constraint subproblem HFTI",
                    8: "Positive directional derivative for linesearch",
                    9: "Iteration limit exceeded"}

    # Wrap func
    feval, func = wrap_function(func, args)

    # Wrap fprime, if provided, or approx_jacobian if not
    if fprime:
        geval, fprime = wrap_function(fprime, args)
    else:
        geval, fprime = wrap_function(approx_jacobian, (func, epsilon))

    # Transform x0 into an array.
    x = asfarray(x0).flatten()

    # Set the parameters that SLSQP will need
    # meq, mieq: number of equality and inequality constraints
    meq = sum(map(len, [atleast_1d(c['fun'](x, *c['args'])) for c in cons['eq']]))
    mieq = sum(map(len, [atleast_1d(c['fun'](x, *c['args'])) for c in cons['ineq']]))
    # m = The total number of constraints
    m = meq + mieq
    # la = The number of constraints, or 1 if there are no constraints
    la = array([1, m]).max()
    # n = The number of independent variables
    n = len(x)

    # Define the workspaces for SLSQP
    n1 = n + 1
    mineq = m - meq + n1 + n1
    len_w = (3*n1+m)*(n1+1)+(n1-meq+1)*(mineq+2) + 2*mineq+(n1+mineq)*(n1-meq) \
            + 2*meq + n1 + ((n+1)*n)//2 + 2*m + 3*n + 3*n1 + 1
    len_jw = mineq
    w = zeros(len_w)
    jw = zeros(len_jw)

    # Decompose bounds into xl and xu
    if bounds is None or len(bounds) == 0:
        xl = np.empty(n, dtype=float)
        xu = np.empty(n, dtype=float)
        xl.fill(np.nan)
        xu.fill(np.nan)
    else:
        bnds = array(bounds, float)
        if bnds.shape[0] != n:
            raise IndexError('SLSQP Error: the length of bounds is not '
                             'compatible with that of x0.')

        bnderr = bnds[:, 0] > bnds[:, 1]
        if bnderr.any():
            raise ValueError('SLSQP Error: lb > ub in bounds %s.' %
                             ', '.join(str(b) for b in bnderr))
        xl, xu = bnds[:, 0], bnds[:, 1]

        # Mark infinite bounds with nans; the Fortran code understands this
        infbnd = ~isfinite(bnds)
        xl[infbnd[:, 0]] = np.nan
        xu[infbnd[:, 1]] = np.nan

    # Clip initial guess to bounds (SLSQP may fail with bounds-infeasible initial point)
    have_bound = np.isfinite(xl)
    x[have_bound] = np.clip(x[have_bound], xl[have_bound], np.inf)
    have_bound = np.isfinite(xu)
    x[have_bound] = np.clip(x[have_bound], -np.inf, xu[have_bound])

    # Initialize the iteration counter and the mode value
    mode = array(0,int)
    acc = array(acc,float)
    majiter = array(iter,int)
    majiter_prev = 0

    # Print the header if iprint >= 2
    if iprint >= 2:
        print("%5s %5s %16s %16s" % ("NIT","FC","OBJFUN","GNORM"))

    while 1:

        if mode == 0 or mode == 1:  # objective and constraint evaluation requird

            # Compute objective function
            fx = func(x)
            try:
                fx = float(np.asarray(fx))
            except (TypeError, ValueError):
                raise ValueError("Objective function must return a scalar")
            # Compute the constraints
            if cons['eq']:
                c_eq = concatenate([atleast_1d(con['fun'](x, *con['args']))
                                     for con in cons['eq']])
            else:
                c_eq = zeros(0)
            if cons['ineq']:
                c_ieq = concatenate([atleast_1d(con['fun'](x, *con['args']))
                                     for con in cons['ineq']])
            else:
                c_ieq = zeros(0)

            # Now combine c_eq and c_ieq into a single matrix
            c = concatenate((c_eq, c_ieq))

        if mode == 0 or mode == -1:  # gradient evaluation required

            # Compute the derivatives of the objective function
            # For some reason SLSQP wants g dimensioned to n+1
            g = append(fprime(x),0.0)

            # Compute the normals of the constraints
            if cons['eq']:
                a_eq = vstack([con['jac'](x, *con['args'])
                               for con in cons['eq']])
            else:  # no equality constraint
                a_eq = zeros((meq, n))

            if cons['ineq']:
                a_ieq = vstack([con['jac'](x, *con['args'])
                                for con in cons['ineq']])
            else:  # no inequality constraint
                a_ieq = zeros((mieq, n))

            # Now combine a_eq and a_ieq into a single a matrix
            if m == 0:  # no constraints
                a = zeros((la, n))
            else:
                a = vstack((a_eq, a_ieq))
            a = concatenate((a,zeros([la,1])),1)

        # Call SLSQP
        slsqp(m, meq, x, xl, xu, fx, c, g, a, acc, majiter, mode, w, jw)

        # call callback if major iteration has incremented
        if callback is not None and majiter > majiter_prev:
            callback(x)

        # Print the status of the current iterate if iprint > 2 and the
        # major iteration has incremented
        if iprint >= 2 and majiter > majiter_prev:
            print("%5i %5i % 16.6E % 16.6E" % (majiter,feval[0],
                                               fx,linalg.norm(g)))

        # If exit mode is not -1 or 1, slsqp has completed
        if abs(mode) != 1:
            break

        majiter_prev = int(majiter)

    # Optimization loop complete.  Print status if requested
    if iprint >= 1:
        print(exit_modes[int(mode)] + "    (Exit mode " + str(mode) + ')')
        print("            Current function value:", fx)
        print("            Iterations:", majiter)
        print("            Function evaluations:", feval[0])
        print("            Gradient evaluations:", geval[0])

    return OptimizeResult(x=x, fun=fx, jac=g[:-1], nit=int(majiter),
                          nfev=feval[0], njev=geval[0], status=int(mode),
                          message=exit_modes[int(mode)], success=(mode == 0))

if __name__ == '__main__':

    # objective function
    def fun(x, r=[4, 2, 4, 2, 1]):
        """ Objective function """
        return exp(x[0]) * (r[0] * x[0]**2 + r[1] * x[1]**2 +
                            r[2] * x[0] * x[1] + r[3] * x[1] +
                            r[4])

    # bounds
    bnds = array([[-inf]*2, [inf]*2]).T
    bnds[:, 0] = [0.1, 0.2]

    # constraints
    def feqcon(x, b=1):
        """ Equality constraint """
        return array([x[0]**2 + x[1] - b])

    def jeqcon(x, b=1):
        """ Jacobian of equality constraint """
        return array([[2*x[0], 1]])

    def fieqcon(x, c=10):
        """ Inequality constraint """
        return array([x[0] * x[1] + c])

    def jieqcon(x, c=10):
        """ Jacobian of Inequality constraint """
        return array([[1, 1]])

    # constraints dictionaries
    cons = ({'type': 'eq', 'fun': feqcon, 'jac': jeqcon, 'args': (1, )},
          {'type': 'ineq', 'fun': fieqcon, 'jac': jieqcon, 'args': (10,)})

    # Bounds constraint problem
    print(' Bounds constraints '.center(72, '-'))
    print(' * fmin_slsqp')
    x, f = fmin_slsqp(fun, array([-1, 1]), bounds=bnds, disp=1,
                      full_output=True)[:2]
    print(' * _minimize_slsqp')
    res = _minimize_slsqp(fun, array([-1, 1]), bounds=bnds,
                          **{'disp': True})

    # Equality and inequality constraints problem
    print(' Equality and inequality constraints '.center(72, '-'))
    print(' * fmin_slsqp')
    x, f = fmin_slsqp(fun, array([-1, 1]),
                      f_eqcons=feqcon, fprime_eqcons=jeqcon,
                      f_ieqcons=fieqcon, fprime_ieqcons=jieqcon,
                      disp=1, full_output=True)[:2]
    print(' * _minimize_slsqp')
    res = _minimize_slsqp(fun, array([-1, 1]), constraints=cons,
                          **{'disp': True})

max_stars_repo_name: SKPsanjeevi/scipy

max_issues_repo_name: SKPsanjeevi/scipy

max_forks_repo_name: SKPsanjeevi/scipy

reason: import numpy,from numpy,from scipy