text:
import sqlite3 as sq
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import pandas as pd
import itertools

class HBdataMonitor():
    '''
    Python class to extract data from a SQLite HB database
    requirements
    Python      = 3.6
    numpy       = 1.15.0
    pandas      = 0.23.4
    matplotlib  = 2.2.2
    sqlite      = 3.24.0
    '''

    def __init__(self):
        ## default inputVar is zero to check wheter get_DB_info has runned
        self.inputVar = 0
        ## colors
        self.colors = []
        ## linestyle
        self.linestyle = plt.Line2D.filled_markers

#    def read_sql_list(self, path, sql_command):
#        # connect withe the myTable database
#        connection = sq.connect(path)
#
#        # cursor object
#        crsr = connection.cursor()
#
#        # execute the command to fetch all the data from the table emp
#        crsr.execute(sql_command)
#
#        # store all the fetched data in the ans variable
#        data= crsr.fetchall()
#        return data

    def read_sql_table(self, path, command):
        '''
        Execute a SQL command and returns the output as DataFrame
        :param str path: path of the databse
        :param str command: SQL command
        :return: DataFrame table: table with the results
        '''

        db = sq.connect(path)
        table = pd.read_sql_query(command, db)
        return table

    def read_variables(self,path, columnId, inputvar=True):
        '''
        Returns the sorted unique values for a given variable.
        :param str path: path of the databse
        :param int columnId: ColumnID of the variabele
        :param logical inputvar: switch between input or output variable (True =  inputvariables)
        :return: array var: array with the unique sorted variables
        '''

        if inputvar==True:
            string = 'Input'
        else:
            string = 'Result'
        ## command
        command1 = 'SELECT Value '
        command2 = 'FROM HydroDynamic{}Data '.format(string)
        command3 = 'INNER JOIN HydroDynamicData '
        command4 = 'ON HydroDynamicData.HydroDynamicDataId = HydroDynamic{}Data.HydroDynamicDataId '.format(string)
        command5 = 'WHERE HRD{}columnId={} '.format(string, columnId)
        ## total command
        total_command = command1 + command2 + command3 + command4 + command5
        ## execute command
        values = self.read_sql_table(path , total_command)
        ## unique and sort
        var = np.sort(np.unique(values))
        return var

    def get_variables(self, path):
        '''
        Get the values of the variables
        :param str path: path of the databse
        '''
        ## Read input/output variables
        inputVar = self.read_sql_table(path, 'SELECT * FROM HRDInputVariables')
        resultVar = self.read_sql_table(path, 'SELECT * FROM HRDResultVariables')

        ## link columnname with columndID
        self.inputVar = {}
        for i, item in enumerate(inputVar.ColumnName):
            self.inputVar[item] = inputVar.HRDInputColumnId[i]
        self.resultVar = {}
        for i, item in enumerate(resultVar.ColumnName):
            self.resultVar[item] = resultVar.HRDResultColumnId[i]

        ## input Variabele
        self.inputValues = {}
        for i, item in enumerate(inputVar.ColumnName):
            self.inputValues[item] = self.read_variables(path, inputVar.HRDInputColumnId[i], inputvar=True)
        ## winddirection
        self.inputValues['WindDir'] = np.squeeze(
            self.read_sql_table(path, 'SELECT Direction FROM HRDWindDirections').values)
        ## location
        self.inputValues['loc'] = np.squeeze(self.read_sql_table(path, 'SELECT Name FROM HRDLocations').values)
        ## Closing situations
        self.inputValues['ClosingSituation'] = np.squeeze(
            self.read_sql_table(path, 'SELECT ClosingSituationId FROM closingSituations').values, axis=1)

    def get_DB_info(self,path, fname='Variables.txt'):
        '''
        Read the info from table general.
        :param str path: path of the databse
        :param str fname: name of the asci file
        '''

        ## show general information
        command = "SELECT * FROM General"
        table = self.read_sql_table(path ,command)
        for i, item in enumerate(table.columns):
            print('{}: {}'.format(item, table.iloc[0,i]))
        ## get the variables
        self.get_variables(path)

        ## show input variables
        print('Number of input Variables={}'.format(len(self.inputValues)))
        for i, item in enumerate(self.inputValues.keys()):
            print(item)

        ## show ouput variables
        print('Number of result Variables={}'.format(len(self.resultVar)))
        for i, item in enumerate(self.resultVar):
            print(item)

        ## write variables to file
        with open(fname, 'w') as f:
            f.write('InputVariables\n')
            ## loop through variables
            for row, item in enumerate(self.inputValues.keys()):
                f.write('{};'.format(item))
                ## loop through values
                string = ''
                for col, item in enumerate(self.inputValues[item]):
                    string = string + '{};'.format(item)
                f.write(string + '\n')

    def HBoutput(self,Y,inputX,fname='result.txt'):
        '''
        input function to execute sql command
        1. create a dictionary which can be used to extract the data
        2. extract the data for every combinations
        3. write data to file
        :param str Y: result variabele
        :param dict inputX: dictionary with the variabele and value in a list {'WINDS':[0],'MEERP':[1,2]}
        :param fname: name of asci file
        '''

        ## when get_get_DB_info is not used
        if self.inputVar ==0:
            ## only first databse is selected!
            self.get_variables(inputX['DB'][0])

        ## 1 Determine constant and varying variables
        var_keys = {}
        con_keys = {}
        for item in inputX.keys():
            ## use all input variables
            if inputX[item] == ['all']:
                inputX[item] = self.inputValues[item]
            if len(inputX[item])>1:
                var_keys[item] = len(inputX[item])
            else:
               con_keys[item] = inputX[item]

        ## 2. checks
        for item in self.inputValues.keys():
            ## check if all variables are given
            assert item in inputX.keys(), 'Error: variable {} is missing'.format(item)
            ## check if values exist
            if not all(x in self.inputValues[item] for x in inputX[item]):
                print('Warning: Not all values are present in the database')

        ## 3. create new dict with empty lists
        input_data = {}
        for item in inputX.keys():
            input_data[item] = []

        ## 4. get all combinations of the variables
        var = list(var_keys.keys())
        comb_input = []
        for item in var:
            comb_input.append(inputX[item])
        combinations = list(itertools.product(*comb_input))
        N = len(combinations)

        ## 5. fill dict with all combinations
        for i in range(N):
            for item in inputX.keys():
                ## variables which vary
                if item in var_keys.keys():
                    for k,_ in enumerate(var):
                        if item==var[k]:
                            tmp = input_data[item]
                            tmp.append(combinations[i][k])
                            input_data[item] = tmp
                ## variable which are constant
                else:
                   tmp = input_data[item]
                   tmp.append(inputX[item][0])
                   input_data[item] = tmp

        ## 6. Write data and execute SQL command
        f = open(fname,'w')
        f.write('constants \n')
        for item in con_keys.keys():
            f.write('{}={}\n'.format(item,con_keys[item][0]))
        f.write('varied \n')
        f.write('{};'.format(Y))
        for item in var:
            f.write('{};'.format(item))
        f.write('\n')
        for i in range(N):
            ##create temporary variable dict of the inputVariables
            tmp = {}
            for item in inputX.keys():
                ## WindDir, loc, DB and closingsituation are not inputVariables
                if item!='WindDir' and item!='loc' and item!='ClosingSituation' and item!='DB':
                    tmp[item] = input_data[item][i]
            y = self.create_sql_command(input_data['DB'][i], Y, input_data['loc'][i], input_data['WindDir'][i], input_data['ClosingSituation'][i], tmp, fig=False)
            if y.empty:
                y.set_value(0, -999)

            if i%10==0:
                print('progress: {} %'.format(np.round(float(i)/N * 100),2))
            ## write data
            for j in range(len(y)):
                f.write('{};'.format(y.iloc[j]))
                for item in var:
                    f.write('{};'.format(input_data[item][i]))
                f.write('\n')
        f.close()

    def create_plot(self,Y=0, X=1, Z1=0, Z2=0, fname='result.txt', path='result'):
        '''
        Create a plot from the data points
        :param str Y: variabele on y-axis
        :param str X: variabele on x-axis
        :param str Z1: variabele with different symbols
        :param str Z2: variabele with different colors
        :param str fname: name of asci file
        :param str path: location where the figure is saved
        '''

        ## get data information
        f = open(fname,'r')
        meta = f.readlines()
        f.close()
        constant = ''
        for i,line in enumerate(meta):
            if line[0:6] == 'varied':
                n = i
                break
            if i>0:
                constant = constant + line
        ## get data
        data = pd.read_csv(fname,sep=';',skiprows=n+1,na_values=-999)

        ## checks
        assert len(data.columns)<=5, 'Error: too many columns to plot'
        if len(data.columns)==5:
            assert Z2!=0 and Z1!=0, 'Error: Z1 and Z2 must be specified'
        if len(data.columns)==4:
            assert Z1!=0 and Z2==0, 'Error: Only Z1 can be specified'
        if len(data.columns)==3:
            assert Z1==0 and Z2==0, 'Error: No Z can be specified'

        ## set column indexes
        if X==1:
            Xind = 1
        else:
            Xind = np.where(data.columns == X)[0][0]
        if Y==0:
            Yind =0
        else:
            Yind = np.where(data.columns==Y)[0][0]

        if Z1 !=0:
            Z1ind = np.where(data.columns==Z1)[0][0]
            z1 = np.unique(data.iloc[:, Z1ind])
        else:
            Z1ind = []
            z1 = np.zeros(1)
        if Z2 !=0:
            Z2ind = np.where(data.columns==Z2)[0][0]
            z2 = np.unique(data.iloc[:, Z2ind])
        else:
            Z2ind = []
            z2 = np.zeros(1)

        ## if self.color is empy, use default colormap
        if not self.colors:
            self.colors = cm.rainbow(np.linspace(0, 1, len(z1)))
        assert len(z2)<=15, 'Error: More values in Z2 than symbols'
        ## plot
        fig = plt.figure(figsize=[10,10])
        ax = plt.subplot(1,1,1)
        for i in range(len(z1)):
            for j in range(len(z2)):
                if Z1!=0 and Z2!=0:
                    ind = np.where(np.logical_and(data.iloc[:,Z1ind]==z1[i], data.iloc[:,Z2ind]==z2[j]))[0]
                elif Z1!=0 and Z2==0:
                    ind = np.where(data.iloc[:, Z1ind] == z1[i])[0]
                elif Z1==0 and Z2==0:
                    ind = np.where(data.iloc[:, 0] == data.iloc[:, 0])[0]
                ax.plot(data.iloc[ind,Xind],data.iloc[ind,Yind],linewidth=2,marker=self.linestyle[j],color=self.colors[i])

        ## create labels for legend
        if Z1 != 0:
            lines = ax.get_lines()
            labels_z1 = []
            ind_z1 = []
            for i,item in enumerate(z1):
                labels_z1.append('{}={}'.format(data.columns[Z1ind],item))
                ind_z1.append(i * len(z2))
            legend1 = plt.legend([lines[i] for i in ind_z1], labels_z1, loc=1,
                                 shadow=True, fancybox=True, edgecolor='inherit', framealpha=0.8)
            ax.add_artist(legend1)
        if Z2 !=0:
            labels_z2 = []
            ind_z2 = []
            for i,item in enumerate(z2):
                labels_z2.append('{}={}'.format(data.columns[Z2ind],item))
                ind_z2.append(i)
            legend2 = plt.legend([lines[i] for i in ind_z2], labels_z2, loc=4,
                             shadow=True, fancybox=True, edgecolor='inherit', framealpha=0.8)
            ax.add_artist(legend2)

        ## set x lim. check if x is a string (e.g. loc)
        if type(min(data.iloc[ind,Xind]))==str:
            ax.set_xlim([0, (len(data.iloc[ind, Xind]) -1)*1.5 ])
            plt.xticks(rotation=30)
            bottom = 0.2
        else:
            ax.set_xlim([min(data.iloc[ind, Xind]), max(data.iloc[ind, Xind]) * 1.5])
            bottom = 0.1

        ## figure style
        ax.grid('on',linestyle=':')
        plt.tick_params(labelsize=15)
        plt.subplots_adjust(top=0.8,bottom=bottom)
        plt.xlabel(data.columns[Xind],fontsize=20)
        plt.ylabel(data.columns[Yind],fontsize=20)
        plt.title(constant, fontsize=18)
        ax.grid('on',linestyle=':')
        plt.tick_params(labelsize=15)
        plt.savefig(path +'.png', dpi=500)

        fig.set_size_inches(11.69, 8.27)
        plt.savefig(path + '.pdf', papertype='a4')

#    def DBinput_old(self,X,Y,loc):
#        ## 1 check location
#        loc = loc
#
#        xlist = []
#        ylist = []
#        labels = []
#        while True:
#            ## 2 input waardes
#            Xinput = {}
#            for item in self.inputValues.keys():
#                if item !=X:
#                    tmp = float(input('Input {}='.format(item)))
#                    while tmp not in self.inputValues[item]:
#                        tmp = float(input('Bestaat niet! Input {}='.format(item)))
#                    Xinput[item] = tmp
#            print(Xinput)
#            ## 3 wind
#            windDir = float(input('WindDir ='))
#
#
#            x,y = self.create_sql_command(X,Y,loc,windDir,Xinput,fig=False)
#
#
#            tmp = input('Volgende reeks (y/n)=')
#            ## label
#            string = 'WindDir={}'.format(windDir)
#            for item in Xinput.keys():
#                string = string + '| {}={}'.format(item,Xinput[item])
#
#            if tmp=='y':
#                xlist.append(x)
#                ylist.append(y)
#                labels.append(string)
#            else:
#                xlist.append(x)
#                ylist.append(y)
#
#                labels.append(string)
#                break
#
#
#        plt.figure()
#        for i, item in enumerate(xlist):
#            plt.plot(xlist[i],ylist[i],label=labels[i])
#        plt.title(loc)
#        plt.legend()
#        plt.xlabel(X)
#        plt.ylabel(Y)
#        ##
#        plt.grid('on')

    def create_sql_command(self, path,Y, loc, windDir,closing, Xinput, fig=False):
        '''
        Make SQL command to extract data
        :param str path: location of databse
        :param str Y: result variabele
        :param str loc: location
        :param float windDir: wind direction
        :param int closing: closingcreterium
        :param dict Xinput: dictionary with values for all input variabele
        :param logical fig: switch beteen plotting results (only for testing)
        :return: Pandas.Series with the result variabele given the input variabele
        '''

        ## 1. Location
        command = "SELECT HRDLocationId FROM HRDLocations WHERE NAME='{}'".format(loc)
        HRDLocationId = self.read_sql_table(path ,command).values
        HRDLocationId = HRDLocationId[0,0]

        ## 2. wind direction
        windDir = str(windDir)
        command = "SELECT HRDWindDirectionId FROM HRDWindDirections WHERE Direction='{}'".format(windDir)
        HRDWindDirectionId = self.read_sql_table(path ,command).values
        HRDWindDirectionId = HRDWindDirectionId[0,0]

        ## 3. Closing
        ClosingsituationId = closing

        ## 4. Input Variable
        #command = "SELECT HRDInputColumnId FROM HRDInputVariables WHERE ColumnName='{}'".format(X)
        #HRDInputColumnId = self.read_sql_table(path ,command).values
        #HRDInputColumnId = HRDInputColumnId[0,0]

        ## 5. Result variable
        command = "SELECT HRDResultColumnId FROM HRDResultVariables WHERE ColumnName='{}'".format(Y)
        HRDResultColumnId = self.read_sql_table(path ,command).values
        HRDResultColumnId = HRDResultColumnId[0,0]

        ## get columnId from ColumnName
        Xvalue = {}
        for item in Xinput.keys():
            ColumnId = self.inputVar[item]
            Xvalue[ColumnId] = Xinput[item]
        N = len(Xvalue.keys())

        ## SQL command
        ## Select columns to return
        command1 = 'SELECT HydroDynamicResultData.Value '
        command2 = 'FROM HydroDynamicResultData '
        ## Join hydrodynamic data
        command3 = 'INNER JOIN HydroDynamicData '
        command4 = 'ON HydroDynamicData.HydroDynamicDataID=HydroDynamicResultData.HydroDynamicDataID '

        command_total = command1 + command2 + command3 + command4

        ## join inputData
        for i in range(N):
            command5 = 'INNER JOIN HydroDynamicInputData I{} '.format(i)
            command6 = 'ON I{}.HydroDynamicDataID=HydroDynamicData.HydroDynamicDataID '.format(i)
            command_total = command_total + command5 + command6

        ##  select result variable
        command7 = 'WHERE HydroDynamicResultData.HRDResultColumnId={} '.format(HRDResultColumnId)
        command8 = 'AND '
        command9 = 'HydroDynamicData.HRDLocationID={} AND  HydroDynamicData.ClosingSituationID={} AND HydroDynamicData.HRDWindDirectionID={} '.format(HRDLocationId,ClosingsituationId,HRDWindDirectionId)
        ## Input variable

        command_total = command_total + command7 + command8 + command9
        ## set contrains
        for i,item in enumerate(Xvalue.keys()):
            command11 = 'AND I{}.HRDInputColumnID={} '.format(i,item)
            command12 = 'AND '
            command13 = 'I{}.Value={} '.format(i,Xvalue[item])
            command_total = command_total + command11 + command12 + command13 #+ 'AND '

        value = self.read_sql_table(path ,command_total)

        ## plot for testing
        if fig==True:
            plt.figure()
            plt.plot(value.iloc[:,1],value.iloc[:,0],'.')
            plt.xlabel(X)
            plt.ylabel(Y)
            plt.title('{} and windDir={}'.format(loc,windDir))
        return value.iloc[:,0]


text:
#!/usr/bin/env python3
#
# Copyright (c) 2021, NVIDIA CORPORATION. All rights reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#

"""
This script demonstrates how to use the Calibrator API provided by Polygraphy
to calibrate a TensorRT engine to run in INT8 precision.
"""
import os

import numpy as np
from polygraphy.backend.trt import (Calibrator, CreateConfig,
                                    EngineFromNetwork, NetworkFromOnnxPath,
                                    TrtRunner)
from polygraphy.logger import G_LOGGER

MODEL = os.path.join(os.path.dirname(__file__), os.path.pardir, os.path.pardir, "models", "identity.onnx")
INPUT_SHAPE = (1, 1, 2, 2)

# The data loader argument to Calibrator can be any iterable or generator that yields `feed_dict`s.
# A feed_dict is just a mapping of input names to corresponding inputs (as NumPy arrays).
# Calibration will continue until our data loader runs out of data (4 batches in this example).
def calib_data():
    for _ in range(4):
        yield {"x": np.ones(shape=INPUT_SHAPE, dtype=np.float32)} # Totally real data

# We can provide a path or file-like object if we want to cache calibration data.
# This lets us avoid running calibration the next time we build the engine.
calibrator = Calibrator(data_loader=calib_data(), cache="identity-calib.cache")
build_engine = EngineFromNetwork(NetworkFromOnnxPath(MODEL), config=CreateConfig(int8=True, calibrator=calibrator))

# When we activate our runner, it will calibrate and build the engine. If we want to
# see the logging output from TensorRT, we can temporarily increase logging verbosity:
with G_LOGGER.verbosity(G_LOGGER.VERBOSE):
    with TrtRunner(build_engine) as runner:
        feed_dict = {"x": np.ones(shape=INPUT_SHAPE, dtype=np.float32)}
        outputs = runner.infer(feed_dict=feed_dict)
        assert np.all(outputs["y"] == feed_dict["x"])

ext: py

max_stars_repo_name: martellz/TensorRT

max_issues_repo_name: martellz/TensorRT

max_forks_repo_name: martellz/TensorRT

reason: import numpy


text:
#!/usr/bin/python3

import torch
from torch.autograd import Variable
from collections import defaultdict, OrderedDict
import numpy as np

def LoadWordVectors(wordVectorFile, dimension):
    word2vector = {}
    with open(wordVectorFile, 'r') as wordVectorsStream:
        for line in wordVectorsStream:
            items = line.strip().split(" ")
            if len(items) < dimension:
                continue
            word2vector[items[0]] = Variable(torch.FloatTensor(np.array([float(elem) for elem in items[1:]], dtype=np.float).reshape(dimension, 1)), requires_grad=True)
    return word2vector

def LoadWord2id(wordIdFile):
    word2id = defaultdict(lambda:0)
    with open(wordIdFile, 'r') as f:
        for line in f:
            items = line.strip().split(" ")
            if len(items) < 2:
                continue
            word2id[items[0]] = int(items[1])
    return word2id

def LoadEntityVectors(entityVectorFile, dimension):
    entity2vector = OrderedDict()
    entity2vector['reserved'] = Variable(torch.FloatTensor(np.array(np.random.uniform(-0.0, 0.0, (dimension, 1)), dtype=np.float)), requires_grad=True)
    id = 0
    with open(entityVectorFile, 'r') as entityVectorsStream:
        for line in entityVectorsStream:
            items = line.strip().split("\t")
            if len(items) < dimension:
                continue
            entity2vector[str(id)] = Variable(torch.FloatTensor(np.array([float(elem) for elem in items], dtype=np.float).reshape(dimension, 1)), requires_grad=True)
            id += 1
    return entity2vector

def LoadSamples(corpusPath):
    samples = []
    with open(corpusPath, 'r') as corpusStream:
        for line in corpusStream:
            items = line.strip().split("\t")
            if len(items) < 2:
                continue
            samples.append(items)
    return samples

def LoadMulInstanceSamples(corpusPath):
    samples = []
    with open(corpusPath, 'r') as corpusStream:
        for line in corpusStream:
            sentenceCount = int(line.strip().split("\t")[0])
            sentences = []
            for i in range(sentenceCount):
                sentences.append(corpusStream.readline().strip().split("\t"))
            samples.append(sentences)
    return samples

def LoadEntity2Id(filename):
    mapping = defaultdict( lambda: str(len(mapping)) )
    with open(filename, 'r') as f:
        # skep the first line, which is a total number.
        f.readline()
        for line in f:
            items = line.strip().split("\t")
            mapping[items[0]] = items[1]
    return mapping

def LoadRelationVectors(filename, dimension, cuda=False, requires_grad=False):
    relation2vector = defaultdict(lambda: unrelated)
    if cuda:
        unrelated = Variable(torch.FloatTensor(np.array(np.random.uniform(-0.0, 0.0, (dimension, 1)), dtype=np.float)).cuda(),  requires_grad=requires_grad)
    else:
        unrelated = Variable(torch.FloatTensor(np.array(np.random.uniform(-0.0, 0.0, (dimension, 1)), dtype=np.float)),  requires_grad=requires_grad)

    id = 0
    with open(filename, 'r') as relationVectorsStream:
        for line in relationVectorsStream:
            items = line.strip().split("\t")
            if len(items) < dimension:
                continue
            if cuda:
                relation2vector[str(id)] = Variable(torch.FloatTensor(np.array([float(elem) for elem in items], dtype=np.float).reshape(dimension, 1)).cuda(), requires_grad=True)
            else:
                relation2vector[str(id)] = Variable(torch.FloatTensor(np.array([float(elem) for elem in items], dtype=np.float).reshape(dimension, 1)), requires_grad=True)
            id += 1
    relation2vector['unknow'] = unrelated
    return relation2vector

def LoadTriples(filename):
    triple = {}
    with open(filename, 'r') as f:
        f.readline()
        for line in f:
            items = line.strip().split("\t")
            triple[items[0] + "_" + items[1]] = items[2]
    return triple

def LoadEntityName2FormatID(filename):
    mapping = defaultdict(lambda key:key)
    with open(filename, 'r') as f:
        for line in f:
            items = line.strip().split("\t")
            mapping[items[1].lower()] = items[0]
    return mapping

max_stars_repo_name: zhuango/KAN

max_issues_repo_name: Huiweizhou/MNM

max_forks_repo_name: Huiweizhou/MNM


text:
# -*- coding: utf-8 -*-
"""json_to_csv.ipynb

This script generates two csv files:
annotation_data.csv - containing bounding box and image category.
image_data.csv - containing width and height of image.

Required arguments
annotation json file and output directory to create csv's.
"""

import numpy as np
import pandas as pd
import json
import csv
import argparse
import os

def create_image_csv(data, output_file_dir):
  image_data = data['images']
  img_data_file = open(output_file_dir + '/image_data.csv', 'w')
  csv_writer = csv.writer(img_data_file)
  csv_writer.writerow(['img_id', 'filename', 'width', 'height'])

  for img_data in image_data:
      width = img_data['width']
      height = img_data['height']
      img_id = img_data['id']
      file_name = img_data['file_name']
      csv_writer.writerow([img_id, file_name, width, height])
  img_data_file.close()

def create_annotation_csv(data, output_file_dir):
  ann_data_file = open(output_file_dir + '/annotation_data.csv', 'w')
  csv_writer = csv.writer(ann_data_file)
  csv_writer.writerow(['img_id', 'img_label', 'xmin', 'ymin', 'xmax', 'ymax'])
  annotations_data = data["annotations"]

  for ann_data in annotations_data:
      img_id = ann_data['image_id']
      img_label = ann_data['category_id']
      xmin, ymin, bb_width, bb_height = ann_data['bbox']
    # xmax will be xmin + boundingbox width
    # ymax will be ymin + boundingbox height
      xmax = xmin + bb_width
      ymax = ymin + bb_height
      csv_writer.writerow([img_id, img_label, xmin, ymin, xmax, ymax])
  ann_data_file.close()

  create_image_csv(data, output_file_dir)

def convert_to_csv(json_file_path, output_file_dir):
  with open(json_file_path) as json_file:
    data = json.load(json_file)

  create_annotation_csv(data, output_file_dir)

def main():
  parser = argparse.ArgumentParser(description="Tensorflow JSON to CSV converter")
  parser.add_argument("-i", "--jsonFilePath",
                      help="Path to folder where json file is stored, including file name with extension",
                      type=str)
  parser.add_argument("-o", "--outputFileDir",
                      help="Output folder path to export csv files",
                      type=str)
  args = parser.parse_args()

  if(args.jsonFilePath is None):
    args.jsonFilePath = os.getcwd()
  if(args.outputFileDir is None):
    args.outputFileDir = args.jsonFilePath

  convert_to_csv(args.jsonFilePath, args.outputFileDir)
  print("JSON to CSV conversion successful")

if __name__ == '__main__':
  main()


text:
#!/usr/bin/python

__author__ = "Matan Lachmish"
__copyright__ = "Copyright 2016, Tel Aviv University"
__version__ = "1.0"
__status__ = "Development"

import numpy as np
import pickle
import sys
import os

import librosa

SOUND_SAMPLE_LENGTH = 30000

HAMMING_SIZE = 100
HAMMING_STRIDE = 40

def die_with_usage():
    """ HELP MENU """
    print 'USAGE: python preproccess.py [path to MSD mp3 data]'
    sys.exit(0)

def update_progress(progress):
    print '\r[{0}] {1}%'.format('#' * (progress / 10), progress)

def rreplace(s, old, new, occurrence):
    li = s.rsplit(old, occurrence)
    return new.join(li)

def prepossessingAudio(audioPath, ppFilePath):
    print 'Prepossessing ' + audioPath

    featuresArray = []
    for i in range(0, SOUND_SAMPLE_LENGTH, HAMMING_STRIDE):
        if i + HAMMING_SIZE <= SOUND_SAMPLE_LENGTH - 1:
            y, sr = librosa.load(audioPath, offset=i / 1000.0, duration=HAMMING_SIZE / 1000.0)

            # Let's make and display a mel-scaled power (energy-squared) spectrogram
            S = librosa.feature.melspectrogram(y, sr=sr, n_mels=128)

            # Convert to log scale (dB). We'll use the peak power as reference.
            log_S = librosa.logamplitude(S, ref_power=np.max)

            mfcc = librosa.feature.mfcc(S=log_S, sr=sr, n_mfcc=13)
            featuresArray.append(mfcc)

            # featuresArray.append(S)

            if len(featuresArray) == 599:
                break

    print 'storing pp file: ' + ppFilePath

    f = open(ppFilePath, 'w')
    f.write(pickle.dumps(featuresArray))
    f.close()

if __name__ == "__main__":

    # help menu
    if len(sys.argv) < 2:
        die_with_usage()

    i = 0.0
    walk_dir = sys.argv[1]

    print('walk_dir = ' + walk_dir)

    for root, subdirs, files in os.walk(walk_dir):
        for filename in files:
            if filename.endswith('.au'):
                file_path = os.path.join(root, filename)
                # print('\t- file %s (full path: %s)' % (filename, file_path))
                ppFileName = rreplace(file_path, ".au", ".pp", 1)

                # if os.path.isfile(ppFileName):  # Skip if pp file already exist
                #     continue

                try:
                    prepossessingAudio(file_path, ppFileName)
                except Exception as e:
                    print "Error accured" + str(e)

            if filename.endswith('au'):
                sys.stdout.write("\r%d%%" % int(i / 7620 * 100))
                sys.stdout.flush()
                i += 1

max_issues_repo_name: shradhit/MusicGenreClassification

max_forks_repo_name: shradhit/MusicGenreClassification


text:
#!/usr/bin/env python3

"""
Name:   utils.py
Author: Joshua Tice
Date:   November 18, 2018

Description
-----------
This file contains classes and fuctions that support both train.py and
predict.py, including:

Classifier          Class used for instantiating the fully-connected classifier
                    portion of the neural network model
build_model         Function that handles the contruction of a neural network
load_checkpoint     Function that contructs a neural network from a checkpoint
load_data           Function that creates an interface for training, validation,
                    and testing data
preprocess_image    Function for processing an image such that it can be
                    classified by the neural network

Required packages
-----------------
numpy
skimage
torch
torchvision
utils.py
"""

import json
import numpy as np
import pprint
from PIL import Image
import skimage
import skimage.io
import skimage.transform
import skimage.util
import torch
from torch import nn
from torch import optim
import torch.nn.functional as F
from torchvision import datasets, models, transforms

class Classifier(nn.Module):
    """Neural network classifier contructor"""

    def __init__(self, input_size, output_size, hidden_layers, drop_p=0.5):
        """
        Parameters
        ----------
        input_size : int
            Number of input nodes
        output_size : int
            Number of output nodes
        hidden_layers : list of ints
            Number of nodes in each hidden layer of the classifier
        drop_p : float, optional
            Dropout of the classifier (between 0.0 and 1.0)
        """

        super().__init__()
        self.hidden_layers = nn.ModuleList([nn.Linear(input_size, hidden_layers[0])])
        layer_sizes = zip(hidden_layers[:-1], hidden_layers[1:])
        self.hidden_layers.extend([nn.Linear(h1, h2) for h1, h2 in layer_sizes])
        self.output = nn.Linear(hidden_layers[-1], output_size)
        self.dropout = nn.Dropout(p=drop_p)

    def forward(self, x):
        """
        Operations defining the forward pass of the neural network

        Parameters
        ----------
        x : torch tensor
            Tensor of images

        Returns
        -------
        array-like
            vector of predictions
        """

        for linear in self.hidden_layers:
            x = F.relu(linear(x))
            x = self.dropout(x)

        x = self.output(x)

        return F.log_softmax(x, dim=1)

def build_model(base_model='vgg16', input_size=25088, output_size=102,
                hidden_layers=[4096, 1000], drop_p=0.5):
    """
    Loads a pre-trained convolutional neural network and configures a fully-
    connected classifier

    Parameters
    ----------
    base_model : str, optional
        The convolutional neural network the model will be based on.
        {'vgg16', 'vgg13'}
    input_size : int, optional
        Number of input nodes
    output_size : int, optional
        Number of output nodes
    hidden_layers : list of ints, optional
        Number of nodes in each hidden layer of the classifier
    drop_p : float, optional
        Dropout of the classifier (between 0.0 and 1.0)

    Returns
    -------
    torch convolutional neural network
    """

    if base_model == 'vgg11':
        model = models.vgg11(pretrained=True)
    elif base_model == 'vgg13':
        model = models.vgg13(pretrained=True)
    elif base_model == 'vgg16':
        model = models.vgg16(pretrained=True)
    else:
        raise ValueError(
            "base_model can only accept values in {'vgg11', 'vgg13', 'vgg16'}")

    for param in model.parameters():
        param.requires_grad = False

    model.classifier = Classifier(input_size, output_size, hidden_layers, drop_p)

    return model

def load_checkpoint(filepath, device='cpu'):
    """
    Instantiates and configures a torch neural network based on the parameters
    of a given checkpoint.

    Parameters
    ----------
    filepath : str
        File path to a saved checkpoint
    device : str
        Device on which to load the model {'cpu', 'cuda:0'}

    Returns
    -------
    torch neural network
        Instantiated network configured with checkpoint parameters
    """

    checkpoint = torch.load(filepath, map_location=device)
    model = build_model(
        base_model=checkpoint['model'],
        input_size=checkpoint['classifier_input_size'],
        output_size=checkpoint['classifier_output_size'],
        hidden_layers=checkpoint['classifier_hidden_layers'],
        drop_p=checkpoint['classifier_dropout'])
    model.load_state_dict(checkpoint['model_state_dict'])
    model.class_to_index = checkpoint['class_to_index']
    model.to(torch.device(device))

    return model

def load_data(data_dir):
    """
    Loads data from data directory for neural network training

    Parameters
    ----------
    data_dir : str
        Path to the data directory containing images for training

    Returns
    -------
        dict
            Dictionary with training, validation, and testing data loaders in
            addition to dictionary with class to index relationship
    """
    training_dir = data_dir + '/train'
    validation_dir = data_dir + '/valid'
    testing_dir = data_dir + '/test'

    means = (0.485, 0.456, 0.406)
    std_devs = (0.229, 0.224, 0.225)

    training_transforms = transforms.Compose([
        transforms.RandomHorizontalFlip(),
        transforms.RandomRotation(degrees=30, expand=True),
        transforms.RandomResizedCrop(size=224),
        transforms.ToTensor(),
        transforms.Normalize(means, std_devs),
    ])
    validation_transforms = transforms.Compose([
        transforms.Resize(224),
        transforms.CenterCrop(size=(224, 224)),
        transforms.ToTensor(),
        transforms.Normalize(means, std_devs),
    ])
    testing_transforms = transforms.Compose([
        transforms.Resize(224),
        transforms.CenterCrop(size=(224, 224)),
        transforms.ToTensor(),
        transforms.Normalize(means, std_devs),
    ])

    training_data = datasets.ImageFolder(
        training_dir, transform=training_transforms)
    validation_data = datasets.ImageFolder(
        validation_dir, transform=validation_transforms)
    testing_data = datasets.ImageFolder(
        testing_dir, transform=testing_transforms)

    training_loader = torch.utils.data.DataLoader(
        training_data, batch_size=64, shuffle=True)
    validation_loader = torch.utils.data.DataLoader(
        validation_data, batch_size=32)
    testing_loader = torch.utils.data.DataLoader(
        testing_data, batch_size=32)

    data = {
        'training_loader': training_loader,
        'validation_loader': validation_loader,
        'testing_loader': testing_loader,
        'class_to_index': training_data.class_to_idx
    }

    return data

def preprocess_image(image_path):
    '''
    Scales, crops, and normalizes a PIL image for a PyTorch model

    Parameters
    ----------
    image_path : str
        Path to the image to be pre-processed

    Returns
    -------
    numpy.array
        Processed numpy image
    '''

    image = Image.open(image_path)

    # Scale image
    scaling_factor = 256 / min(image.size)
    scaled_size = (int(image.size[0] * scaling_factor),
                   int(image.size[1] * scaling_factor))
    image = image.resize(scaled_size, Image.ANTIALIAS)

    # Crop image
    box = ((int(image.size[1] / 2 - 112),
            int(image.size[0] / 2 - 112),
            int(image.size[1] / 2 + 112),
            int(image.size[0] / 2 + 112),))
    image = image.crop(box)

    # Normalize image
    np_image = np.array(image, dtype=float)
    means = np.array([0.485, 0.456, 0.406])
    std_devs = np.array([0.229, 0.224, 0.225])
    np_image /= 255
    np_image = (np_image - means) / std_devs

    # Transpose image for pytorch model
    np_image = np_image.transpose((2, 0, 1))

    return np_image


text:
# --------------------------------------------------------
# Tensorflow VCL
# Licensed under The MIT License [see LICENSE for details]
# Written by Zhi Hou, based on code from Chen Gao, Zheqi he and Xinlei Chen
# --------------------------------------------------------

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
import os
import pickle
import _init_paths
import numpy as np
from torch.utils.data.dataloader import DataLoader
from torch.utils.data.dataset import Dataset
from torchvision.transforms import transforms, ToTensor

from networks.ResNet50_HICO_torch import HICO_HOI
from ult.timer import Timer

os.environ['DATASET'] = 'HICO'

import argparse

from ult.config import cfg
from models.test_HICO import obtain_test_dataset1, test_net_data_api1

def parse_args():
    parser = argparse.ArgumentParser(description='Test VCL on HICO')
    parser.add_argument('--num_iteration', dest='iteration',
            help='Specify which weight to load',
            default=10, type=int)
    parser.add_argument('--model', dest='model',
            help='Select model',
            default='VCL_union_l2_rew_aug5_3_x5new_res101', type=str)
    parser.add_argument('--object_thres', dest='object_thres',
            help='Object threshold',
            default=0.3, type=float)
    parser.add_argument('--human_thres', dest='human_thres',
                        help='Human threshold',
                        default=0.8, type=float)
    # TODO For better object detector, the object_thres and human_thres should also be changed accordingly.
    #  e.g. in our fine-tuned detector, object_thres and human_thres is 0.1 and 0.3 respectively.
    parser.add_argument('--debug', dest='debug',
                        help='Human threshold',
                        default=0, type=int)
    parser.add_argument('--type', dest='test_type',
                        help='Human threshold',
                        default='vcl', type=str)
    args = parser.parse_args()

    return args

class HicoDataset(Dataset):

    def __init__(self, Pos_augment=15, Neg_select=60, augment_type=0, with_pose=False, zero_shot_type=0,
                 large_neg_for_ho=False, isalign=False, epoch=0, transform=None):
        self.transform = transform

        self.generator = obtain_test_dataset1(args.object_thres, args.human_thres,
                         stride=stride, test_type=args.test_type, model_name=args.model)

    def __len__(self):
        return 800000

    def __getitem__(self, idx):
        # im_orig, image_id, num_pos, Human_augmented, Object_augmented, action_HO, Pattern = next(self.generator)
        im_orig, blobs, image_id  = next(self.generator)
        # im_orig = im_orig.transpose([0, 3, 1, 2])
        # Pattern = Pattern.transpose([0, 3, 1, 2]).astype(np.float32)
        # Human_augmented = Human_augmented.astype(np.float32)
        # Object_augmented = Object_augmented.astype(np.float32)
        # Human_augmented = Human_augmented.astype(np.float32)
        # print(im_orig.dtype, Pattern.dtype)
        # print(im_orig)
        # print(im_orig.shape, im_orig)
        blobs['sp'] = blobs['sp'].transpose([0, 3, 1, 2]).astype(np.float32)
        if self.transform:
            im_orig = self.transform(im_orig)
        # print('after', im_orig)
        print(np.asarray(im_orig).shape)
        return im_orig, image_id, blobs

if __name__ == '__main__':

    args = parse_args()
    print(args)
    # test detections result
    from sys import version_info

    weight = cfg.ROOT_DIR + '/Weights/' + args.model

    import os

    print('weight:', weight)
    print ('Human thres = ' + str(args.human_thres) + ', Object thres = ' + str(args.object_thres) + ', iter = ' + str(args.iteration) + ', path = ' + weight )
    output_file = cfg.LOCAL_DATA + '/Results/' + str(args.iteration) + '_' + args.model + '_tin.pkl'

    HICO_dir = cfg.ROOT_DIR + '/Results/HICO/' + str(args.iteration) + '_' + args.model + '/'
    stride = 200
    # test_generator = obtain_test_dataset1(args.object_thres, args.human_thres,
    #                      stride=stride, test_type=args.test_type, model_name=args.model)
    import torch
    model = HICO_HOI(args.model)
    checkpoint = torch.load(weight + '/' + '{}_checkpoint.pth.tar'.format(args.iteration))
    model.load_state_dict(checkpoint['state_dict'], strict=True)
    for name, p in model.named_parameters():
        print(name, p.mean())

    dataset = HicoDataset(
                          transform=transforms.Compose([ToTensor(), transforms.Normalize(mean=[0.485, 0.456, 0.406],
                                                                                std=[0.229, 0.224, 0.225])]),
                          )
    dataloader_test = DataLoader(dataset, 1,
                                  shuffle=False, num_workers=1)
    # transform = transforms.Compose([ToTensor(), transforms.Normalize(mean=[0.485, 0.456, 0.406],
    #                                                                  std=[0.229, 0.224, 0.225])])
    device = torch.device("cuda")
    model.to(device)
    model.eval()
    detection = {}
    count = 0
    last_update_value = {}
    _t = {'im_detect': Timer(), 'misc': Timer()}
    for image, image_id, blobs in dataloader_test:
        _t['im_detect'].tic()
        print(image.shape, image.mean(), image.std())
        print(image.shape, image.mean(), image.std())

        model(image.to(device), image_id[0], blobs['H_num'][0],
              blobs['H_boxes'][0].to(device),
              blobs['O_boxes'][0].to(device), None,
              blobs['sp'][0].to(device), False)
        sp = model.predictions["cls_prob_sp"]
        hoi = model.predictions["cls_prob_hoi"]
        spHOI = torch.mul(sp, hoi)
        import ipdb
        ipdb.set_trace()
        #
        sp = sp.detach().cpu().numpy()
        hoi = hoi.detach().cpu().numpy()

        # print([h[:10] for h in hoi])
        # print([h[:10] for h in sp])
        # print(image.mean(), image.std())

        spHOI = spHOI.detach().cpu().numpy()

        if count == 10:
            for k in model.state_dict().keys():
                tmp = model.state_dict()[k].type(torch.float32).mean().detach().cpu().numpy()
                if k in last_update_value:
                    print(tmp)
                    if abs(last_update_value[k] - tmp) > 0:
                        print(k, last_update_value[k], tmp, last_update_value[k] - tmp)

                last_update_value[k] = tmp
                # print(k, model.state_dict()[k].type(torch.float32).mean())
            print('-'*80)

        last_img_id = -1
        # ipdb.set_trace()
        # print(len(blobs['H_boxes']), )
        temp = [[blobs['H_boxes'][i][1:], blobs['O_boxes'][i][1:], blobs['O_cls'][i],
                 0, blobs['H_score'][i], blobs['O_score'][i], 0, 0, sp[i], hoi[i],
                 0] for i in range(len(blobs['H_boxes']))]
        if image_id in detection:
            detection[image_id].extend(temp)
        else:
            detection[image_id] = temp

        _t['im_detect'].toc()
        count += 1

        print('\rmodel: {} im_detect: {:d}/{:d}  {:d}, {:.3f}s'.format(model.model_name, count, 10566, image_id,
                                                                       _t['im_detect'].average_time), end='',
              flush=True)
        # if count > 100:
        #     break

    # TODO remove
    # pickle.dump(detection, open('test_new.pkl', 'wb'))
    pickle.dump(detection, open(output_file, "wb"))


text:
#
# Author: Cunren Liang
# Copyright 2015-present, NASA-JPL/Caltech
#

import os
import logging
import numpy as np

import isceobj
from isceobj.Alos2Proc.runGeocode import geocode
from isceobj.Alos2Proc.Alos2ProcPublic import getBboxGeo

logger = logging.getLogger('isce.alos2insar.runGeocodeSd')

def runGeocodeSd(self):
    '''geocode final products
    '''
    catalog = isceobj.Catalog.createCatalog(self._insar.procDoc.name)
    self.updateParamemetersFromUser()

    masterTrack = self._insar.loadTrack(master=True)
    #slaveTrack = self._insar.loadTrack(master=False)

    demFile = os.path.abspath(self._insar.demGeo)

    sdDir = 'sd'
    if not os.path.exists(sdDir):
        os.makedirs(sdDir)
    os.chdir(sdDir)

    if self.geocodeListSd == None:
        geocodeList = self._insar.multilookCoherenceSd + self._insar.azimuthDeformationSd + self._insar.maskedAzimuthDeformationSd
    else:
        geocodeList = self.geocodeListSd

    if self.bbox == None:
        bbox = getBboxGeo(masterTrack)
    else:
        bbox = self.bbox
    catalog.addItem('geocode bounding box', bbox, 'runGeocodeSd')

    numberRangeLooks = self._insar.numberRangeLooks1 * self._insar.numberRangeLooksSd
    numberAzimuthLooks = self._insar.numberAzimuthLooks1 * self._insar.numberAzimuthLooksSd

    for inputFile in geocodeList:
        if self.geocodeInterpMethodSd == None:
            img = isceobj.createImage()
            img.load(inputFile + '.xml')
            if img.dataType.upper() == 'CFLOAT':
                interpMethod = 'sinc'
            else:
                interpMethod = 'bilinear'
        else:
            interpMethod = self.geocodeInterpMethodSd.lower()

        geocode(masterTrack, demFile, inputFile, bbox, numberRangeLooks, numberAzimuthLooks, interpMethod, 0, 0)

    os.chdir('../')

    catalog.printToLog(logger, "runGeocodeSd")
    self._insar.procDoc.addAllFromCatalog(catalog)


text:
from __future__ import print_function
from __future__ import unicode_literals

import scipy.stats
import sklearn.ensemble

"""
Helpers for feature selection. Not entirely generalized from Mars Express Challenge code.
"""

import numpy
import sklearn.feature_selection
import sklearn.linear_model

def score_features_elastic_net(X, Y, alpha=0.001):
    """ElasticNet coefs against summed Y to simplify multivariate"""
    model = sklearn.linear_model.ElasticNet(alpha)
    model.fit(X, Y.sum(axis=1))
    return sklearn.feature_selection.from_model._get_feature_importances(model)

def score_features_ridge(X, Y, alpha=0.01):
    model = sklearn.linear_model.Ridge(alpha)
    model.fit(X, Y.sum(axis=1))
    return sklearn.feature_selection.from_model._get_feature_importances(model)

class BaseFeatureSelectorCV(object):
    """Base class for cross-validated feature selection. It's here because sklearn doesn't have any such functionality"""
    def __init__(self, std_weight=-0.05):
        self.std_weight = std_weight

    def _score(self, X_train, Y_train, X_test, Y_test):
        pass

    def score(self, X, Y, splits):
        score_matrix = []
        for train, test in splits:
            score_matrix.append(self._score(X[train], Y[train], X[test], Y[test]))

        score_matrix = numpy.asarray(score_matrix)
        feature_scores = score_matrix.mean(axis=0) + self.std_weight * score_matrix.std(axis=0)
        return feature_scores

class LoiSelector(BaseFeatureSelectorCV):
    """Feature selection using leave-one-in, which measures the quality of the model with each feature independently."""
    def __init__(self, model, scorer, std_weight=-0.05):
        super(LoiSelector, self).__init__(std_weight=std_weight)
        self.model = model
        self.scorer = scorer

    def _score(self, X_train, Y_train, X_test, Y_test):
        scores = [0 for _ in range(X_train.shape[1])]
        for i in range(X_train.shape[1]):
            self.model.fit(X_train[:, [i]], Y_train)
            scores[i] = self.scorer(self.model, X_test[:, [i]], Y_test)

        return scores

class LooSelector(BaseFeatureSelectorCV):
    """Feature selection using leave-one-out, which measures the impact of removing each feature from the total set.
    In comparison to leave-one-in, if the model handles combinations of features this can assess each feature's contrib."""
    def __init__(self, model, scorer, std_weight=-0.05):
        super(LooSelector, self).__init__(std_weight=std_weight)
        self.model = model
        self.scorer = scorer

    def _score(self, X_train, Y_train, X_test, Y_test):
        scores = [0 for _ in range(X_train.shape[1])]

        self.model.fit(X_train, Y_train)
        baseline_score = self.scorer(self.model, X_test, Y_test)

        for i in range(X_train.shape[1]):
            included = numpy.asarray([j for j in range(X_train.shape[1]) if j != i])
            self.model.fit(X_train[:, included], Y_train)
            scores[i] = baseline_score - self.scorer(self.model, X_test[:, included], Y_test)

        return scores

class DeformSelector(BaseFeatureSelectorCV):
    def __init__(self, model, scorer, std_weight=-0.05):
        super(DeformSelector, self).__init__(std_weight=std_weight)
        self.model = model
        self.scorer = scorer

    @staticmethod
    def deform_feature(X, i, adjustment=0.1):
        X_def = X.copy()
        X_def[:, i] *= 1 + adjustment
        return X_def

    def _score(self, X_train, Y_train, X_test, Y_test):
        self.model.fit(X_train, Y_train)

        baseline_score = self.scorer(self.model, X_test, Y_test)
        deform_scores_increased = numpy.asarray([self.scorer(self.model, self.deform_feature(X_test, i, adjustment=0.1), Y_test) for i in range(X_train.shape[1])])
        deform_scores_decreased = numpy.asarray([self.scorer(self.model, self.deform_feature(X_test, i, adjustment=-0.1), Y_test) for i in range(X_train.shape[1])])

        return numpy.vstack((baseline_score - deform_scores_decreased, baseline_score - deform_scores_increased)).mean(axis=0)

def score_features_random_forest(X, Y, n_jobs=-1):
    """Random forest feature importances, supports multivariate"""
    model = sklearn.ensemble.RandomForestRegressor(n_estimators=100, n_jobs=n_jobs)
    model.fit(X, Y)
    return sklearn.feature_selection.from_model._get_feature_importances(model)

def inverse_rank_order(weights, base=0.9):
    """Inverse rank order, can be used to bring different feature scores onto the same scale but it didn't work well for me"""
    return base ** scipy.stats.rankdata(weights)

def multivariate_select(X, Y, univariate_feature_scoring_function, weight_outputs=False):
    """
    Handle multivariate feature selection by running univariate on each output and averaging.
    :param weight_outputs: Use the mean + stddev of outputs to weight output importance
    """
    output_weights = Y.mean(axis=0) + Y.std(axis=0)

    scores = []
    for output_index in range(Y.shape[1]):
        output = Y[:, output_index]
        scores.append(univariate_feature_scoring_function(X, output))

    score_matrix = numpy.vstack(scores) # M outputs x N features

    if weight_outputs:
        # should be 1 x N
        return output_weights.dot(score_matrix)
    else:
        return score_matrix.mean(axis=0)

max_stars_repo_name: ktrnka/helpers

max_issues_repo_name: ktrnka/helpers

max_forks_repo_name: ktrnka/helpers

reason: import numpy,import scipy


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

# # Pooling and Blending
#
# This notebook presents an overview of bilinear pooling and blending problems in the context of a simple milk blending operation. The essential non-convex nature of the problem is demonstrated, and the two basic formulations (the P- and Q- parameterizations) are shown.

# In[1]:

# Install Pyomo and solvers for Google Colab
import sys
if "google.colab" in sys.modules:
    get_ipython().system('wget -N -q https://raw.githubusercontent.com/jckantor/MO-book/main/tools/install_on_colab.py ')
    get_ipython().run_line_magic('run', 'install_on_colab.py')

# ## Bibliographic Notes
#
# The pooling and blending is a large scale, high value, fundamental problem of logistics for the process and refining industries. The prototypical examples are the pooling and blending crude oils to meet the feed stock constraints of refineries, and for the pooling of refinery products for pipeline delivery to distribution terminals. Un
#
# Haverly (1978) is a commonly cited small benchmark problem for the pooling and blending of sulfurous fuels.
#
# > Haverly, C. A. (1978). Studies of the behavior of recursion for the pooling problem. Acm sigmap bulletin, (25), 19-28. https://dl.acm.org/doi/pdf/10.1145/1111237.1111238
#
# There is an extensive literature on pooling and blending. The following encyclopedia entry explains the history of the pooling problem, how it leads to multiple local minima and other pathological behaviors, and approaches to finding practical solutions.
#
# > Visweswaran, V. (2009). MINLP: Applications in Blending and Pooling Problems. https://link.springer.com/referenceworkentry/10.1007/978-0-387-74759-0_375
#
# Recent research overviews include
#
#
# > Misener, R., & Floudas, C. A. (2009). Advances for the pooling problem: Modeling, global optimization, and computational studies. Applied and Computational Mathematics, 8(1), 3-22. https://www.researchgate.net/profile/Ruth-Misener/publication/242290955_Advances_for_the_pooling_problem_Modeling_global_optimization_and_computational_studies_Survey/links/0046352e7d1dfeb40f000000/Advances-for-the-pooling-problem-Modeling-global-optimization-and-computational-studies-Survey.pdf
#
# > Gupte, A., Ahmed, S., Dey, S. S., & Cheon, M. S. (2013). Pooling problems: relaxations and discretizations. School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, GA. and ExxonMobil Research and Engineering Company, Annandale, NJ. http://www.optimization-online.org/DB_FILE/2012/10/3658.pdf
#
# The current state-of-the-art appears to be a formulation of the pooling problem is a mixed-integer quadratically-constrained quadratic optimization on a given network.
#
# > Ceccon, F., & Misener, R. (2022). Solving the pooling problem at scale with extensible solver GALINI. Computers & Chemical Engineering, 107660. https://arxiv.org/pdf/2105.01687.pdf
#
# Applications for pooling and blending are probably underappreciated. In particular, what role might pooling and blending problems have in projects like the World Food Programme (WFP)?

# ## Example: Pooling Milk for Wholesale Blending and Distribution
#
# A bulk distributor supplies several customers with deliveries of raw milk purchased from local farms. The deliveries may be blended to meet the minimum milk fat content specified by each customer.
#
# The distributor has found lower cost sources of raw milk from several farms located a long distance away. The distant farms produce grades of milk that may be blended with local sources. The distributor, however, has just one truck with a single tank available to transport raw milk from the distant farms to the blending facility. Milk purchased from the distant farms would be mixed together in the tank before transport to the blending station. The contents of the tank form a "pool" of uniform composition that can be combined with local sources to meet the specific needs of each customer.
#
# The situation is illustrated in the following diagram where arrows show the flow and blending of milk supplies.
#
# ![](milk-pooling.dio.png)
#
# What should the distributor do?
#
# * Option 1. Do nothing, continue operating the business as usual with local suppliers.
#
# * Option 2. Buy a second truck so raw milk can be transported from the remote farms to the blending facility without pooling.
#
# * Option 3. Pool raw milk from the remote farms into a single tank for truck transport to the blending facility.
#

# In[2]:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

customers = pd.DataFrame({
    "Customer 1": {"fat": 0.045, "price": 52.0, "demand": 6000.0},
    "Customer 2": {"fat": 0.030, "price": 48.0, "demand": 2500.0},
    "Customer 3": {"fat": 0.040, "price": 50.0, "demand": 4000.0},
}).T

suppliers = pd.DataFrame({
    "Farm A": {"fat": 0.045, "cost": 45.0, "location": "local"},
    "Farm B": {"fat": 0.030, "cost": 42.0, "location": "local"},
    "Farm C": {"fat": 0.033, "cost": 37.0, "location": "remote"},
    "Farm D": {"fat": 0.050, "cost": 45.0, "location": "remote"}},
    ).T

local_suppliers = suppliers[suppliers["location"]=="local"]
remote_suppliers = suppliers[suppliers["location"]=="remote"]

print("\nCustomers")
display(customers)

print("\nLocal Suppliers")
display(local_suppliers)

print("\nRemote Suppliers")
display(remote_suppliers)

# ## Option 1. Business as usual
#
# Business as usual for the milk distributor is to blend supplies from the two local farms to meet customer requirements. Let $L$ designate the set of local suppliers, and let $C$ designate the set of customers. The decision variable $z_{l, c}$ is the amount of milk from local supplier $l\in L$ mixed into the blend for customer $c\in C$.
#
# The distributor's objectives is to maximize profit
#
# $$
# \begin{align*}
# \text{Profit} & = \sum_{(l, c)\ \in\ L \times C} (\text{price}_c - \text{cost}_l) z_{l,c}
# \end{align*}
# $$
#
# where the notation $(l, c)\ \in\ L\times C$ indicates a summation over the cross-product of two sets. A useful interpretation is that $(l, c)\ \in\ L\times C$ describes all ways of delivering milk from $l$ to $c$. Each term $(\text{price}_c - \text{cost}_l)$ is then the profit earned by delivering one unit of milk from $l\in L$ to $c\in C$.
#
# The amount of milk delivered to each customer $c\in C$ can not exceed the demand.
#
# $$
# \begin{align*}
# \sum_{l\in L} z_{l, c} & \leq \text{demand}_{c} & \forall c\in C
# \end{align*}
# $$
#
# The milk blend delivered to each customer $c\in C$ must meet the minimum product quality requirement for milk fat. Assuming linear blending, the model becomes
#
# $$
# \begin{align*}
# \sum_{(l,c)\ \in\ L \times C} \text{conc}_{l} z_{l,c} & \geq \text{conc}_{c} \sum_{l\in L} z_{l, c} & \forall c \in C
# \end{align*}
# $$
#
# This is a standard linear blending problem.

# In[3]:

import pyomo.environ as pyo

q = "fat"

m = pyo.ConcreteModel()

# define sources and customers
m.L = pyo.Set(initialize=local_suppliers.index)
m.C = pyo.Set(initialize=customers.index)

# set of quality components
m.Q = pyo.Set(initialize=["fat"])

# define local -> customer flowrates
m.z = pyo.Var(m.L * m.C, domain=pyo.NonNegativeReals)

@m.Objective(sense=pyo.maximize)
def profit(m):
    return sum(m.z[l, c]*(customers.loc[c, "price"] - suppliers.loc[l, "cost"]) for l, c in m.L * m.C)

@m.Constraint(m.C)
def demand(m, c):
    return sum(m.z[l, c] for l in m.L) <= customers.loc[c, "demand"]

@m.Constraint(m.C)
def quality(m, c):
    return sum(m.z[l, c] * suppliers.loc[l, q] for l in m.L) >= sum(m.z[l, c] for l in m.L) * customers.loc[c, q]

pyo.SolverFactory('cbc').solve(m)

# report results
print(f"\nprofit = {m.profit():0.2f}\n")
Z = pd.DataFrame([[l, c, round(m.z[l, c](), 1)] for l, c in m.L * m.C], columns = ["supplier", "customer", "blend"])
Z = Z.pivot_table(index="customer", columns="supplier")
Z["fat"]  = sum(Z.loc[:,('blend',l)] * suppliers.loc[l, "fat"] for l in m.L)/Z.sum(axis=1)
Z["Total"] = Z.sum(axis=1)
display(Z)

# ## Option 2. Buy an additional truck
#
# The current profit earned is 81,000 using only local suppliers. Before considering pooling, however, the distributor may wish to know the maximum profit that could be earned if the remote suppliers could be blended as if they were local suppliers. This would require acquiring and operating a separate transport for each remote supplier, but it is worth knowing if the additional profit earned would justify the additional expense.
#
# The model used above applies by extending the set of suppliers to include both local and remote farms. The required changes are noted in the cell below.

# In[4]:

import pyomo.environ as pyo

q = "fat"

m = pyo.ConcreteModel()

# define sources and customers
m.R = pyo.Set(initialize=remote_suppliers.index)       # <== set of remote suppliers
m.L = m.R | pyo.Set(initialize=local_suppliers.index)  # <== treat remote suppliers same a local suppliers
m.C = pyo.Set(initialize=customers.index)

# set of quality components
m.Q = pyo.Set(initialize=["fat"])

# define local -> customer flowrates
m.z = pyo.Var(m.L * m.C, domain=pyo.NonNegativeReals)

@m.Objective(sense=pyo.maximize)
def profit(m):
    return sum(m.z[l, c]*(customers.loc[c, "price"] - suppliers.loc[l, "cost"]) for l, c in m.L * m.C)

@m.Constraint(m.C)
def demand(m, c):
    return sum(m.z[l, c] for l in m.L) <= customers.loc[c, "demand"]

@m.Constraint(m.C)
def quality(m, c):
    return sum(m.z[l, c] * suppliers.loc[l, q] for l in m.L) >= sum(m.z[l, c] for l in m.L) * customers.loc[c, q]

pyo.SolverFactory('cbc').solve(m)

# report results
print(f"\nprofit = {m.profit():0.2f}\n")
Z = pd.DataFrame([[l, c, round(m.z[l, c](), 1)] for l, c in m.L * m.C], columns = ["supplier", "customer", "blend"])
Z = Z.pivot_table(index="customer", columns="supplier")
Z["fat"]  = sum(Z.loc[:,('blend',l)] * suppliers.loc[l, "fat"] for l in m.L)/Z.sum(axis=1)
Z["Total"] = Z.sum(axis=1)
display(Z)

# The calculation shows that including raw milk from the remote farms could significantly increases profits. This blending requires transporting raw milk from farms C and  D to the blending station. This calculation assumed the milk would be delivered at the original concentration, which will require at least two trucks to keep the sources of milk separated until they reach the blending facility. The plan does not allow pooling these sources prior to transport, otherwise they would need to be "unblended" to meet the differing quality requirements of the customers.

# ## Option 3. Pool delivery from remote farms
#
# Comparing Option 1 with Option 2 shows there is significantly more profit to be earned by purchasing raw milk from the remote farms. But that option requires an additional truck to keep the supplies separated during transport. Otherwise mixing raw milk from the remote farms would result in a uniform mixture that would not meet the requirements of both customers.
#
# Because only one truck with a single tank is available for transport from the remote farms, the new problem is to combine purchases from the remote farms into a single pool of uniform composition that could be blended milk from local farms to meet individual customer requirements. Compared to option 2, the profit potential may be reduced due to pooling, but at least it does not require an additional truck.
#
# This is the pooling problem. There are a several of mathematical formulations in the literature of this problem. The first analysis here uses a formulation called the "p-parameterization" where the pool composition is a new decision variable $p$. Other decision new variables are $x_r$ which are the amounts of raw milk purchased from remote farms $r\in R$, and $y_c$ which are the amounts delivered to customer $c\in C$ from the pool.
#
# The profit objective now includes the cost of purchasing raw milk from the remote farms and the income received for selling material from the pool.
#
# $$
# \begin{align*}
# \text{Profit} & = \sum_{(l,c)\ \in\ L \times C} (\text{price}_c - \text{cost}_l)\ z_{l,c}
# + \sum_{c\in C} \text{price}_c y_{c} - \sum_{r\in R} \text{cost}_r x_{r}
# \end{align*}
# $$
#
# The product delivered to each customer from local farms and the pool can not exceed demand.
#
# $$
# \begin{align*}
# \sum_{l\in L} z_{l, c} + y_{c} & \leq \text{demand}_{c} & \forall c\in C
# \end{align*}
# $$
#
# Purchases from the remote farms and the amounts delivered to customers from the pool must balance.
#
# $$
# \begin{align*}
# \sum_{r\in R}x_{r} & = \sum_{c\in C} y_{c} \\
# \end{align*}
# $$
#
# The average milk fat composition of the pool, $p$, must satisfy an overall balance on milk fat entering the pool from the remote farms and the milk fat delivered to customers.
#
# $$
# \begin{align*}
# \sum_{r\in R}\text{conc}_{r}\ x_{r}  & = \underbrace{p \sum_{c\in C} y_{c}}_{\text{bilinear}}
# \end{align*}
# $$
#
# Finally, the milk fat required by each customer $c\in C$ satisfies a blending constraint.
#
# $$
# \begin{align*}
# \underbrace{p y_{c}}_{\text{bilinear}}  + \sum_{(l,c)\ \in\ L \times C} \text{conc}_{l}\ z_{l,c}
# & \geq \text{conc}_{c}\ (\sum_{l\in L} z_{l, c} + y_{c})
# & \forall c \in C
# \end{align*}
# $$
#
# The last two constraints include bilinear terms from the project of decision variable $p$ with decision variables $y_c$ for all $c\in C$.
#
# The bilinear terms have a profound consequence on the nature of the optimization problem. To demonstrate, the following cell creates a linear program to maximize profit as a function of $p$, then explores how profit changes as a function of parameter $p$.

# In[5]:

import pyomo.environ as pyo

def milk_pooling(p=0, q="fat"):

    m = pyo.ConcreteModel('Milk Pooling Model')

    # define sources
    m.L = pyo.Set(initialize=local_suppliers.index)
    m.R = pyo.Set(initialize=remote_suppliers.index)

    # define customers
    m.C = pyo.Set(initialize=customers.index)

    # define flowrates
    m.x = pyo.Var(m.R, domain=pyo.NonNegativeReals)
    m.y = pyo.Var(m.C, domain=pyo.NonNegativeReals)
    m.z = pyo.Var(m.L * m.C, domain=pyo.NonNegativeReals)

    m.p = pyo.Param(default=p)

    @m.Objective(sense=pyo.maximize)
    def profit(m):
        return + sum(m.z[l, c]*(customers.loc[c, "price"] - suppliers.loc[l, "cost"]) for l, c in m.L * m.C)                + sum(m.y[c]*customers.loc[c, "price"] for c in m.C)                - sum(m.x[r]*suppliers.loc[r, "cost"] for r in m.R)

    @m.Constraint(m.C)
    def customer_demand(m, c):
        return sum(m.z[l, c] for l in m.L) + m.y[c] <= customers.loc[c, "demand"]

    @m.Constraint()
    def pool_balance(m,):
        return sum(m.x[r] for r in m.R) == sum(m.y[c] for c in m.C)

    @m.Constraint()
    def pool_quality(m):
        return sum(suppliers.loc[r, q] * m.x[r] for r in m.R) == m.p * sum(m.x[r] for r in m.R)

    @m.Constraint(m.C)
    def customer_quality(m, c):
        return m.p * m.y[c] + sum(suppliers.loc[l, q] * m.z[l, c] for l in m.L)                  >= customers.loc[c, q] * (sum(m.z[l, c] for l in m.L) + m.y[c])

    pyo.SolverFactory('cbc').solve(m)

    return m

p = 0.04
m = milk_pooling(p)
print(f"\nprofit = {m.profit():0.2f}")

# The result shows the profit if the pool of milk transported from the remote farms has a milk fat content $p = 0.04$. The profit of 100,000, is better than 81,000 earned for business as usual with just local suppliers, but falls short of the 122,441 earned if the remote milk supply could be transported without pooling.
#
# The following cell presents a full report of the solution.

# In[6]:

def report_solution(m):

    # Supplier report
    S = suppliers.copy()
    for l in m.L:
        for c in m.C:
            S.loc[l, c] = m.z[l, c]()
    for r in m.R:
        S.loc[r, "Pool"] = m.x[r]()
    S = S.fillna(0)
    S["Amount"] = S[m.C].sum(axis=1) + S["Pool"]
    S["Expense"] = S["Amount"]*S["cost"]

    # Customer report
    C = customers.copy()
    for c in m.C:
        for l in m.L:
            C.loc[c, l] = m.z[l, c]()
    for c in m.C:
        C.loc[c, "Pool"] = m.y[c]()
    C = C.fillna(0)
    C["Amount"] = C[m.L].sum(axis=1) + C["Pool"]
    C["fat delivered"] = (sum(C[l]*S.loc[l, "fat"] for l in m.L) + C["Pool"] * m.p())/C["Amount"]
    C["Income"] = C["Amount"] * C["price"]

    print(m)
    print(f"\npool composition = {m.p()}")
    print(f"profit = {m.profit():0.2f}")
    print(f"\nSupplier Report\n")
    display(S.round(4))
    print(f"\nCustomer Report\n")
    display(C.round(4))

report_solution(m)

# ## Profit depends on pool composition $p$

# As this stage the calculations find the maximum profit for a given value of $p$. The challenge, of course, is that the optimal value of $p$ is unknown. The following cell computes profits over a range of $p$.

# In[7]:

p_plot = np.linspace(0.025, 0.055, 200)
profit_plot = [milk_pooling(_).profit() for _ in p_plot]

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(p_plot, profit_plot)
ax.set_title("Milk Pooling")
ax.set_xlabel("Pool composition p")
ax.set_ylabel("Profit")
ax.grid(True)

# The results show the maximum achievable profit with pooling is over than 102,000. What is needed is an optimization technique that can solve for the optimal pool composition and profit. This plot demonstrates how the non-convex bilinear constraints can result in multiple local maxima and a saddle points.

# ## Convex Approximation
#
# Two of the constraints in the p-parameterized milk pooling model have bilinear terms $p y_c$ that are the product of decision variables $p$ and $y_c$ for $c\in C$. Designating the value of these terms as $w_c = p y_c$ the balance on milk fat in the inflow and outflow from the pool becomes a linear constraint
#
# $$
# \begin{align*}
# \sum_{r\in R}\text{conc}_{r}\ x_{r}  & = \sum_{c\in C} w_{c}
# \end{align*}
# $$
#
# and the blending constraint to meet the milk fat requirement for each customer becomes a set of linear constraints
#
# $$
# \begin{align*}
# w_c + \sum_{(l,c)\ \in\ L \times C} \text{conc}_{l}\ z_{l,c}
# & \geq \text{conc}_{c}\ (\sum_{l\in L} z_{l, c} + y_{c})
# & \forall c \in C
# \end{align*}
# $$
#
# What remain, of course, are the bilinear terms $w_c = p y_c$. The values of $y_c$ are bounded between 0 and the demand of customer $c$, and the value of $p$ is bounded between the minimum and maximum milk fat concentrations of the remote farms. The bilinear terms form a set of equality constraints and bounded decision variables.
#
# $$
# \begin{align*}
# & w_c = p y_c & \forall c\in C \\
# \\
# 0 \leq\ & y_c \leq \text{demand}_c\ & \forall c\in C \\
# \min_{r\in R} \text{conc}_r \leq\ & p \leq \max_{r\in R} \text{conc}_r \\
# \end{align*}
# $$
#
# One strategy to create to convert this bilinear problem into a linear problem is to replace the equality constraints $w_c = p y_c$ with a set of linear inequality constraints that keeps the difference $w_c - y p_c$ difference small. This process is a "relaxation" of the bilinear constraint with a goal of creating a problem with better computational characteristics. The process increases the space of admissible solutions which will result in an overly optimistic estimate of profit.
#
# Representing the bounds on $p$ and $y_c$ as
#
# $$
# \begin{align*}
# \underline{p} & \leq p \leq \bar{p} \\
# \underline{y}_c & \leq y_c \leq \bar{y}_c & \forall c\in C
# \end{align*}
# $$
#
# the McCormick envelope on $w_c$ is given by a system of four inequalities. For each $c\in C$,
#
# $$
# \begin{align*}
# w_c & \geq \underline{y}_c p + \underline{p} y_c - \underline{p}\underline{y}_c \\
# w_c & \geq \bar{y}_c p + \bar{p} y_c - \bar{p}\bar{y}_c \\
# w_c & \leq \bar{y}_c p + \underline{p} y_c - \bar{p}\underline{y}_c \\
# w_c & \leq \underline{y}_c p + \bar{p} y_c - \underline{p}\bar{y}_c \\
# \end{align*}
# $$
#
# The following cell implements these constraint in Pyomo. The features to note are:
#
# * Use of a rule to specify bounds on the decision variables `m.y[c]`.
# * Creating a new decision variable `m.p` with bounds.
# * Creating a new set of decision variables `m.w[c]` to replace the bilinear terms.
# * An indexed Pyomo `Block` to bound `m.w[c]` using McCormick envelopes.
#
# The result of these operations is a linear model that will provide an upper bound on the profit. Hopefully the resulting solution and bound will be a close enough approximation to be useful.

# In[8]:

import pyomo.environ as pyo

def milk_pooling_convex(q="fat"):

    m = pyo.ConcreteModel('Milk Pooling Model - Convex Approximation')

    # define sources
    m.L = pyo.Set(initialize=local_suppliers.index)
    m.R = pyo.Set(initialize=remote_suppliers.index)

    # define customers
    m.C = pyo.Set(initialize=customers.index)

    # define flowrates
    m.x = pyo.Var(m.R, domain=pyo.NonNegativeReals)
    m.y = pyo.Var(m.C, bounds=lambda m, c: (0, customers.loc[c, "demand"]))
    m.z = pyo.Var(m.L * m.C, domain=pyo.NonNegativeReals)

    # composition of the pool
    m.p = pyo.Var(bounds=(remote_suppliers["fat"].min(), remote_suppliers["fat"].max()))

    # w[c] = p * y[c]
    m.w = pyo.Var(m.C, domain=pyo.NonNegativeReals)

    @m.Objective(sense=pyo.maximize)
    def profit(m):
        return + sum(m.z[l, c]*(customers.loc[c, "price"] - suppliers.loc[l, "cost"]) for l, c in m.L * m.C)                + sum(m.y[c]*customers.loc[c, "price"] for c in m.C)                - sum(m.x[r]*suppliers.loc[r, "cost"] for r in m.R)

    @m.Block(m.C)
    def mccormick(b, c):
        m = b.model()
        b.ll = pyo.Constraint(expr=m.w[c] >= m.y[c].lower * m.p + m.p.lower * m.y[c] - m.p.lower * m.y[c].lower)
        b.hh = pyo.Constraint(expr=m.w[c] >= m.y[c].upper * m.p + m.p.upper * m.y[c] - m.p.upper * m.y[c].upper)
        b.lh = pyo.Constraint(expr=m.w[c] <= m.y[c].upper * m.p + m.p.lower * m.y[c] - m.p.lower * m.y[c].upper)
        b.hl = pyo.Constraint(expr=m.w[c] <= m.y[c].lower * m.p + m.p.upper * m.y[c] - m.p.upper * m.y[c].lower)

    @m.Constraint(m.C)
    def customer_demand(m, c):
        return sum(m.z[l, c] for l in m.L) + m.y[c] <= customers.loc[c, "demand"]

    @m.Constraint()
    def pool_balance(m,):
        return sum(m.x[r] for r in m.R) == sum(m.y[c] for c in m.C)

    @m.Constraint()
    def pool_quality(m):
        return sum(suppliers.loc[r, q] * m.x[r] for r in m.R) == sum(m.w[c] for c in m.C)

    @m.Constraint(m.C)
    def customer_quality(m, c):
        return m.w[c] + sum(suppliers.loc[l, q] * m.z[l, c] for l in m.L)                  >= customers.loc[c, q] * (sum(m.z[l, c] for l in m.L) + m.y[c])

    pyo.SolverFactory('cbc').solve(m)

    return m

m_convex = milk_pooling_convex()
report_solution(m_convex)

# The convex approximation of the milk pooling model estimates an upper bound on profit of 111,412 for a pool composition $p = 0.040$. The plot below compares this solution to what was found by in an exhaustive search over values of $p$.

# In[9]:

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(p_plot, profit_plot)
ax.set_title("Milk Pooling")
ax.set_xlabel("Pool composition p")
ax.set_ylabel("Profit")
ax.grid(True)

ax.plot(m_convex.p(), m_convex.profit(), 'ro', ms=10)
ax.axhline(m_convex.profit(), color='r', linestyle='--')
ax.axvline(m.p(), color='r', linestyle='--')
ax.annotate("convex approximation", xy=(m_convex.p(), m_convex.profit()),
            xytext=(0.035, 108000), ha="right", fontsize=14,
            arrowprops=dict(shrink=0.1, width=1, headwidth=5))

# The convex approximation is clearly misses the market in the estimate of profit and pool composition $p$. Without the benefit of the full scan of profit as a function of $p$, the only check on the profit estimate would be to compute the solution to model for the reported value of $p$. This is done below.

# In[10]:

m_est = milk_pooling(m_convex.p())

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(p_plot, profit_plot)
ax.set_title("Milk Pooling")
ax.set_xlabel("Pool composition p")
ax.set_ylabel("Profit")
ax.grid(True)

ax.plot(m_convex.p(), m_convex.profit(), 'ro', ms=10)
ax.axhline(m_convex.profit(), color='r', linestyle='--')
ax.axvline(m.p(), color='r', linestyle='--')
ax.annotate("convex approximation", xy=(m_convex.p(), m_convex.profit()),
            xytext=(0.035, 108000), ha="right", fontsize=14,
            arrowprops=dict(shrink=0.1, width=1, headwidth=5))

ax.plot(m_est.p(), m_est.profit(), 'go', ms=10)
ax.axhline(m_est.profit(), color='g', linestyle='--')
ax.annotate("local maxima", xy=(m_convex.p(), m_est.profit()),
            xytext=(0.045, 105000), ha="left", fontsize=14,
            arrowprops=dict(shrink=0.1, width=1, headwidth=5))

# In[11]:

report_solution(m_est)

# With regard to the practical impact, the results of using this particular convex approximation are mixed. The approximation did successfully produce a value for the pool composition $p$ which would produce a profit of over 100,000. However, the reported value for $p$ was actually the smallest of the three local maxima for this problem. This large discrepancy may have large consequences regarding the choice of suppliers.

# ## Nonlinear optimization
#
# The final version of this milk pooling model returns to the bilinear formulation with pool composition $p$ as a decision variable. The following Pyomo implementation needs to specify a solver capable of solving the resulting problem. This has been tested with nonlinear solvers [`ipopt`](https://github.com/coin-or/Ipopt) and [`couenne`](https://github.com/coin-or/Couenne). [Pre-compiled binaries for these solvers can be downloaded from AMPL](https://ampl.com/products/solvers/open-source/).

# In[12]:

import pyomo.environ as pyo

def milk_pooling_bilinear(q="fat"):

    m = pyo.ConcreteModel('Milk Pooling Model - Bilinear')

    # define sources
    m.L = pyo.Set(initialize=local_suppliers.index)
    m.R = pyo.Set(initialize=remote_suppliers.index)

    # define customers
    m.C = pyo.Set(initialize=customers.index)

    # define flowrates
    m.x = pyo.Var(m.R, domain=pyo.NonNegativeReals)
    m.y = pyo.Var(m.C, bounds=lambda m, c: (0, customers.loc[c, "demand"]))
    m.z = pyo.Var(m.L * m.C, domain=pyo.NonNegativeReals)

    m.p = pyo.Var(bounds=(remote_suppliers["fat"].min(), remote_suppliers["fat"].max()))

    @m.Objective(sense=pyo.maximize)
    def profit(m):
        return + sum(m.z[l, c]*(customers.loc[c, "price"] - suppliers.loc[l, "cost"]) for l, c in m.L * m.C)                + sum(m.y[c]*customers.loc[c, "price"] for c in m.C)                - sum(m.x[r]*suppliers.loc[r, "cost"] for r in m.R)

    @m.Constraint(m.C)
    def customer_demand(m, c):
        return sum(m.z[l, c] for l in m.L) + m.y[c] <= customers.loc[c, "demand"]

    @m.Constraint()
    def pool_balance(m,):
        return sum(m.x[r] for r in m.R) == sum(m.y[c] for c in m.C)

    @m.Constraint()
    def pool_quality(m):
        return sum(suppliers.loc[r, q] * m.x[r] for r in m.R) == m.p * sum(m.x[r] for r in m.R)

    @m.Constraint(m.C)
    def customer_quality(m, c):
        return m.p * m.y[c] + sum(suppliers.loc[l, q] * m.z[l, c] for l in m.L)                  >= customers.loc[c, q] * (sum(m.z[l, c] for l in m.L) + m.y[c])

    pyo.SolverFactory('couenne').solve(m)

    return m

# In[13]:

m_global = milk_pooling_bilinear()

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(p_plot, profit_plot)
ax.set_title("Milk Pooling")
ax.set_xlabel("Pool composition p")
ax.set_ylabel("Profit")
ax.grid(True)

ax.plot(m_convex.p(), m_convex.profit(), 'ro', ms=10)
ax.axhline(m_convex.profit(), color='r', linestyle='--')
ax.axvline(m.p(), color='r', linestyle='--')
ax.annotate("convex approximation", xy=(m_convex.p(), m_convex.profit()),
            xytext=(0.035, 108000), ha="right", fontsize=14,
            arrowprops=dict(shrink=0.1, width=1, headwidth=5))

ax.plot(m_est.p(), m_est.profit(), 'go', ms=10)
ax.axhline(m_est.profit(), color='g', linestyle='--')
ax.annotate("local maxima", xy=(m_convex.p(), m_est.profit()),
            xytext=(0.045, 105000), ha="left", fontsize=14,
            arrowprops=dict(shrink=0.1, width=1, headwidth=5))

ax.plot(m_global.p(), m_global.profit(), 'bo', ms=10)
ax.axhline(m_global.profit(), color='b', linestyle='--')
ax.annotate("global maxima", xy=(m_global.p(), m_global.profit()),
            xytext=(0.025, 95000), ha="left", fontsize=14,
            arrowprops=dict(shrink=0.1, width=1, headwidth=5))

report_solution(m_global)

# ## Concluding Remarks
#
# The solution for the bilinear pooling model reveals several features of the problem.
#
# * For the given parameters, pooling raw materials for shipment from remote suppliers yields the most profitable solution, but that solution is only possible because there are local suppliers to augment the pool blend to meet individual customer requirements.
#
# * Customer 2 receives a blend that is 3.3% exceeding the requirement of 3%. This results in some "give away" of product quality in return for the economic benefits of pooling.

# ## Suggested Exercises
#
# This simple model demonstrates practical issues that arise in modeling the non-convex problems. Take time explore the behavior of the model to parameter changes and the nature of the solution.
#
# 1. Examine the model data and explain why the enhanced profits are observed only for a particular range of values in $p$.
#
# 2. Think carefully about the non-convex behavior observed in the plot of profit versus parameter $p$.  Why are the local maxima located where they are, and how are they related to problem data? Test your hypothesis by changing the problem data.  What happens when the product specification for Customer A is set equal to 0.045? What happens when it is set to 0.04?
#
# 3. Revise the Pyomo model using 'cbc', 'gurobi_direct', 'ipopt', and 'bonmin' to find a solution. Did you find a solver that could solve this nonlinear problem?
#
# 4. The above analysis assumed unlimited transport. If the truck turns out to have a limit of 4,000 units of milk, write the mathematical constraint necessary to introduce that limit into the model. Add that constraint to the Pyomo model and discuss the impact on profit.

# In[ ]:

max_stars_repo_name: jckantor/MO-book

max_issues_repo_name: jckantor/MO-book

max_forks_repo_name: jckantor/MO-book

reason: import numpy,import pyomo


text:
"""Creating pedestal results over time"""

import datetime
import os

import matplotlib.dates as md
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from panter.base.mapPerkeo import MapPerkeo
from panter.config import conf_path
from panter.data.dataMisc import FilePerkeo
from panter.data.dataRootPerkeo import RootPerkeo
from panter.data.dataloaderPerkeo import DLPerkeo
from panter.eval.pedPerkeo import PedPerkeo

output_path = os.getcwd()

class PedMapPerkeo(MapPerkeo):
    """Class for creating and handling of drift correction factors.

    Based on base master class MapPerkeo. Can either create from scratch or import a map
    of pedestal results for each pmt over time.

    Parameters
    ----------
    fmeas: np.array(MeasPerkeo)
        Array of data loader output (DLPerkeo).
    bimp_ped: bool
        To import existing map for final pedestal results.
    outfile_flag: str
        Flag which will be added to output file name. Should be chosen according to
        filter in data loader, i.e. fmeas input.

    Attributes
    ----------
    cache: np.array
        Used for storing relevant outputs, besides resulting maps.
        In this case, it is not used.
    maps: list of pd.DataFrame
        map[0]: pedestal map
        Pandas DataFrame with pedestal position and sigma, fit error and rCh2 for each
        PMT with a time stamp. Needs to be imported or calculated.

    Examples
    --------
    Importing existing map with beam flag and plotting result:

    >>> ppm = PedMapPerkeo(filt_meas, bimp_ped=True, outfile_flag="beam")
    >>> ppm()
    >>> ppm.plot_ped_map()

    Starting from scratch:

    >>> file_dir = "/mnt/sda/PerkeoDaten1920/cycle201/cycle201/"
    >>> dataloader = DLPerkeo(file_dir)
    >>> dataloader.auto()
    >>> filt_meas = dataloader.ret_filt_meas(["tp", "src"], [0, 5])
    >>> ppm = PedMapPerkeo(filt_meas, bimp_ped=False, outfile_flag="beam")
    >>> ppm()
    >>> ppm.plot_ped_map()
    """

    def __init__(
        self,
        fmeas: np.array = np.asarray([]),
        bimp_ped: bool = True,
        outfile_flag: str = None,
    ):
        super().__init__(fmeas=fmeas, level=1, bimport=[bimp_ped])
        self.outfile_flag = outfile_flag
        if outfile_flag is None:
            self._outfile = f"ped_map.p"
        else:
            self._outfile = f"ped_map_{outfile_flag}.p"

    def _get_level(self, level: int = 0, bimp: bool = True) -> bool:
        """Try to import and/or calculate given level. Return True/False"""

        if level == 0 and bimp:
            # try to import pedestal map
            impfile = FilePerkeo(f"{conf_path}/{self._outfile}")
            self.maps[level], self.cache = impfile.imp()
            assert self.maps[level].shape[0] > 0, "ERROR: Pedestal map empty."

            return True

        elif level == 0 and not bimp:
            # calc ped map
            self.maps[level] = pd.DataFrame(
                columns=[
                    "time",
                    "ped_list",
                    "ped_err_list",
                    "sig_list",
                    "sig_err_list",
                    "rchi2",
                ]
            )
            self._calc_pedestals()

            return True

        else:
            return False

    def _calc_pedestals(self):
        """Calculate pedestals from measurements"""

        for i, meas in enumerate(self._fmeas):
            print(f"Meas No: {i}")
            if i % 100 == 0 and i > 0:
                self._write_map2file(map_ind=0, fname=self._outfile)

            time = meas.date_list[0]

            data = RootPerkeo(meas.file_list[0])
            pedtest = PedPerkeo(
                dataclass=data,
                bplot_res=False,
                bplot_fit=False,
                bplot_log=False,
            )
            ped_list = pedtest.ret_pedestals().T[0]
            ped_err_list = pedtest.ret_pedestals().T[1]
            sig_list = pedtest.ret_pedestals().T[2]
            sig_err_list = pedtest.ret_pedestals().T[3]
            rchi2 = pedtest.ret_pedestals().T[4]

            meas_dict = {
                "time": time,
                "ped_list": ped_list,
                "ped_err_list": ped_err_list,
                "sig_list": sig_list,
                "sig_err_list": sig_err_list,
                "rchi2": rchi2,
            }
            self.maps[0] = self.maps[0].append(meas_dict, ignore_index=True)

        assert (
            self._write_map2file(map_ind=0, fname=self._outfile) == 0
        ), "ERROR: Export of pedestal map failed."

        return 0

    def plot_ped_map(self, bsave: bool = False):
        """Plot Pedestal map for all PMTs"""

        fig, axs = plt.subplots(1, 2, figsize=(20, 12))
        fig.suptitle(f"Pedestal Map {self.outfile_flag}")

        axs[0].set_title("Pedestals over time Det 0")
        axs[1].set_title("Peak pos over time Det 1")
        axs.flat[0].set(xlabel="Time [D - M ]", ylabel="Ped pos [ch]")
        axs.flat[1].set(xlabel="Time [D - M ]", ylabel="Ped pos [ch]")

        xfmt = md.DateFormatter("%m-%d\n%H:%M")
        axs[0].xaxis.set_major_formatter(xfmt)
        axs[1].xaxis.set_major_formatter(xfmt)

        ped_df = self.maps[0]["ped_list"].apply(pd.Series)
        ped_err_df = self.maps[0]["ped_err_list"].apply(pd.Series)

        dates_plot = [datetime.datetime.fromtimestamp(t) for t in self.maps[0]["time"]]

        for PMT in range(8):
            axs[0].errorbar(
                dates_plot,
                ped_df[PMT],
                yerr=ped_err_df[PMT],
                fmt=".",
                label=f"PMT{PMT}",
            )
            axs[1].errorbar(
                dates_plot,
                ped_df[PMT + 8],
                yerr=ped_err_df[PMT + 8],
                fmt=".",
                label=f"PMT{PMT + 8}",
            )
        axs[0].legend()
        axs[1].legend()
        if bsave:
            plt.savefig(output_path + "/" + self._outfile[:-1] + "png", dpi=300)
        plt.show()

        return 0

if False:
    # file_dir = "/mnt/sda/PerkeoDaten1920/cycle201/cycle201/"
    # dataloader = DLPerkeo(file_dir)
    # dataloader.auto()
    # filt_meas = dataloader.ret_filt_meas(["tp", "src"], [1, 4])
    # filt_meas = dataloader.ret_filt_meas(["tp", "src", "nomad_no"], [1, 3, 67732])
    sources = np.asarray(["beam", "ce", "cd", "cs", "bi", "sn"])
    for src in sources:
        ppm = PedMapPerkeo([], bimp_ped=True, outfile_flag=src)
        ppm()
        ppm.plot_ped_map(bsave=True)

if False:
    ind = 1
    file_dir = "/mnt/sda/PerkeoDaten1920/cycle201/cycle201/"
    dataloader = DLPerkeo(file_dir)
    dataloader.auto()
    filt_meas = dataloader.ret_filt_meas(["tp", "src"], [1, ind])
    sources = np.asarray(["cd", "ce", "bi", "sn", "cs", "beam"])
    src = sources[ind]
    ppm = PedMapPerkeo(filt_meas, bimp_ped=False, outfile_flag=src)
    ppm()
    ppm.plot_ped_map(bsave=True)

max_stars_repo_name: maxlampe/panter

max_issues_repo_name: maxlampe/panter

max_forks_repo_name: maxlampe/panter


text:
#!/usr/bin/env python
import sys
import numpy
import numpy.linalg
import time

sys.path = ['../..'] + sys.path
from vector_forward_mode import *

def speelpenning(x):
	return numpy.prod(x)

N = 100
P = N*(N+1)/2
D = 3

x = numpy.zeros(N)
S = numpy.zeros((N,P))

ax = double_to_adouble(x,S,D)

start_time = time.time()
speelpenning(ax)
end_time =time.time()

print('needed time was %0.6f seconds!'%(end_time-start_time))