text:
from __future__ import division

import unittest

import numpy as np

from wyrm.types import Data
from wyrm.processing import append_cnt
from wyrm.processing import swapaxes

class TestAppendCnt(unittest.TestCase):

    def setUp(self):
        ones = np.ones((10, 5))
        # cnt with 1, 2, 3
        cnt = np.append(ones, ones*2, axis=0)
        cnt = np.append(cnt, ones*3, axis=0)
        channels = ['ca1', 'ca2', 'cb1', 'cb2', 'cc1']
        time = np.linspace(0, 3000, 30, endpoint=False)
        self.dat = Data(cnt, [time, channels], ['time', 'channel'], ['ms', '#'])
        self.dat.markers = [[0, 'a'], [1, 'b']]
        self.dat.fs = 10

    def test_append_cnt(self):
        """append_cnt."""
        dat = append_cnt(self.dat, self.dat)
        self.assertEqual(dat.data.shape[0], 2*self.dat.data.shape[0])
        self.assertEqual(len(dat.axes[0]), 2*len(self.dat.axes[0]))
        np.testing.assert_array_equal(dat.data, np.concatenate([self.dat.data, self.dat.data], axis=0))
        np.testing.assert_array_equal(dat.axes[0], np.linspace(0, 6000, 60, endpoint=False))
        self.assertEqual(dat.markers, self.dat.markers + [[x[0] + 3000, x[1]] for x in self.dat.markers])

    def test_append_cnt_with_extra(self):
        """append_cnt with extra must work with list and ndarrays."""
        self.dat.a = list(range(10))
        self.dat.b = np.arange(10)
        dat = append_cnt(self.dat, self.dat, extra=['a', 'b'])
        self.assertEqual(dat.a, list(range(10)) + list(range(10)))
        np.testing.assert_array_equal(dat.b, np.concatenate([np.arange(10), np.arange(10)]))
        self.assertEqual(dat.markers, self.dat.markers + [[x[0] + 3000, x[1]] for x in self.dat.markers])

    def test_append_cnt_swapaxes(self):
        """append_cnt must work with nonstandard timeaxis."""
        dat = append_cnt(swapaxes(self.dat, 0, 1), swapaxes(self.dat, 0, 1), timeaxis=1)
        dat = swapaxes(dat, 0, 1)
        dat2 = append_cnt(self.dat, self.dat)
        self.assertEqual(dat, dat2)

    def test_append_cnt_copy(self):
        """append_cnt means must not modify argument."""
        cpy = self.dat.copy()
        append_cnt(self.dat, self.dat)
        self.assertEqual(self.dat, cpy)

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

max_stars_repo_name: jscastanoc/wyrm

max_issues_repo_name: jscastanoc/wyrm

max_forks_repo_name: jscastanoc/wyrm


text:
#!/usr/bin/env/python
##============================================================

import os
import glob
import numpy
import scipy
import pylab
import Image
import ImageEnhance
import pyfits
import getopt
from sys import argv
from OSMethods import file_seek

##============================================================

def main():
	"""
	NAME
	  ImageMethods.py

	PURPOSE
	  A variety of image manipulation methods -- filetype conversions,
	  cropping function, image anhancement, and stretching.

	COMMENTS

	USAGE
	  ImageMethods.py [flags] function args

	  * args depend on function

	FLAGS
	  -h            Print this message

	INPUTS
	  Directory and/or file name depending on function

	OPTIONAL INPUTS

	OUTPUTS
	  stdout        Useful information
	  outfile       When converting/modifying images
	  array					When modifying image data

	EXAMPLES
	  ImageMethods.py F2P ../Images smileyimage?
	  ImageMethods.py cropwhite ./whiteborderimg.png

	BUGS
		No bugs per se, but a lot of usability improvements / tidying up
		to do.

	HISTORY
	  2012-06-12 started Sandford (NYU)
	"""

	## -------------------------------------------------------------------
	## Options and help-string

	try:
		opts, args = getopt.getopt(argv, "h",["help"])
	except getopt.GetoptError, err:
		## Print help information and exit:
		print str(err)
		print main.__doc__
		return

	for o,a in opts:
		print o,a
		if o in ("-h", "--help"):
			print main.__doc__
			return

	## -------------------------------------------------------------------

	if args[1]=="F2P": FITStoPNG(args)
	elif args[1]=="P2P": PSFtoPNG_all(args)
	elif args[1]=="cropwhite": pngcropwhite(args[2])
	elif args[1]=="brighten": bright_all()
	elif args[1]=="stretch": stretch(args[2])
	elif args[1]=="stretchall": stretch_all(0)
	elif args[1]=="Stretchall": stretch_all(1)
	else: print "ImageScript.py: select valid option:\
							\n\tF2P\n\tP2P\n\tbrighten\n\tstretch\n\tcropwhite"

	return

##============================================================

##============================================================
##============================================================
## IMAGE CONVERSION
##============================================================
##============================================================

## Converts FITS images in subdirectories of argv[2] to PNG format
  ## argv[3] is string (e.g. "CFHTLS_")
def FITStoPNG(args):
  i=0
  for f_img in file_seek(args[2], args[3]+"*.fits"):
		## Prepare destination directory
		outdir=os.path.dirname(f_img)+"/Enhanced_Original/" ## Needs to be commented when brightening
		if os.path.isdir(outdir)==False: os.mkdir(outdir)
		## Outfile name
		p_img = outdir+"/"+os.path.basename(f_img)[:-5]+".png"
		## File conversion
		#os.system("ds9 "+f_img+" -minmax -colorbar no -saveimage png "+p_img+" -exit")
		try:
			hdulist = pyfits.open(f_img)
			scipy.misc.imsave(f_img[:-4]+"png", hdulist[0].data)
			i+=1
		except:
			pass
  print "ImageScript.py: FITStoPNG: Converted",i,"files to png format."
  return

##============================================================

## Converts all .psf to .png in subdirectories of argv[2]
def PSFtoPNG_all(args):
	for psffile in file_seek(args[2], "*.psf"):
		PSFtoPNG(psffile)
	return

##============================================================

## Convert .psf to a number of .png

def PSFtoPNG(imag_in):

	## Prepare directory for output
	out_dir=os.path.dirname(imag_in)+"/PSF_PNG_Breakdown/"
	if os.path.isdir(out_dir)==False: os.mkdir(out_dir)

	## Load image
	IMG = pyfits.open(imag_in)

	## Extract image data
	psfdata = IMG[1].data[0][0]
	Nim = IMG[1].header["psfaxis3"]
	Xdim = IMG[1].header["psfaxis1"]
	Ydim = IMG[1].header["psfaxis2"]

	## Reshape to separate images
	psfdata = numpy.reshape(psfdata, [Nim,Xdim,Ydim])#, dtype="unit8"

	## Load data into new image to save individually
	for i in range (Nim):
		## Outfile name
		imag_out = out_dir+os.path.basename(imag_in)[:-4]+"_fit"+str(i)+".png"
		#image = Image.new("L",[Xdim,Ydim],None)
		#image.putdata(psfdata[i,:,:],1.0)
		#image = Image.frombuffer('L',(Xdim,Ydim),psfdata,'raw','L',0,1)
		#image.save(imag_out)
		scipy.misc.imsave(imag_out, psfdata[i,:,:])

	print "ImageScript.py: PSFtoPNG: made",Nim,"PNG files from",os.path.basename(imag_in)

	return

##============================================================

## Convert from .ext to .jpg
	## Call from Python
def EXTtoJPG(infile):
  f,e=os.path.splitext(infile)
  outfile=f+".jpg"
  if infile!=outfile:
    try:
      Image.open(infile).save(outfile)
    except IOError:
      print "ImageScript.py: EXTtoJPG: cannot convert", infile
  return

##============================================================

## Crop white pixels from image border
	## Expects a pixel-array with 255-pixels around central image

def pngcropwhite(imgarr, level=254.9):

	## Check whether imgarr is filename or array
	if isinstance(imgarr, numpy.ndarray):
		returntype="array"
	elif isinstance(imgarr, str):
		filename = imgarr
		imgarr = png_pix(imgarr)
		returntype="pngfile"

	## Original dimensions and centre
	img_h,img_w = imgarr.shape
	crow,ccol = int(0.5*img_h), int(0.5*img_w)

	## Find dimensions of the central image by stepping inwards from periphery
	for i in range (img_h):
		if imgarr[i,ccol]<level:
			top=i
			break
	for i in range (img_h-1,0,-1):
		if imgarr[i,ccol]<level:
			bottom=i+1
			break
	for i in range (img_w):
		if imgarr[crow,i]<level:
			left=i
			break
	for i in range (img_w-1,0,-1):
		if imgarr[crow,i]<level:
			right=i+1
			break

	## Retain slice of interest
	imgarr = imgarr[top:bottom, left:right]

	## Output depends on input -- do we want an array or a file?
	if returntype is "array":
		return imgarr[top:bottom, left:right]
	elif returntype is "pngfile":
		scipy.misc.imsave(filename, imgarr)
		return

##============================================================
## Image to array conversion

## Read in .fits and return array of (decimal) pixel values
	## Works for PSFEx snap - other fits have different heirarchy
def fits_pix(filename):
	im = pyfits.open(filename)
	return numpy.array(im[0].data)

def fits_hdr(filename):
	im = pyfits.open(filename)
	return im[0].header

## Convert png file to an array of (0-255) pixel values
def png_pix(filename):
	img = Image.open(filename)
	dim = [img.size[1],img.size[0]]	## Because of course they need to be reversed
	img = img.convert("L")
	dat = numpy.array(img.getdata())
	dat = numpy.reshape(dat, dim)
	return dat

##============================================================

## Make FITS image from array of pixels
def makeFITS(outfile, array, head=None):
	hdu = pyfits.PrimaryHDU(array[::-1])
	if head is not None:
		hdu.header = head
	hdu.writeto(outfile)
	return

##============================================================
## Pixel file to PNG
## Greyscale

def pixfiletopng(inname):
	DAT = numpy.loadtxt(inname)
	outname = os.path.splitext(inname)[0]+".png"
	scipy.misc.imsave(outname,DAT)
	return

##============================================================
##============================================================
## IMAGE ENHANCEMENT
##============================================================
##============================================================

## Automate image stretching
	## sel=0: Stretch original images only
	## sel=1: Stretch original images, and apply the same stretch to diagnostic images

def stretch_all(sel=0):

	for string in ["CFHTLS_*_sci.fits"]:#
		for fits in file_seek("../Data/",string):
			if sel==0: stretch(fits)
			else: stretch_plus(fits)

	return

##============================================================

## Stretches a file and then its children

def stretch_plus(fits_original):

	## Stretching parameters
	vmin,vmax=stretch(fits_original)

	## Apply to diagnostic files
	homedir=os.path.dirname(fits_original)
	for string in ["snap_*_sci.fits","samp_*_sci.fits","resi_*_sci.fits"]:#"CFHTLS_*_sci.fits"
		for fits in file_seek(homedir,string):
			stretch(fits, vmin,vmax)

	return

##============================================================

def stretch(fits_original, vmin=None,vmax=None):

	## Prepare destination directory
	outdir=os.path.dirname(fits_original)+"/Enhanced_Original/"
	if os.path.isdir(outdir)==False: os.mkdir(outdir)

	## Outfile name
	outfile = outdir+os.path.basename(fits_original)[:-5]+"_rescaled.fits"

	## Load image
	IMG = pyfits.open(fits_original)
	## Array of pixel values
	im_data = IMG[0].data

	## If we haven't provided stretch parameters, calculate them
	if vmin==None or vmax==None: vmin,vmax = stretch_params(im_data)

	## Rescale
	im_data = linear_rescale(im_data, vmin,vmax)

	## Update FITS data
	IMG[0].data = im_data
	## Delete pre-existing files
	if os.path.isfile(outfile): os.remove(outfile)
	## Save
	IMG.writeto(outfile)
	IMG.close()

	return vmin,vmax

##============================================================
## Linear scaling of input array
def linear_rescale(inputArray, scale_min=None, scale_max=None):

	imageData=numpy.array(inputArray, copy=True)

	if scale_min == None:
		scale_min = imageData.min()
	if scale_max == None:
		scale_max = imageData.max()

	imageData = imageData.clip(min=scale_min, max=scale_max)
	imageData = (imageData-scale_min)/(scale_max-scale_min)
	imageData[imageData<0.] = 0.0
	imageData[imageData>1.] = 1.0

	return imageData

##============================================================

## Parameters for re-scale/stretch
def stretch_params(array):
	## Find median and 97.5% value
	a=numpy.median(array)
	b=numpy.sort(array.flatten())[0.975*len(array.flatten())]
	## Stretch parameters
	sigma = 0.5*(b-a)
	vmin = a-5.0*sigma
	vmax = a+10.0*sigma
	return vmin,vmax

##============================================================
## Standard enhancements

## Contrast
def cont(imagefile, factor=10):
	im=Image.open(imagefile)
	enhancer = ImageEnhance.Contrast(im)
	enhancer.enhance(factor).show("Contrast %f" % factor)
	return
## Brightness
def brig(imagefile,factor=10,brightfile=None):
	im=Image.open(imagefile)
	enhancer = ImageEnhance.Brightness(im)
	if brightfile==None: brightfile=imagefile[:-4]+"_br"+str(factor)+".png"
	else: brightfile=os.path.dirname(imagefile)+"/Enhanced_Original/"+os.path.basename(imagefile)[:-4]+"_br"+str(factor)+".png"
	enhancer.enhance(factor).save(brightfile)
	return

##============================================================

## Brighten ALL originals

def bright_all(factor=8):
	for original in file_seek("../Data/","CFHTLS_*sci.png"):
		brig(original,factor,1)
	return

##============================================================
##============================================================
## Miscellanious
##============================================================
##============================================================

## Write 2D array to a file in easy-parsed fashion
def writearray(f,arr,close=True,dp=5):
	for i in range (arr.shape[0]):
		for j in range (arr.shape[1]):
			f.write(str(round(arr[i,j],dp))+"\t")
		f.write("\n")
	if close == True: f.close()
	return

##============================================================
if __name__=="__main__":
  main()

ext: py

max_stars_repo_name: ccld/HumVI

max_issues_repo_name: ccld/HumVI

max_forks_repo_name: ccld/HumVI

reason: import numpy,import scipy


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

import click

import datetime
import pandas as pd

import numpy as np
from matplotlib import pyplot as plt
import matplotlib.cm as cm

import windrose  # noqa

pd.set_option('max_rows', 10)

def get_by_func(by=None, by_func=None):
    if by is None and by_func is None:
        by = 'MS'

    if by in ['year', 'yearly', 'Y']:
        return lambda dt: dt.year
    elif by in ['month', 'monthly', 'MS']:  # MS: month start
        return lambda dt: (dt.year, dt.month)
    elif by in ['day', 'daily', 'D']:
        return lambda dt: (dt.year, dt.month, dt.day)
    elif by is None and by_func is not None:
        return by_func
    else:
        raise NotImplementedError("'%s' is not an allowed 'by' parameter" % by)

def tuple_position(i, nrows, ncols):
    i_sheet, sheet_pos = divmod(i, ncols * nrows)
    i_row, i_col = divmod(sheet_pos, ncols)
    return i_sheet, i_row, i_col

@click.command()
@click.option("--filename", default="samples/sample_wind_poitiers.csv", help="Input filename")
@click.option("--year", default=2014, help="Year")
def main(filename, year):
    df_all = pd.read_csv(filename, parse_dates=['Timestamp'])
    df_all = df_all.set_index('Timestamp')

    f_year = get_by_func('year')
    df_all['by_page'] = df_all.index.map(f_year)
    f_month = get_by_func('month')
    df_all['by'] = df_all.index.map(f_month)

    df_all = df_all.reset_index().set_index(['by_page', 'by', 'Timestamp'])

    print(df_all)

    nrows, ncols = 3, 4
    fig = plt.figure()
    bins = np.arange(0.01, 8, 1)

    fig.suptitle("Wind speed - %d" % year)
    for month in range(1, 13):
        ax = fig.add_subplot(nrows, ncols, month, projection='windrose')
        title = datetime.datetime(year, month, 1).strftime("%b")
        ax.set_title(title)
        try:
            df = df_all.loc[year].loc[(year, month)]
        except KeyError:
            continue
        direction = df['direction'].values
        var = df['speed'].values
        # ax.contour(direction, var, bins=bins, colors='black', lw=3)
        ax.contourf(direction, var, bins=bins, cmap=cm.hot)
        ax.contour(direction, var, bins=bins, colors='black')

    plt.show()

if __name__ == "__main__":
    main()

max_stars_repo_name: jokva/windrose

max_issues_repo_name: jokva/windrose

max_forks_repo_name: jokva/windrose

reason: import numpy


text:
from copro import machine_learning, conflict, utils, evaluation, data
import pandas as pd
import numpy as np
import pickle
import os, sys

def all_data(X, Y, config, scaler, clf, out_dir, run_nr):
    """Main model workflow when all XY-data is used.
    The model workflow is executed for each classifier.

    Args:
        X (array): array containing the variable values plus IDs and geometry information.
        Y (array): array containing merely the binary conflict classifier data.
        config (ConfigParser-object): object containing the parsed configuration-settings of the model.
        scaler (scaler): the specified scaling method instance.
        clf (classifier): the specified model instance.
        out_dir (str): path to output folder.

    Returns:
        dataframe: containing the test-data X-array values.
        datatrame: containing model output on polygon-basis.
        dict: dictionary containing evaluation metrics per simulation.
    """
    if config.getboolean('general', 'verbose'): print('DEBUG: using all data')

    # split X into training-set and test-set, scale training-set data
    X_train, X_test, y_train, y_test, X_train_geom, X_test_geom, X_train_ID, X_test_ID = machine_learning.split_scale_train_test_split(X, Y, config, scaler)

    # convert to dataframe
    X_df = pd.DataFrame(X_test)

    # fit classifier and make prediction with test-set
    y_pred, y_prob = machine_learning.fit_predict(X_train, y_train, X_test, clf, config, out_dir, run_nr)
    y_prob_0 = y_prob[:, 0] # probability to predict 0
    y_prob_1 = y_prob[:, 1] # probability to predict 1

    # evaluate prediction and save to dict
    eval_dict = evaluation.evaluate_prediction(y_test, y_pred, y_prob, X_test, clf, config)

    # aggregate predictions per polygon
    y_df = conflict.get_pred_conflict_geometry(X_test_ID, X_test_geom, y_test, y_pred, y_prob_0, y_prob_1)

    return X_df, y_df, eval_dict

def leave_one_out(X, Y, config, scaler, clf, out_dir):
    """Model workflow when each variable is left out from analysis once.
    Output is limited to the metric scores.
    Output is stored to sub-folders of the output directory, with each sub-folder containing output for a n-1 variable combination.
    After computing metric scores per prediction (i.e. n-1 variables combinations), model exit is forced.
    Not tested yet for more than one simulation!

    Args:
        X (array): array containing the variable values plus unique identifer and geometry information.
        Y (array): array containing merely the binary conflict classifier data.
        config (ConfigParser-object): object containing the parsed configuration-settings of the model.
        scaler (scaler): the specified scaling method instance.
        clf (classifier): the specified model instance.
        out_dir (str): path to output folder.

    Raises:
        DeprecationWarning: this function will most likely be deprecated due to lack of added value and applicability.
    """

    raise DeprecationWarning('WARNING: the leave-one-out model is not supported anymore and will be deprecated in a future release')

    X_train, X_test, y_train, y_test, X_train_geom, X_test_geom, X_train_ID, X_test_ID = machine_learning.split_scale_train_test_split(X, Y, config, scaler)

    for i, key in zip(range(X_train.shape[1]), config.items('data')):

        print('INFO: removing data for variable {}'.format(key[0]))

        X_train_loo = np.delete(X_train, i, axis=1)
        X_test_loo = np.delete(X_test, i, axis=1)

        sub_out_dir = os.path.join(out_dir, '_only_'+str(key[0]))
        if not os.path.isdir(sub_out_dir):
            os.makedirs(sub_out_dir)

        y_pred, y_prob = machine_learning.fit_predict(X_train_loo, y_train, X_test_loo, clf, config)

        eval_dict = evaluation.evaluate_prediction(y_test, y_pred, y_prob, X_test_loo, clf, config)

        utils.save_to_csv(eval_dict, sub_out_dir, 'evaluation_metrics')

    sys.exit('INFO: leave-one-out model execution stops here.')

def single_variables(X, Y, config, scaler, clf, out_dir):
    """Model workflow when the model is based on only one single variable.
    Output is limited to the metric scores.
    Output is stored to sub-folders of the output directory, with each sub-folder containing output for a 1 variable combination.
    After computing metric scores per prediction (i.e. per variable), model exit is forced.
    Not tested yet for more than one simulation!

    Args:
        X (array): array containing the variable values plus unique identifer and geometry information.
        Y (array): array containing merely the binary conflict classifier data.
        config (ConfigParser-object): object containing the parsed configuration-settings of the model.
        scaler (scaler): the specified scaling method instance.
        clf (classifier): the specified model instance.
        out_dir (str): path to output folder.

    Raises:
        DeprecationWarning: this function will most likely be deprecated due to lack of added value and applicability.
    """

    raise DeprecationWarning('WARNING: the single-variable model is not supported anymore and will be deprecated in a future release')

    X_train, X_test, y_train, y_test, X_train_geom, X_test_geom, X_train_ID, X_test_ID = machine_learning.split_scale_train_test_split(X, Y, config, scaler)

    for i, key in zip(range(X_train.shape[1]), config.items('data')):

        print('INFO: single-variable model with variable {}'.format(key[0]))

        X_train_svmod = X_train[:, i].reshape(-1, 1)
        X_test_svmod = X_test[:, i].reshape(-1, 1)

        sub_out_dir = os.path.join(out_dir, '_excl_'+str(key[0]))
        if not os.path.isdir(sub_out_dir):
            os.makedirs(sub_out_dir)

        y_pred, y_prob = machine_learning.fit_predict(X_train_svmod, y_train, X_test_svmod, clf, config)

        eval_dict = evaluation.evaluate_prediction(y_test, y_pred, y_prob, X_test_svmod, clf, config)

        utils.save_to_csv(eval_dict, sub_out_dir, 'evaluation_metrics')

    sys.exit('INFO: single-variable model execution stops here.')

def dubbelsteen(X, Y, config, scaler, clf, out_dir):
    """Model workflow when the relation between variables and conflict is based on randomness.
    Thereby, the fraction of actual conflict is equal to observations, but the location in array is randomized by shuffling.
    The model workflow is executed for each model simulation.

    Args:
        X (array): array containing the variable values plus unique identifer and geometry information.
        Y (array): array containing merely the binary conflict classifier data.
        config (ConfigParser-object): object containing the parsed configuration-settings of the model.
        scaler (scaler): the specified scaling method instance.
        clf (classifier): the specified model instance.
        out_dir (str): path to output folder.

    Returns:
        dataframe: containing the test-data X-array values.
        datatrame: containing model output on polygon-basis.
        dict: dictionary containing evaluation metrics per simulation.
    """

    print('INFO: dubbelsteenmodel running')
    raise DeprecationWarning('WARNING: the dubbelsteenmodel model is not supported anymore and will be deprecated in a future release')

    Y = utils.create_artificial_Y(Y)

    X_train, X_test, y_train, y_test, X_train_geom, X_test_geom, X_train_ID, X_test_ID = machine_learning.split_scale_train_test_split(X, Y, config, scaler)

    y_pred, y_prob = machine_learning.fit_predict(X_train, y_train, X_test, clf, config)

    eval_dict = evaluation.evaluate_prediction(y_test, y_pred, y_prob, X_test, clf, config)

    y_df = conflict.get_pred_conflict_geometry(X_test_ID, X_test_geom, y_test, y_pred)

    X_df = pd.DataFrame(X_test)

    return X_df, y_df, eval_dict

def predictive(X, clf, scaler, config):
    """Predictive model to use the already fitted classifier to make annual projections for the projection period.
    As other models, it reads data which are then scaled and used in conjuction with the classifier to project conflict risk.

    Args:
        X (array): array containing the variable values plus unique identifer and geometry information.
        clf (classifier): the fitted specified classifier instance.
        scaler (scaler): the fitted specified scaling method instance.
        config (ConfigParser-object): object containing the parsed configuration-settings of the model for a projection run.

    Returns:
        datatrame: containing model output on polygon-basis.
    """

    # splitting the data from the ID and geometry part of X
    X_ID, X_geom, X_data = conflict.split_conflict_geom_data(X.to_numpy())

    # transforming the data
    # fitting is not needed as already happend before
    if config.getboolean('general', 'verbose'): print('DEBUG: transforming the data from projection period')
    X_ft = scaler.transform(X_data)

    # make projection with transformed data
    if config.getboolean('general', 'verbose'): print('DEBUG: making the projections')
    y_pred = clf.predict(X_ft)

    # predict probabilites of outcomes
    y_prob = clf.predict_proba(X_ft)
    y_prob_0 = y_prob[:, 0] # probability to predict 0
    y_prob_1 = y_prob[:, 1] # probability to predict 1

    # stack together ID, gemoetry, and projection per polygon, and convert to dataframe
    arr = np.column_stack((X_ID, X_geom, y_pred, y_prob_0, y_prob_1))
    y_df = pd.DataFrame(arr, columns=['ID', 'geometry', 'y_pred', 'y_prob_0', 'y_prob_1'])

    return y_df

max_issues_repo_name: JannisHoch/copro


text:
import argparse
import os, time
import torch
torch.cuda.empty_cache()
import shutil
import numpy as np
import torch.nn as nn
import torchvision.transforms as transforms
import torchvision.datasets as datasets
import torch.nn.functional as F
import torch.optim as optim
from ptflops import get_model_complexity_info
from torch.autograd import Variable
import models
import optim
import torch.backends.cudnn as cudnn
from torch_poly_lr_decay import PolynomialLRDecay
import torchsummary
from deepshift.convert import convert_to_shift, round_shift_weights, count_layer_type
import distutils
import matplotlib.pyplot as plt
from adder import adder as adder_fast
import deepshift
from pathlib import Path
import collections
from collections import OrderedDict

# Training settings
parser = argparse.ArgumentParser(description='PyTorch AdderNet Trainning')
parser.add_argument('--dataset', type=str, default='imagenet', help='training dataset')
parser.add_argument('--data_path', type=str, default=None, help='path to dataset')
parser.add_argument('--batch_size', type=int, default=256, metavar='N', help='batch size for training')
parser.add_argument('--test_batch_size', type=int, default=256, metavar='N', help='batch size for testing')
parser.add_argument('--epochs', type=int, default=5, metavar='N', help='number of epochs to train')
parser.add_argument('--start_epoch', type=int, default=0, metavar='N', help='restart point')
parser.add_argument('--schedule', type=int, nargs='+', default=[80, 120], help='learning rate schedule')
parser.add_argument('--lr', type=float, default=0.1, metavar='LR', help='learning rate')
parser.add_argument('--lr-sign', default=None, type=float, help='separate initial learning rate for sign params')
parser.add_argument('--lr_decay', default='stepwise', type=str, choices=['stepwise', 'poly'])
parser.add_argument('--optimizer', type=str, default='sgd', help='used optimizer')
parser.add_argument('--momentum', type=float, default=0.9, metavar='M', help='SGD momentum')
parser.add_argument('--weight-decay', '--wd', default=1e-4, type=float, metavar='W', help='weight decay')
parser.add_argument('--resume', default='', type=str, metavar='PATH', help='path to latest checkpoint')
parser.add_argument('--seed', type=int, default=1, metavar='S', help='random seed')
parser.add_argument('--save', default='./logs', type=str, metavar='PATH', help='path to save prune model')
parser.add_argument('--save_suffix', default='retrain_', type=str, help='identify which retrain')
parser.add_argument('--arch', default='resnet20', type=str, help='architecture to use')
parser.add_argument('--no-cuda', action='store_true', default=False, help='disables CUDA training')
parser.add_argument('--log_interval', type=int, default=100, metavar='N', help='how many batches to wait before logging training status')
# multi-gpus
parser.add_argument('--gpu_ids', type=str, default='0', help='gpu ids: e.g. 0  0,1,2, 0,2. use -1 for CPU')
# shift hyper-parameters
parser.add_argument('--shift_depth', type=int, default=0, help='how many layers to convert to shift')
parser.add_argument('--shift_type', type=str, choices=['Q', 'PS'], help='shift type for representing weights')
parser.add_argument('--rounding', default='deterministic', choices=['deterministic', 'stochastic'])
parser.add_argument('--weight_bits', type=int, default=5, help='number of bits to represent the shift weights')
parser.add_argument('--use-kernel', type=lambda x:bool(distutils.util.strtobool(x)), default=False, help='whether using custom shift kernel')
# add hyper-parameters
parser.add_argument('--add_quant', type=bool, default=False, help='whether to quantize adder layer')
parser.add_argument('--add_bits', type=int, default=8, help='number of bits to represent the adder filters')
# visualization
parser.add_argument('--visualize', action="store_true", default=False, help='if use visualization')
# reinit
parser.add_argument('--reinit', type=str, default=None, help='whether reinit or finetune')

args = parser.parse_args()
args.cuda = not args.no_cuda and torch.cuda.is_available()

torch.manual_seed(args.seed)
if args.cuda:
    torch.cuda.manual_seed(args.seed)

if not os.path.exists(args.save):
    os.makedirs(args.save)

cudnn.benchmark = True

gpu = args.gpu_ids
gpu_ids = args.gpu_ids.split(',')
args.gpu_ids = []
for gpu_id in gpu_ids:
    id = int(gpu_id)
    args.gpu_ids.append(id)
print(args.gpu_ids)
# if len(args.gpu_ids) > 0:
#    torch.cuda.set_device(args.gpu_ids[0])

kwargs = {'num_workers': 1, 'pin_memory': True} if args.cuda else {}
if args.dataset == 'cifar10':
    train_loader = torch.utils.data.DataLoader(
        datasets.CIFAR10('./data.cifar10', train=True, download=True,
                       transform=transforms.Compose([
                           transforms.Pad(4),
                           transforms.RandomCrop(32),
                           transforms.RandomHorizontalFlip(),
                           transforms.ToTensor(),
                           transforms.Normalize((0.4914, 0.4822, 0.4465), (0.2023, 0.1994, 0.2010))
                       ])),
        batch_size=args.batch_size, shuffle=True, **kwargs)
    test_loader = torch.utils.data.DataLoader(
        datasets.CIFAR10('./data.cifar10', train=False, transform=transforms.Compose([
                           transforms.ToTensor(),
                           transforms.Normalize((0.4914, 0.4822, 0.4465), (0.2023, 0.1994, 0.2010))
                       ])),
        batch_size=args.test_batch_size, shuffle=True, **kwargs)
elif args.dataset == 'cifar100':
    train_loader = torch.utils.data.DataLoader(
        datasets.CIFAR100('./data.cifar100', train=True, download=True,
                       transform=transforms.Compose([
                           transforms.Pad(4),
                           transforms.RandomCrop(32),
                           transforms.RandomHorizontalFlip(),
                           transforms.ToTensor(),
                           transforms.Normalize((0.4914, 0.4822, 0.4465), (0.2023, 0.1994, 0.2010))
                       ])),
        batch_size=args.batch_size, shuffle=True, **kwargs)
    test_loader = torch.utils.data.DataLoader(
        datasets.CIFAR100('./data.cifar100', train=False, transform=transforms.Compose([
                           transforms.ToTensor(),
                           transforms.Normalize((0.4914, 0.4822, 0.4465), (0.2023, 0.1994, 0.2010))
                       ])),
        batch_size=args.test_batch_size, shuffle=True, **kwargs)
elif args.dataset == 'mnist':
    trainset = datasets.MNIST('../MNIST', download=True, train=True,
        transform=transforms.Compose([
            transforms.ToTensor(),
            transforms.Normalize((0.1307,), (0.3081,))]
        )
    )
    train_loader = torch.utils.data.DataLoader(trainset, batch_size=args.batch_size, shuffle=True, num_workers=4)
    testset = datasets.MNIST('../MNIST', download=True, train=False,
        transform=transforms.Compose([
            transforms.ToTensor(),
            transforms.Normalize((0.1307,), (0.3081,))]
        )
    )
    test_loader = torch.utils.data.DataLoader(testset, batch_size=args.test_batch_size, shuffle=True, num_workers=4)
else:
    DATA = Path("D:/datasets/imagenet-mini")
    normalize = transforms.Normalize(mean=[0.485, 0.456, 0.406],
                                     std=[0.229, 0.224, 0.225])
    transform_train = transforms.Compose([
        transforms.Resize(64),
        transforms.RandomResizedCrop(64),
        transforms.RandomHorizontalFlip(),
        transforms.ToTensor(),
        normalize,
    ])

    transform_test = transforms.Compose([
        transforms.Resize(64),
        transforms.CenterCrop(64),
        transforms.ToTensor(),
        normalize,
    ])
    train_dataset = datasets.ImageFolder(DATA / 'train', transform_train)
    test_dataset = datasets.ImageFolder(DATA / 'val', transform_test)
    train_loader = torch.utils.data.DataLoader(
        train_dataset, batch_size=args.batch_size, shuffle=True,
        num_workers=6, pin_memory=True)
    test_loader = torch.utils.data.DataLoader(
        test_dataset, batch_size=args.test_batch_size, shuffle=False,
        pin_memory=True, num_workers=6)
    # Data loading code
    #traindir = os.path.join(args.data, 'train')
    #valdir = os.path.join(args.data, 'val')
    #normalize = transforms.Normalize(mean=[0.485, 0.456, 0.406],
    #                                 std=[0.229, 0.224, 0.225])

    #train_dataset = datasets.ImageFolder(
    #    traindir,
    #    transforms.Compose([
    #        transforms.RandomResizedCrop(224),
    #       transforms.RandomHorizontalFlip(),
    #        transforms.ToTensor(),
    #        normalize,
    #    ]))

    #train_loader = torch.utils.data.DataLoader(
    #    train_dataset, batch_size=args.batch_size, shuffle=True,
    #    num_workers=16, pin_memory=True)

    #test_loader = torch.utils.data.DataLoader(
    #    datasets.ImageFolder(valdir, transforms.Compose([
    #        transforms.Resize(256),
    #        transforms.CenterCrop(224),
    #        transforms.ToTensor(),
    #        normalize,
    #    ])),
    #    batch_size=args.test_batch_size, shuffle=False,
    #    num_workers=16, pin_memory=True)

if args.dataset == 'imagenet':
    num_cls = 1000
elif args.dataset == 'cifar10':
    num_cls = 10
elif args.dataset == 'cifar100':
    num_cls = 100
elif args.dataset == 'mnist':
    num_cls = 10
else:
    raise NotImplementedError('No such dataset!')

model = models.__dict__['resnet50_add'](num_classes=num_cls, quantize=args.add_quant, weight_bits=args.add_bits)
if args.reinit is not None:
    model_ref = models.__dict__[args.arch](num_classes=num_cls, quantize=args.add_quant, weight_bits=args.add_bits)

if 'shift' in args.arch: # no pretrain
    model, _ = convert_to_shift(model, args.shift_depth, args.shift_type, convert_weights=False, use_kernel=args.use_kernel, rounding=args.rounding, weight_bits=args.weight_bits)
    if args.reinit is not None:
        model_ref, _ = convert_to_shift(model_ref, args.shift_depth, args.shift_type, convert_weights=False, use_kernel=args.use_kernel, rounding=args.rounding, weight_bits=args.weight_bits)

# create optimizer
model_other_params = []
model_sign_params = []
model_shift_params = []

for name, param in model.named_parameters():
    if(name.endswith(".sign")):
        model_sign_params.append(param)
    elif(name.endswith(".shift")):
        model_shift_params.append(param)
    else:
        model_other_params.append(param)

params_dict = [
    {"params": model_other_params},
    {"params": model_sign_params, 'lr': args.lr_sign if args.lr_sign is not None else args.lr, 'weight_decay': 0},
    {"params": model_shift_params, 'lr': args.lr, 'weight_decay': 0}
    ]

# optimizer = optim.SGD(model.parameters(), lr=args.lr, momentum=args.momentum, weight_decay=args.weight_decay)
optimizer = None
if (args.optimizer.lower() == "sgd"):
    optimizer = torch.optim.SGD(params_dict, args.lr, momentum=args.momentum, weight_decay=args.weight_decay)
elif (args.optimizer.lower() == "adadelta"):
    optimizer = torch.optim.Adadelta(params_dict, args.lr, weight_decay=args.weight_decay)
elif (args.optimizer.lower() == "adagrad"):
    optimizer = torch.optim.Adagrad(params_dict, args.lr, weight_decay=args.weight_decay)
elif (args.optimizer.lower() == "adam"):
    optimizer = torch.optim.Adam(params_dict, args.lr, weight_decay=args.weight_decay)
elif (args.optimizer.lower() == "rmsprop"):
    optimizer = torch.optim.RMSprop(params_dict, args.lr, weight_decay=args.weight_decay)
elif (args.optimizer.lower() == "radam"):
    optimizer = optim.RAdam(params_dict, args.lr, weight_decay=args.weight_decay)
elif (args.optimizer.lower() == "ranger"):
    optimizer = optim.Ranger(params_dict, args.lr, weight_decay=args.weight_decay)
else:
    raise ValueError("Optimizer type: ", args.optimizer, " is not supported or known")

scheduler_poly_lr_decay = PolynomialLRDecay(optimizer, max_decay_steps=args.epochs, end_learning_rate=0.0001, power=0.9)

def save_checkpoint(state, is_best, epoch, filepath):
    if epoch == 'init':
        filepath = os.path.join(filepath, 'init.pth.tar')
        torch.save(state, filepath)
    else:
        # filename = os.path.join(filepath, 'ckpt'+str(epoch)+'.pth.tar')
        # torch.save(state, filename)
        filename = os.path.join(filepath, 'ckpt.pth.tar')
        torch.save(state, filename)
        if is_best:
            shutil.copyfile(filename, os.path.join(filepath, 'model_best.pth.tar'))

if args.resume:
    if os.path.isfile(args.resume):
        print("=> loading checkpoint '{}'".format(args.resume))
        checkpoint = torch.load(args.resume, map_location='cpu')
        # args.start_epoch = checkpoint['epoch']
        # best_prec1 = checkpoint['best_prec1']
        model.load_state_dict(checkpoint['state_dict'])
        # optimizer.load_state_dict(checkpoint['optimizer'])
        print("=> loaded checkpoint '{}' (epoch {})"
              .format(args.resume, checkpoint['epoch']))
    else:
        print("=> no checkpoint found at '{}'".format(args.resume))
else:
    save_checkpoint({'state_dict': model.state_dict()}, False, epoch='init', filepath=args.save)

def change_name(state_dict):
    new_state_dict = OrderedDict()
    for key, value in state_dict.items():
        if 'conv' in key and '.1.weight' in key:
            new_key = key.replace('weight', 'adder')
        elif 'downsample' in key and '.1.weight' in key:
            new_key = key.replace('weight', 'adder')
        else:
            new_key = key
        new_state_dict[new_key] = value
    return new_state_dict

if args.cuda:
    model.cuda()
    if args.reinit is not None:
        model_ref.cuda()
if len(args.gpu_ids) > 1:
    model = torch.nn.DataParallel(model, device_ids=args.gpu_ids)

if args.reinit is not None:
    print('use reinit!')
    # checkpoint = torch.load(args.reinit)
    # checkpoint['state_dict'] = change_name(checkpoint['state_dict'])
    # model_ref.load_state_dict(checkpoint['state_dict'])

    for m, m_ref in zip(model.modules(), model_ref.modules()):
        if isinstance(m, deepshift.modules_q.Conv2dShiftQ): # Q shift
            weight_copy = m.weight.data.abs().clone()
            # mask = weight_copy.gt(0).float().cuda()
            mask = torch.ne(weight_copy, 1.).float().cuda()
            m_ref.weight.data = m_ref.weight.data.mul_(mask) + 1 - mask # no shift
        elif isinstance(m, deepshift.modules.Conv2dShift): # PS shift
            weight_copy = m.shift.data.abs().clone()
            mask = weight_copy.gt(0).float().cuda()
            m_ref.shift.data.mul_(mask)
            m_ref.sign.data.mul_(mask)
        elif isinstance(m, adder_fast.Adder2D):
            adder_copy = m.adder.data.abs().clone()
            mask = adder_copy.gt(0).float().cuda()
            m_ref.adder.data.mul_(mask)
    model = model_ref

# print model summary
#model_summary = None
#try:
#model_summary, model_params_info = torchsummary.summary(model, (3,32,32))
#print(model_summary)
#print("WARNING: The summary function reports duplicate parameters for multi-GPU case")
#except:
#    print("WARNING: Unable to obtain summary of model")

# save name
# name model sub-directory "shift_all" if all layers are converted to shift layers
conv2d_layers_count = count_layer_type(model, nn.Conv2d) #+ count_layer_type(model, unoptimized.UnoptimizedConv2d)
linear_layers_count = count_layer_type(model, nn.Linear) #+ count_layer_type(model, unoptimized.UnoptimizedLinear)
print('{:<30}  {:<8}'.format('conv2d counts: ', conv2d_layers_count))
print('{:<30}  {:<8}'.format('linear counts: ', linear_layers_count))

macs, params = get_model_complexity_info(model, (3, 32, 32), as_strings=True, print_per_layer_stat=True)
print('{:<30}  {:<8}'.format('Computational complexity: ', macs))
print('{:<30}  {:<8}'.format('Number of parameters: ', params))

max_stars_repo_name: JIABI/GhostShiftAddNet

max_issues_repo_name: selkerdawy/GhostShiftAddNet

max_forks_repo_name: selkerdawy/GhostShiftAddNet


text:
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats as stats

BetaAgent = None

class Environment(object):
    """
    Where agents and bandits interact.
    """
    def __init__(self, bandit, agents, label='Multi-Armed Bandit'):
        self.bandit = bandit
        self.agents = agents
        self.label = label

    def reset(self):
        self.bandit.reset()
        for agent in self.agents:
            agent.reset()

    def run(self, trials=100, experiments=1):
        scores = np.zeros((trials, len(self.agents)))
        optimal = np.zeros_like(scores)

        for _ in range(experiments):
            self.reset()
            for t in range(trials):
                for i, agent in enumerate(self.agents):
                    action = agent.choose()
                    reward, is_optimal = self.bandit.pull(action)
                    agent.observe(reward)

                    scores[t, i] += reward
                    if is_optimal:
                        optimal[t, i] += 1

        return scores / experiments, optimal / experiments

    def plot_results(self, scores, optimal):
        sns.set_style('white')
        sns.set_context('talk')
        plt.subplot(2, 1, 1)
        plt.title(self.label)
        plt.plot(scores)
        plt.ylabel('Average Reward')
        plt.legend(self.agents, loc=4)
        plt.subplot(2, 1, 2)
        plt.plot(optimal * 100)
        plt.ylim(0, 100)
        plt.ylabel('% Optimal Action')
        plt.xlabel('Time Step')
        plt.legend(self.agents, loc=4)
        sns.despine()
        plt.show()

    def plot_beliefs(self):
        sns.set_context('talk')
        pal = sns.color_palette("cubehelix", n_colors=len(self.agents))
        plt.title(self.label + ' - Agent Beliefs')

        rows = 2
        cols = int(self.bandit.k_arms / 2)

        axes = [plt.subplot(rows, cols, i+1) for i in range(self.bandit.k_arms)]
        for i, val in enumerate(self.bandit.action_values):
            color = 'r' if i == self.bandit.optimal else 'k'
            axes[i].vlines(val, 0, 1, colors=color)

        for i, agent in enumerate(self.agents):
            if type(agent) is not BetaAgent:
                for j, val in enumerate(agent.value_estimates):
                    axes[j].vlines(val, 0, 0.75, colors=pal[i], alpha=0.8)
            else:
                x = np.arange(0, 1, 0.001)
                y = np.array([stats.beta.pdf(x, a, b) for a, b in
                             zip(agent.alpha, agent.beta)])
                y /= np.max(y)
                for j, _y in enumerate(y):
                    axes[j].plot(x, _y, color=pal[i], alpha=0.8)

        min_p = np.argmin(self.bandit.action_values)
        for i, ax in enumerate(axes):
            ax.set_xlim(0, 1)
            if i % cols != 0:
                ax.set_yticklabels([])
            if i < cols:
                ax.set_xticklabels([])
            else:
                ax.set_xticks([0, 0.25, 0.5, 0.75, 1.0])
                ax.set_xticklabels(['0', '', '0.5', '', '1'])
            if i == int(cols/2):
                title = '{}-arm Bandit - Agent Estimators'.format(self.bandit.k_arms)
                ax.set_title(title)
            if i == min_p:
                ax.legend(self.agents)

        sns.despine()
        plt.show()

max_stars_repo_name: yngtodd/tanuki

max_issues_repo_name: yngtodd/tanuki

max_forks_repo_name: yngtodd/tanuki


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
#    Project: Sift implementation in Python + OpenCL
#             https://github.com/silx-kit/silx
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation
# files (the "Software"), to deal in the Software without
# restriction, including without limitation the rights to use,
# copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following
# conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
# OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
# WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
# OTHER DEALINGS IN THE SOFTWARE.

"""
Simple test of an addition
"""

from __future__ import division, print_function

__authors__ = ["Henri Payno, Jérôme Kieffer"]
__contact__ = "jerome.kieffer@esrf.eu"
__license__ = "MIT"
__copyright__ = "2013 European Synchrotron Radiation Facility, Grenoble, France"
__date__ = "09/01/2018"

import logging
import numpy

import unittest
from ..common import ocl, _measure_workgroup_size
if ocl:
    import pyopencl.array
from ..utils import get_opencl_code
from ...test.utilstest import UtilsTest
logger = logging.getLogger(__name__)

@unittest.skipIf(UtilsTest.opencl is False, "User request to skip OpenCL tests")
@unittest.skipUnless(ocl, "PyOpenCl is missing")
class TestAddition(unittest.TestCase):

    @classmethod
    def setUpClass(cls):
        super(TestAddition, cls).setUpClass()
        if ocl:
            cls.ctx = ocl.create_context()
            if logger.getEffectiveLevel() <= logging.INFO:
                cls.PROFILE = True
                properties = pyopencl.command_queue_properties.PROFILING_ENABLE
                cls.queue = pyopencl.CommandQueue(cls.ctx, properties=properties)
            else:
                cls.PROFILE = False
                cls.queue = pyopencl.CommandQueue(cls.ctx)
            cls.max_valid_wg = 0

    @classmethod
    def tearDownClass(cls):
        super(TestAddition, cls).tearDownClass()
        print("Maximum valid workgroup size %s on device %s" % (cls.max_valid_wg, cls.ctx.devices[0]))
        cls.ctx = None
        cls.queue = None

    def setUp(self):
        if ocl is None:
            return
        self.shape = 4096
        self.data = numpy.random.random(self.shape).astype(numpy.float32)
        self.d_array_img = pyopencl.array.to_device(self.queue, self.data)
        self.d_array_5 = pyopencl.array.zeros_like(self.d_array_img) - 5
        self.program = pyopencl.Program(self.ctx, get_opencl_code("addition")).build()

    def tearDown(self):
        self.img = self.data = None
        self.d_array_img = self.d_array_5 = self.program = None

    @unittest.skipUnless(ocl, "pyopencl is missing")
    def test_add(self):
        """
        tests the addition  kernel
        """
        maxi = int(round(numpy.log2(self.shape)))
        for i in range(maxi):
            d_array_result = pyopencl.array.empty_like(self.d_array_img)
            wg = 1 << i
            try:
                evt = self.program.addition(self.queue, (self.shape,), (wg,),
                                            self.d_array_img.data,
                                            self.d_array_5.data,
                                            d_array_result.data,
                                            numpy.int32(self.shape))
                evt.wait()
            except Exception as error:
                max_valid_wg = self.program.addition.get_work_group_info(pyopencl.kernel_work_group_info.WORK_GROUP_SIZE, self.ctx.devices[0])
                msg = "Error %s on WG=%s: %s" % (error, wg, max_valid_wg)
                self.assertLess(max_valid_wg, wg, msg)
                break
            else:
                res = d_array_result.get()
                good = numpy.allclose(res, self.data - 5)
                if good and wg > self.max_valid_wg:
                    self.__class__.max_valid_wg = wg
                self.assert_(good, "calculation is correct for WG=%s" % wg)

    @unittest.skipUnless(ocl, "pyopencl is missing")
    def test_measurement(self):
        """
        tests that all devices are working properly ...
        """
        for platform in ocl.platforms:
            for _did, device in enumerate(platform.devices):
                meas = _measure_workgroup_size((platform.id, device.id))
                self.assertEqual(meas, device.max_work_group_size,
                                 "Workgroup size for %s/%s: %s == %s" % (platform, device, meas, device.max_work_group_size))

def suite():
    loader = unittest.defaultTestLoader.loadTestsFromTestCase
    testSuite = unittest.TestSuite()
    testSuite.addTest(loader(TestAddition))
    return testSuite

if __name__ == '__main__':
    unittest.main(defaultTest="suite")

max_stars_repo_name: vallsv/pyFAI

max_issues_repo_name: vallsv/pyFAI

max_forks_repo_name: vallsv/pyFAI


text:
import neural_network_lyapunov.mip_utils as mip_utils

import unittest

import torch
import numpy as np
from itertools import chain, combinations
import itertools
import gurobipy

import neural_network_lyapunov.utils as utils

class TestStrengthenReLUMipWIndices(unittest.TestCase):
    def special_cases_tester(self, c):
        dtype = torch.float64
        w = torch.tensor([3., -2.], dtype=dtype)
        b = torch.tensor(2., dtype=dtype)
        lo = torch.tensor([-2., -3.], dtype=dtype)
        up = torch.tensor([2., 4.], dtype=dtype)
        # Test the case when the index set is empty.
        # This should corresponds to
        # y <= c*(wx + b) + (1-c) * max(w*x+b) * beta
        x_coeff, binary_coeff, constant = \
            mip_utils.strengthen_relu_mip_w_indices(
                c, w, b, lo, up, set())
        self.assertEqual(constant.item(), c * b.item())
        np.testing.assert_allclose(x_coeff.detach().numpy(),
                                   c * w.detach().numpy())
        relu_input_lo, relu_input_up = mip_utils.compute_range_by_IA(
            w.reshape((1, -1)), b.reshape((-1, )), lo, up)
        self.assertAlmostEqual(binary_coeff.item(),
                               (1 - c) * relu_input_up[0].item())

        # Test with index set equals to the whole set. This should corresponds
        # to y <= w*x+b + (1-c) * relu_input_lo * beta + (c-1)*relu_input_lo
        x_coeff, binary_coeff, constant = \
            mip_utils.strengthen_relu_mip_w_indices(
                c, w, b, lo, up, {0, 1})
        self.assertEqual(constant.item(), (b + (c - 1) * relu_input_lo).item())
        np.testing.assert_allclose(x_coeff.detach().numpy(),
                                   w.detach().numpy())
        np.testing.assert_allclose(binary_coeff.item(),
                                   ((1 - c) * relu_input_lo).item())

    def test_special_cases(self):
        # Test ReLU
        self.special_cases_tester(0.)
        # Test leaky ReLU
        self.special_cases_tester(0.1)

    def general_case_tester(self, c):
        dtype = torch.float64
        w = torch.tensor([3., -2.], dtype=dtype)
        b = torch.tensor(2., dtype=dtype)
        lo = torch.tensor([-2., -3.], dtype=dtype)
        up = torch.tensor([2., 4.], dtype=dtype)
        # Test index set equals to {0}
        x_coeff, binary_coeff, constant = \
            mip_utils.strengthen_relu_mip_w_indices(
                c, w, b, lo, up, {0})
        self.assertAlmostEqual(constant.item(),
                               (b * c - (1 - c) * w[0] * lo[0]).item())
        self.assertAlmostEqual(binary_coeff.item(),
                               (b * (1 - c) + (1 - c) * w[0] * lo[0] +
                                (1 - c) * w[1] * lo[1]).item())
        np.testing.assert_allclose(x_coeff.detach().numpy(),
                                   np.array([w[0], c * w[1]]))

        # Test index set equals to {1}
        x_coeff, binary_coeff, constant = \
            mip_utils.strengthen_relu_mip_w_indices(
                c, w, b, lo, up, {1})
        self.assertAlmostEqual(constant.item(),
                               (b * c - (1 - c) * w[1] * up[1]).item())
        self.assertAlmostEqual(binary_coeff.item(),
                               (b * (1 - c) + (1 - c) * w[1] * up[1] +
                                (1 - c) * w[0] * up[0]).item())
        np.testing.assert_allclose(x_coeff.detach().numpy(),
                                   np.array([c * w[0], w[1]]))

    def test_general_case(self):
        self.general_case_tester(0.)
        self.general_case_tester(0.1)

class TestFindIndexSetToStrengthen(unittest.TestCase):
    def maximal_separation_tester(self, c, w, b, lo, up, xhat, beta_hat):
        # The maximal separated plane has the maximal violation of
        # y <= bc + b(1-c)β + ∑ i∈ℑ (wᵢxᵢ−(1−c)(1−β)wᵢL̅ᵢ)
        #      + ∑i∉ℑ(cwᵢxᵢ+(1−c)βwᵢU̅ᵢ)
        # Namely the right-hand side should be minimized
        indices = mip_utils.find_index_set_to_strengthen(
            w, lo, up, xhat, beta_hat)
        min_rhs = np.inf
        max_separation_set = None
        # Loop through all subsets of {0, 1, ..., nx-1}, make sure `indices`
        # returns the maximal separation plane.
        nx = w.shape[0]
        for candidate_index in chain.from_iterable(
                combinations(list(range(nx)), r) for r in range(nx + 1)):
            x_coeff, binary_coeff, constant = \
                mip_utils.strengthen_relu_mip_w_indices(
                    c, w, b, lo, up, set(candidate_index))
            rhs = x_coeff @ xhat + binary_coeff * beta_hat + constant
            if rhs < min_rhs:
                min_rhs = rhs
                max_separation_set = set(candidate_index)
        self.assertSetEqual(indices, max_separation_set)

    def test1(self):
        dtype = torch.float64
        self.maximal_separation_tester(c=0.1,
                                       w=torch.tensor([2., -3.], dtype=dtype),
                                       b=torch.tensor(-2., dtype=dtype),
                                       lo=torch.tensor([-2., -4.],
                                                       dtype=dtype),
                                       up=torch.tensor([3., 1.], dtype=dtype),
                                       xhat=torch.tensor([-2., 1.],
                                                         dtype=dtype),
                                       beta_hat=torch.tensor(0.3, dtype=dtype))

    def test2(self):
        dtype = torch.float64
        self.maximal_separation_tester(c=0.1,
                                       w=torch.tensor([2., -3.], dtype=dtype),
                                       b=torch.tensor(-2., dtype=dtype),
                                       lo=torch.tensor([-2., -4.],
                                                       dtype=dtype),
                                       up=torch.tensor([3., 1.], dtype=dtype),
                                       xhat=torch.tensor([-2., -4.],
                                                         dtype=dtype),
                                       beta_hat=torch.tensor(0.8, dtype=dtype))

    def test3(self):
        dtype = torch.float64
        self.maximal_separation_tester(c=0.1,
                                       w=torch.tensor([2., -3.], dtype=dtype),
                                       b=torch.tensor(2., dtype=dtype),
                                       lo=torch.tensor([-2., -4.],
                                                       dtype=dtype),
                                       up=torch.tensor([3., 1.], dtype=dtype),
                                       xhat=torch.tensor([3., -4.],
                                                         dtype=dtype),
                                       beta_hat=torch.tensor(0.8, dtype=dtype))

    def test4(self):
        dtype = torch.float64
        self.maximal_separation_tester(c=0.1,
                                       w=torch.tensor([2., -3.], dtype=dtype),
                                       b=torch.tensor(2., dtype=dtype),
                                       lo=torch.tensor([-2., -4.],
                                                       dtype=dtype),
                                       up=torch.tensor([3., 1.], dtype=dtype),
                                       xhat=torch.tensor([3., 1.],
                                                         dtype=dtype),
                                       beta_hat=torch.tensor(0.8, dtype=dtype))

class TestStrengthenReLUMipGivenPts(unittest.TestCase):
    def test1(self):
        dtype = torch.float64
        c = 0.
        w = torch.tensor([1, 1], dtype=dtype)
        b = torch.tensor(-1.5, dtype=dtype)
        lo = torch.tensor([0, 0], dtype=dtype)
        up = torch.tensor([1, 1], dtype=dtype)
        linear_inputs = [
            torch.tensor([0, 1], dtype=dtype),
            torch.tensor([1, 0], dtype=dtype)
        ]
        relu_outputs = [
            torch.tensor(0.25, dtype=dtype),
            torch.tensor(0.25, dtype=dtype)
        ]
        relu_activations = [
            torch.tensor(0.5, dtype=dtype),
            torch.tensor(0.5, dtype=dtype)
        ]
        x_coeffs, binary_coeffs, constants = \
            mip_utils.strengthen_relu_mip_given_pts(
                c, w, b, lo, up, linear_inputs, relu_outputs, relu_activations)
        self.assertEqual(x_coeffs.shape, (2, 2))
        self.assertEqual(binary_coeffs.shape, (2, ))
        self.assertEqual(constants.shape, (2, ))
        for i in range(2):
            # Each additional constraint should rule out each point.
            self.assertGreater(
                relu_outputs[i], x_coeffs[i] @ linear_inputs[i] +
                binary_coeffs[i] * relu_activations[i] + constants[i])
            # The point should pass (x, y, beta) = ([0, 0], 0, 0) and
            # ([1, 1], 0.5, 1)
            self.assertAlmostEqual(constants[i].item(), 0)
            self.assertAlmostEqual(
                (x_coeffs[i] @ torch.tensor([1, 1], dtype=dtype) +
                 binary_coeffs[i] * 1 + constants[i]).item(), 0.5)
        # additional constraint 1 should pass (x, y, beta) = [0, 1], 0, 0
        self.assertAlmostEqual(
            (x_coeffs[0] @ torch.tensor([0, 1], dtype=dtype) +
             constants[0]).item(), 0)
        # additional constraint 2 should pass (x, y, beta) = [1, 0], 0, 0
        self.assertAlmostEqual(
            (x_coeffs[1] @ torch.tensor([1, 0], dtype=dtype) +
             constants[1]).item(), 0)
        for i in range(2):
            x_coeff_expected, binary_coeff_expected, constant_expected =\
                mip_utils.strengthen_relu_mip_w_indices(
                    c, w, b, lo, up,
                    mip_utils.find_index_set_to_strengthen(
                        w, lo, up, linear_inputs[i], relu_activations[i]))
            np.testing.assert_allclose(x_coeffs[i].detach().numpy(),
                                       x_coeff_expected.detach().numpy())
            self.assertEqual(binary_coeffs[i].item(),
                             binary_coeff_expected.item())
            self.assertEqual(constants[i].item(), constant_expected.item())

    def test2(self):
        # This test adds no additional constraint, since the candidate
        # constraint is not violated at the evaluation point.
        dtype = torch.float64
        c = 0.
        w = torch.tensor([1, 1], dtype=dtype)
        b = torch.tensor(-1.5, dtype=dtype)
        lo = torch.tensor([0, 0], dtype=dtype)
        up = torch.tensor([1, 1], dtype=dtype)
        linear_inputs = [
            torch.tensor([0, 1], dtype=dtype),
            torch.tensor([1, 0], dtype=dtype)
        ]
        relu_outputs = [
            torch.tensor(0., dtype=dtype),
            torch.tensor(0., dtype=dtype)
        ]
        relu_activations = [
            torch.tensor(0., dtype=dtype),
            torch.tensor(0., dtype=dtype)
        ]
        x_coeffs, binary_coeffs, constants = \
            mip_utils.strengthen_relu_mip_given_pts(
                c, w, b, lo, up, linear_inputs, relu_outputs, relu_activations)
        self.assertIsNone(x_coeffs)
        self.assertIsNone(binary_coeffs)
        self.assertIsNone(constants)

class TestComputeBetaRange(unittest.TestCase):
    def empty_constraint_tester(self, c):
        # Test when x_coeffs = beta_coeffs = constants = None
        dtype = torch.float64
        w = torch.tensor([2., -3.], dtype=dtype)
        b = torch.tensor([-.5], dtype=dtype)
        beta_range = mip_utils._compute_beta_range(
            c, w, b, [], [], [], torch.tensor([-2, 3], dtype=dtype))
        self.assertEqual(beta_range[0].item(), 0.)
        self.assertEqual(beta_range[1].item(), 1.)

    def test_empty_constraint(self):
        self.empty_constraint_tester(0.)
        self.empty_constraint_tester(0.1)

    def nonempty_constraint_tester(self, c, w, b, lo, up, index_sets, xhat):
        relu_input_lo, relu_input_up = mip_utils.compute_range_by_IA(
            w.reshape((1, -1)), b.reshape((-1, )), lo, up)
        assert (relu_input_lo[0].item() < 0)
        assert (relu_input_up[0].item() > 0)
        x_coeffs = []
        beta_coeffs = []
        constants = []
        for index_set in index_sets:
            x_coeff, beta_coeff, constant =\
                mip_utils.strengthen_relu_mip_w_indices(
                    c, w, b, lo, up, index_set)
            x_coeffs.append(x_coeff)
            beta_coeffs.append(beta_coeff)
            constants.append(constant)

        beta_lo, beta_up = mip_utils._compute_beta_range(
            c, w, b, x_coeffs, beta_coeffs, constants, xhat)
        self.assertGreaterEqual(beta_lo.item(), 0)
        self.assertLessEqual(beta_up.item(), 1)
        self.assertGreaterEqual(beta_up.item(), beta_lo.item())
        x_coeffs_torch = torch.cat([v.reshape((1, -1)) for v in x_coeffs],
                                   dim=0)
        beta_coeffs_torch = torch.stack(beta_coeffs)
        constants_torch = torch.stack(constants)
        lhs = torch.max(c * (w @ xhat + b), w @ xhat + b)
        np.testing.assert_array_less(lhs.item() - 1E-6,
                                     (x_coeffs_torch @ xhat +
                                      beta_coeffs_torch * beta_up +
                                      constants_torch).detach().numpy())
        np.testing.assert_array_less(lhs.item() - 1E-6,
                                     (x_coeffs_torch @ xhat +
                                      beta_coeffs_torch * beta_lo +
                                      constants_torch).detach().numpy())
        np.testing.assert_array_less(
            lhs.item() - 1E-6,
            (x_coeffs_torch @ xhat + beta_coeffs_torch *
             (beta_lo + beta_up) / 2 + constants_torch).detach().numpy())
        # Make sure that if beta > beta_up, then it violates the constraints.
        if beta_up.item() < 1:
            self.assertFalse(
                torch.all(x_coeffs_torch @ xhat + beta_coeffs_torch *
                          (beta_up + 1E-5) + constants_torch >= lhs.item()))
        if beta_lo.item() > 0:
            self.assertFalse(
                torch.all(x_coeffs_torch @ xhat + beta_coeffs_torch *
                          (beta_lo - 1E-5) + constants_torch >= lhs.item()))

    def test_nonempty_constraint1(self):
        dtype = torch.float64
        w = torch.tensor([2., -3.], dtype=dtype)
        b = torch.tensor(-0.5, dtype=dtype)
        lo = torch.tensor([-1., -2.], dtype=dtype)
        up = torch.tensor([-0.5, 1.], dtype=dtype)
        for c in (0., 0.1):
            self.nonempty_constraint_tester(c, w, b, lo, up, [set()], lo)
            self.nonempty_constraint_tester(c, w, b, lo, up, [set()], up)
            self.nonempty_constraint_tester(
                c, w, b, lo, up, [set(), {0}],
                torch.tensor([lo[0], up[1]], dtype=dtype))
            self.nonempty_constraint_tester(
                c, w, b, lo, up, [set(), {0, 1}],
                torch.tensor([lo[0], up[1]], dtype=dtype))

    def test_nonempty_constraint2(self):
        dtype = torch.float64
        w = torch.tensor([2., -3., 1.5], dtype=dtype)
        b = torch.tensor(-0.5, dtype=dtype)
        lo = torch.tensor([-1., -2., 0.2], dtype=dtype)
        up = torch.tensor([-0.5, 1., 0.5], dtype=dtype)
        for c in (0., 0.1):
            self.nonempty_constraint_tester(c, w, b, lo, up, [set()], lo)
            self.nonempty_constraint_tester(c, w, b, lo, up, [set()], up)
            self.nonempty_constraint_tester(
                c, w, b, lo, up, [set(), {0}],
                torch.tensor([lo[0], up[1], lo[2]], dtype=dtype))
            self.nonempty_constraint_tester(
                c, w, b, lo, up, [set(), {0, 1}],
                torch.tensor([lo[0], up[1], up[2]], dtype=dtype))

class TestGetLinearInputVertices(unittest.TestCase):
    def test1(self):
        # Test with relu_input_lo, relu_input_up computed from interval
        # arithmetics.
        dtype = torch.float64
        lo = torch.tensor([-1, -3], dtype=dtype)
        up = torch.tensor([-0.5, 2], dtype=dtype)
        w = torch.tensor([-2, 3.], dtype=dtype)
        b = torch.tensor(0.4, dtype=dtype)
        relu_input_lo_ia, relu_input_up_ia = mip_utils.compute_range_by_IA(
            w.reshape((1, -1)), b.reshape((-1, )), lo, up)
        vertices = mip_utils._get_linear_input_vertices(
            lo, up, w, b, relu_input_lo_ia[0], relu_input_up_ia[0])
        box_vertices = itertools.product(*[(lo[i], up[i]) for i in range(2)])
        self.assertEqual(len(vertices), 4)
        i = 0
        for box_vertex in box_vertices:
            vertex_expect = torch.stack(box_vertex)
            np.testing.assert_allclose(vertices[i].detach().numpy(),
                                       vertex_expect.detach().numpy())
            i += 1

    def test2(self):
        # Test with relu_input_lo, relu_input_up tighter than that computed
        # from interval arithmetics.
        dtype = torch.float64
        lo = torch.tensor([-1, -3], dtype=dtype)
        up = torch.tensor([-0.5, 2], dtype=dtype)
        w = torch.tensor([-2, 3.], dtype=dtype)
        b = torch.tensor(0.4, dtype=dtype)
        relu_input_lo_ia, relu_input_up_ia = mip_utils.compute_range_by_IA(
            w.reshape((1, -1)), b.reshape((-1, )), lo, up)
        vertices = mip_utils._get_linear_input_vertices(
            lo, up, w, b, relu_input_lo_ia[0], relu_input_lo_ia[0] + 1E-5)
        self.assertEqual(len(vertices), 1)
        np.testing.assert_allclose(vertices[0], np.array([-0.5, -3]))

class TestMaxYGivenLinearInput(unittest.TestCase):
    def y_range_tester(self, c, w, b, x_coeffs, beta_coeffs, constants, x_hat):
        y_max, beta_val = mip_utils._max_y_given_linear_input(
            c, w, b, x_coeffs, beta_coeffs, constants, x_hat)
        # First make sure that this y_max and beta_val satisfy the constraint
        lhs = torch.max(c * (w @ x_hat + b), w @ x_hat + b).item()
        self.assertGreaterEqual(y_max, lhs - 1E-6)
        self.assertGreaterEqual(beta_val, 0)
        self.assertLessEqual(beta_val, 1)
        for x_coeff, beta_coeff, constant in zip(x_coeffs, beta_coeffs,
                                                 constants):
            self.assertLessEqual(
                y_max,
                (x_coeff @ x_hat + beta_coeff * beta_val + constant).item() +
                1E-6)
        x_coeffs_stack = torch.cat([v.reshape((1, -1)) for v in x_coeffs],
                                   dim=0)
        beta_coeffs_stack = torch.stack(beta_coeffs)
        constants_stack = torch.stack(constants)
        # Now sample many beta in the range [0, 1]. Either they violate the
        # constraint, or the corresponding y is smaller than y_max.
        for beta in np.linspace(0, 1, 51).tolist():
            rhs = x_coeffs_stack @ x_hat + beta_coeffs_stack * beta +\
                constants_stack
            self.assertTrue(
                torch.any(rhs < lhs) or torch.min(rhs).item() <= y_max)

    def test1(self):
        # Test when relu_input_lo, relu_input_up is computed from interval
        # arithmetic of linear_input_lo and linear_input_up
        dtype = torch.float64
        c = 0.
        w = torch.tensor([1, 1], dtype=dtype)
        b = torch.tensor(-1.5, dtype=dtype)
        lo = torch.tensor([0, 0], dtype=dtype)
        up = torch.tensor([1, 1], dtype=dtype)
        relu_input_lo, relu_input_up = mip_utils.compute_range_by_IA(
            w.reshape((1, -1)), b.reshape((-1, )), lo, up)
        x_coeffs = []
        beta_coeffs = []
        constants = []
        for index_set in (set(), {0, 1}):
            x_coeff, beta_coeff, constant = \
                mip_utils.strengthen_relu_mip_w_indices(
                    c, w, b, lo, up, index_set)
            x_coeffs.append(x_coeff)
            beta_coeffs.append(beta_coeff)
            constants.append(constant)

        for x_hat in (torch.tensor([0, 0], dtype=dtype),
                      torch.tensor([0, 1], dtype=dtype),
                      torch.tensor([1, 0], dtype=dtype),
                      torch.tensor([1, 1], dtype=dtype)):
            self.y_range_tester(c, w, b, x_coeffs, beta_coeffs, constants,
                                x_hat)
        x_coeff, beta_coeff, constant = \
            mip_utils.strengthen_relu_mip_w_indices(
                c, w, b, lo, up, {0})
        x_coeffs.append(x_coeff)
        beta_coeffs.append(beta_coeff)
        constants.append(constant)
        for x_hat in (torch.tensor([0, 0], dtype=dtype),
                      torch.tensor([0, 1], dtype=dtype),
                      torch.tensor([1, 0], dtype=dtype),
                      torch.tensor([1, 1], dtype=dtype)):
            self.y_range_tester(c, w, b, x_coeffs, beta_coeffs, constants,
                                x_hat)

    def test2(self):
        # Test when relu_input_lo, relu_input_up is tighter than those
        # computed from interval arithmetic of linear_input_lo and
        # linear_input_up
        dtype = torch.float64
        c = 0.1
        w = torch.tensor([2, -1], dtype=dtype)
        b = torch.tensor(-1.5, dtype=dtype)
        lo = torch.tensor([-1, 0.5], dtype=dtype)
        up = torch.tensor([1.5, 1], dtype=dtype)
        relu_input_lo_ia, relu_input_up_ia = mip_utils.compute_range_by_IA(
            w.reshape((1, -1)), b.reshape((-1, )), lo, up)
        relu_input_lo = 0.95 * relu_input_lo_ia[0] + 0.05 * relu_input_up_ia[0]
        relu_input_up = 0.1 * relu_input_lo_ia[0] + 0.9 * relu_input_up_ia[0]
        x_coeffs = [c * w, w]
        beta_coeffs = [(1 - c) * relu_input_up, (1 - c) * relu_input_lo]
        constants = [c * b, b - (1 - c) * relu_input_lo]

        for x_hat in mip_utils._get_linear_input_vertices(
                lo, up, w, b, relu_input_lo, relu_input_up):
            if relu_input_lo <= w @ x_hat + b <= relu_input_up:
                self.y_range_tester(c, w, b, x_coeffs, beta_coeffs, constants,
                                    x_hat)
        x_coeff, beta_coeff, constant = \
            mip_utils.strengthen_relu_mip_w_indices(
                c, w, b, lo, up, {0})
        x_coeffs.append(x_coeff)
        beta_coeffs.append(beta_coeff)
        constants.append(constant)
        for x_hat in mip_utils._get_linear_input_vertices(
                lo, up, w, b, relu_input_lo, relu_input_up):
            if relu_input_lo <= w @ x_hat + b <= relu_input_up:
                self.y_range_tester(c, w, b, x_coeffs, beta_coeffs, constants,
                                    x_hat)

class TestStrengthenReLUMip(unittest.TestCase):
    def construct_strengthened_mip(self, c, w, b, lo, up, relu_input_lo,
                                   relu_input_up, use_binary, selective):
        x_coeffs, binary_coeffs, constants = mip_utils.strengthen_relu_mip(
            c, w, b, lo, up, relu_input_lo, relu_input_up, selective)
        # Now setup an optimization problem with these constraints
        prog = gurobipy.Model()
        x = []
        for i in range(w.shape[0]):
            x.append(
                prog.addVar(lb=lo[i].item(), ub=up[i].item(), name=f"x_{i}"))
        beta_type = gurobipy.GRB.BINARY if use_binary else\
            gurobipy.GRB.CONTINUOUS
        beta = prog.addVar(lb=0., ub=1., name="b", vtype=beta_type)
        y = prog.addVar(lb=-gurobipy.GRB.INFINITY, name="y")
        # Add the lower bound constraint y >= c*(w'x+b), y >= w'x+b
        prog.addLConstr(gurobipy.LinExpr([1] + (-w).detach().numpy().tolist(),
                                         [y] + x),
                        rhs=b.item(),
                        sense=gurobipy.GRB.GREATER_EQUAL)
        prog.addLConstr(gurobipy.LinExpr(
            [1] + (-c * w).detach().numpy().tolist(), [y] + x),
                        rhs=c * b.item(),
                        sense=gurobipy.GRB.GREATER_EQUAL)
        # Now add the constraint
        # y <= c(wᵀx+b)+(1−c)m⁺β
        # y <= wᵀx+b−(1−c)m⁻(1−β)
        prog.addLConstr(gurobipy.LinExpr(
            [1] + (-c * w).detach().numpy().tolist() +
            [-(1 - c) * relu_input_up.item()], [y] + x + [beta]),
                        sense=gurobipy.GRB.LESS_EQUAL,
                        rhs=c * b.item())
        prog.addLConstr(gurobipy.LinExpr([1] + (-w).detach().numpy().tolist() +
                                         [-(1 - c) * relu_input_lo.item()],
                                         [y] + x + [beta]),
                        sense=gurobipy.GRB.LESS_EQUAL,
                        rhs=(b - (1 - c) * relu_input_lo).item())
        # Now add the strengthened constraints.
        # y <= x_coeffs * x + beta_coeffs * beta + constants
        for i in range(x_coeffs.shape[0]):
            prog.addLConstr(gurobipy.LinExpr(
                x_coeffs[i].detach().numpy().tolist() +
                [binary_coeffs[i].item()] + [-1], x + [beta] + [y]),
                            sense=gurobipy.GRB.GREATER_EQUAL,
                            rhs=-constants[i].item())
        return prog, x, y, beta

    def strengthened_constraint_tester(self, c, w, b, lo, up, relu_input_lo,
                                       relu_input_up, ideal, selective):
        # Make sure that after strengthening the constraint, the condition
        # y = relu(w*x+b) is satisfied. If we expect to get an ideal
        # formualtion, then also check that the vertices of the constrained
        # polytopes are integral solutions.

        # Now sample many different vector obj_coeff, and maximize
        # obj_coeff * [x;y]. We check if the optimal solution is integral.
        prog, x, y, beta = self.construct_strengthened_mip(c,
                                                           w,
                                                           b,
                                                           lo,
                                                           up,
                                                           relu_input_lo,
                                                           relu_input_up,
                                                           use_binary=False,
                                                           selective=selective)
        if ideal:
            for _ in range(40):
                obj_coeff = np.random.uniform(low=-1.,
                                              high=1.,
                                              size=(w.shape[0] + 1, ))
                prog.setObjective(gurobipy.LinExpr(obj_coeff.tolist(),
                                                   x + [y]),
                                  sense=gurobipy.GRB.MAXIMIZE)
                prog.setParam(gurobipy.GRB.Param.OutputFlag, False)
                prog.optimize()
                self.assertEqual(prog.status, gurobipy.GRB.Status.OPTIMAL)
                self.assertTrue(
                    np.abs(beta.x) < 1E-6 or np.abs(beta.x - 1) < 1E-6)

        # Now sample many points satisfying
        # lo <= x <= up
        # relu_input_lo <= w'*x+b <= relu_input_up
        # And check if (x, y) satisfies the strengthened constraint.
        prog, x, y, beta = self.construct_strengthened_mip(c,
                                                           w,
                                                           b,
                                                           lo,
                                                           up,
                                                           relu_input_lo,
                                                           relu_input_up,
                                                           use_binary=False,
                                                           selective=selective)
        x_samples = utils.uniform_sample_in_box(lo, up, 1000)
        for i in range(x_samples.shape[0]):
            if relu_input_lo <= w @ x_samples[i] + b <= relu_input_up:
                for j in range(w.shape[0]):
                    x[j].lb = x_samples[i, j].item()
                    prog.update()
                    x[j].ub = x_samples[i, j].item()
                    prog.update()
                y_val = torch.max(c * (w @ x_samples[i] + b),
                                  w @ x_samples[i] + b)
                y.lb = y_val.item()
                y.ub = y_val.item()
                prog.setParam(gurobipy.GRB.Param.OutputFlag, False)
                prog.optimize()
                self.assertEqual(prog.status, gurobipy.GRB.Status.OPTIMAL)

        # Now sample many points satisfying
        # lo <= x <= up
        # relu_input_lo <= w'*x+b <= relu_input_up
        # Make sure that if beta is a binary variable, then the solution
        # (y, beta) satisfies the relu condition.
        prog, x, y, beta = self.construct_strengthened_mip(c,
                                                           w,
                                                           b,
                                                           lo,
                                                           up,
                                                           relu_input_lo,
                                                           relu_input_up,
                                                           use_binary=True,
                                                           selective=selective)
        x_samples = utils.uniform_sample_in_box(lo, up, 1000)
        for i in range(x_samples.shape[0]):
            if relu_input_lo <= w @ x_samples[i] + b <= relu_input_up:
                for j in range(w.shape[0]):
                    x[j].lb = x_samples[i, j].item()
                    prog.update()
                    x[j].ub = x_samples[i, j].item()
                    prog.update()
                prog.setParam(gurobipy.GRB.Param.OutputFlag, False)
                prog.optimize()
                self.assertEqual(prog.status, gurobipy.GRB.Status.OPTIMAL)
                self.assertAlmostEqual(y.x,
                                       torch.max(c * (w @ x_samples[i] + b),
                                                 w @ x_samples[i] + b).item(),
                                       places=7)

    def test_2d_case1(self):
        # The example given in Strong mixed-integer programming formulations
        # for trained neural networks.
        c = 0.
        dtype = torch.float64
        w = torch.tensor([1, 1], dtype=dtype)
        b = torch.tensor(-1.5, dtype=dtype)
        lo = torch.tensor([0, 0], dtype=dtype)
        up = torch.tensor([1, 1], dtype=dtype)
        relu_input_lo_ia, relu_input_up_ia = mip_utils.compute_range_by_IA(
            w.reshape((1, -1)), b.reshape((-1, )), lo, up)
        relu_input_lo = relu_input_lo_ia[0]
        relu_input_up = relu_input_up_ia[0]
        for selective in (True, False):
            x_coeffs, binary_coeffs, constants = mip_utils.strengthen_relu_mip(
                c, w, b, lo, up, relu_input_lo, relu_input_up, selective)
            self.assertEqual(x_coeffs.shape, (2, 2))
            self.strengthened_constraint_tester(c,
                                                w,
                                                b,
                                                lo,
                                                up,
                                                relu_input_lo,
                                                relu_input_up,
                                                ideal=True,
                                                selective=selective)

        # relu_input_lo and relu_input_up are tighter than those computed from
        # interval arithematics.
        relu_input_lo = torch.tensor(-0.6, dtype=dtype)
        relu_input_up = relu_input_up_ia[0]
        for selective in (True, False):
            x_coeffs, binary_coeffs, constants = mip_utils.strengthen_relu_mip(
                c, w, b, lo, up, relu_input_lo, relu_input_up, selective)
            self.assertEqual(x_coeffs.shape, (2, 2))
            self.strengthened_constraint_tester(c,
                                                w,
                                                b,
                                                lo,
                                                up,
                                                relu_input_lo,
                                                relu_input_up,
                                                ideal=False,
                                                selective=selective)

        # Now tighten relu_input_lo further, such that the vertices
        # [x, y] = [0, 1, 0], [1, 0, 0] are not in the admissible region.
        relu_input_lo = torch.tensor(-0.45, dtype=dtype)
        for (selective, size) in ((True, 0), (False, 2)):
            x_coeffs, binary_coeffs, constants = mip_utils.strengthen_relu_mip(
                c, w, b, lo, up, relu_input_lo, relu_input_up, selective)
            self.assertEqual(x_coeffs.shape, (size, 2))
            self.assertEqual(binary_coeffs.shape, (size, ))
            self.assertEqual(constants.shape, (size, ))
            self.strengthened_constraint_tester(c,
                                                w,
                                                b,
                                                lo,
                                                up,
                                                relu_input_lo,
                                                relu_input_up,
                                                ideal=True,
                                                selective=selective)

    def test_2d_case2(self):
        # Test with leaky ReLU.
        c = 0.1
        dtype = torch.float64
        w = torch.tensor([2, -1], dtype=dtype)
        b = torch.tensor(0.5, dtype=dtype)
        lo = torch.tensor([-1, -1.5], dtype=dtype)
        up = torch.tensor([-0.5, 2], dtype=dtype)
        relu_input_lo_ia, relu_input_up_ia = mip_utils.compute_range_by_IA(
            w.reshape((1, -1)), b.reshape((-1, )), lo, up)
        assert (relu_input_lo_ia[0].item() < 0)
        assert (relu_input_up_ia[0].item() > 0)
        relu_input_lo = relu_input_lo_ia[0]
        relu_input_up = relu_input_up_ia[0]
        self.strengthened_constraint_tester(c,
                                            w,
                                            b,
                                            lo,
                                            up,
                                            relu_input_lo,
                                            relu_input_up,
                                            ideal=True,
                                            selective=False)
        self.strengthened_constraint_tester(c,
                                            w,
                                            b,
                                            lo,
                                            up,
                                            relu_input_lo,
                                            relu_input_up,
                                            ideal=False,
                                            selective=True)

    def test_3d_case(self):
        c = 0.1
        dtype = torch.float64
        w = torch.tensor([2, -1, 0.5], dtype=dtype)
        b = torch.tensor(-1.5, dtype=dtype)
        lo = torch.tensor([-1, 0.5, -1.5], dtype=dtype)
        up = torch.tensor([2, 1.5, 1], dtype=dtype)
        relu_input_lo_ia, relu_input_up_ia = mip_utils.compute_range_by_IA(
            w.reshape((1, -1)), b.reshape((-1, )), lo, up)
        assert (relu_input_lo_ia[0].item() < 0)
        assert (relu_input_up_ia[0].item() > 0)
        relu_input_lo = relu_input_lo_ia[0]
        relu_input_up = relu_input_up_ia[0]
        self.strengthened_constraint_tester(c,
                                            w,
                                            b,
                                            lo,
                                            up,
                                            relu_input_lo,
                                            relu_input_up,
                                            ideal=True,
                                            selective=False)
        self.strengthened_constraint_tester(c,
                                            w,
                                            b,
                                            lo,
                                            up,
                                            relu_input_lo,
                                            relu_input_up,
                                            ideal=False,
                                            selective=True)

        # Test a tighter relu input bound.
        relu_input_lo = 0.8 * relu_input_lo_ia[0] + 0.2 * relu_input_up_ia[0]
        relu_input_up = 0.1 * relu_input_lo_ia[0] + 0.9 * relu_input_up_ia[0]
        self.strengthened_constraint_tester(c,
                                            w,
                                            b,
                                            lo,
                                            up,
                                            relu_input_lo,
                                            relu_input_up,
                                            ideal=True,
                                            selective=False)
        self.strengthened_constraint_tester(c,
                                            w,
                                            b,
                                            lo,
                                            up,
                                            relu_input_lo,
                                            relu_input_up,
                                            ideal=False,
                                            selective=True)

class TestComputeRangeByLP(unittest.TestCase):
    def test_x_bounds(self):
        # Test with only bounds on x.
        x_lb = np.array([1., 2.])
        x_ub = np.array([3., 4.])
        A = np.array([[0., 1.], [1., 1.], [2., -1.]])
        b = np.array([0., 1., 3.])
        y_lb, y_ub = mip_utils.compute_range_by_lp(A, b, x_lb, x_ub, None,
                                                   None)
        y_lb_expected = np.array([2., 4, 1.])
        y_ub_expected = np.array([4, 8., 7.])
        np.testing.assert_allclose(y_lb, y_lb_expected, atol=1E-7)
        np.testing.assert_allclose(y_ub, y_ub_expected, atol=1E-7)

        # Some of the bounds should be infinity.
        x_lb = np.array([1., -np.inf])
        y_lb, y_ub = mip_utils.compute_range_by_lp(A, b, x_lb, x_ub, None,
                                                   None)
        y_lb_expected = np.array([-np.inf, -np.inf, 1.])
        y_ub_expected = np.array([4, 8., np.inf])
        np.testing.assert_allclose(y_lb, y_lb_expected, atol=1E-7)
        np.testing.assert_allclose(y_ub, y_ub_expected, atol=1E-7)

        # Now the problem is infeasible as x_lb > x_ub
        x_lb = np.array([4., 2])
        y_lb, y_ub = mip_utils.compute_range_by_lp(A, b, x_lb, x_ub, None,
                                                   None)
        np.testing.assert_allclose(y_lb, np.full((3, ), np.inf))
        np.testing.assert_allclose(y_ub, np.full((3, ), -np.inf))

    def test_ineq_bounds(self):
        # x has both bounds x_lb <= x <= x_ub and inequality bounds C * x <= d
        x_lb = np.array([-1., -2.])
        x_ub = np.array([3., 4.])
        C = np.array([[1., 1.], [1., -1], [-1, 1], [-1, -1.]])
        d = np.array([2, 2, 2, 2])
        A = np.array([[0., 1.], [1., 1.], [2., -1.]])
        b = np.array([0., 1., 3.])
        y_lb, y_ub = mip_utils.compute_range_by_lp(A, b, x_lb, x_ub, C, d)
        y_lb_expected = np.array([-2, -1, 0.])
        y_ub_expected = np.array([2., 3, 7.])
        np.testing.assert_allclose(y_lb, y_lb_expected, atol=1E-7)
        np.testing.assert_allclose(y_ub, y_ub_expected, atol=1E-7)

class TestComputeRangeByIA(unittest.TestCase):
    def test(self):
        dtype = torch.float64
        A = torch.tensor([[1., 2., -3.], [2., -1., -4.]], dtype=dtype)
        b = torch.tensor([2., 3.], dtype=dtype)
        x_lb = torch.tensor([-2., 3., -4.], dtype=dtype)
        x_ub = torch.tensor([1., 5., 7.], dtype=dtype)
        output_lb, output_ub = mip_utils.compute_range_by_IA(A, b, x_lb, x_ub)
        np.testing.assert_allclose(output_lb.detach().numpy(),
                                   np.array([-15, -34]))
        np.testing.assert_allclose(output_ub.detach().numpy(),
                                   np.array([25, 18]))

    def test_gradient(self):
        def test_fun(A_np, b_np, x_lb_np, x_ub_np):
            output_lb, output_ub = mip_utils.compute_range_by_IA(
                torch.from_numpy(A_np.reshape((2, 3))), torch.from_numpy(b_np),
                torch.from_numpy(x_lb_np), torch.from_numpy(x_ub_np))
            return torch.sum(output_lb * output_ub).item()

        dtype = torch.float64
        A = torch.tensor([[1., 2., -3.], [2., -1., -4.]],
                         dtype=dtype,
                         requires_grad=True)
        b = torch.tensor([2., 3.], dtype=dtype, requires_grad=True)
        x_lb = torch.tensor([-2., 3., -4.], dtype=dtype, requires_grad=True)
        x_ub = torch.tensor([1., 5., 7.], dtype=dtype, requires_grad=True)
        output_lb, output_ub = mip_utils.compute_range_by_IA(A, b, x_lb, x_ub)
        torch.sum(output_lb * output_ub).backward()
        A_grad = A.grad.detach().numpy()
        b_grad = b.grad.detach().numpy()
        x_lb_grad = x_lb.grad.detach().numpy()
        x_ub_grad = x_ub.grad.detach().numpy()

        A_grad_numerical, b_grad_numerical, x_lb_grad_numerical,\
            x_ub_grad_numerical = utils.compute_numerical_gradient(
                test_fun, A.detach().numpy().reshape((-1)), b.detach().numpy(),
                x_lb.detach().numpy(), x_ub.detach().numpy())
        np.testing.assert_allclose(A_grad,
                                   A_grad_numerical.reshape((2, 3)),
                                   atol=1E-6)
        np.testing.assert_allclose(b_grad, b_grad_numerical, atol=1E-6)
        np.testing.assert_allclose(x_lb_grad, x_lb_grad_numerical, atol=1E-6)
        np.testing.assert_allclose(x_ub_grad, x_ub_grad_numerical, atol=1E-6)

class TestPropagateBounds(unittest.TestCase):
    def test_relu(self):
        layer = torch.nn.ReLU()
        input_lo = torch.tensor([-2., 2., -3.], dtype=torch.float64)
        input_up = torch.tensor([3., 6., -1.], dtype=torch.float64)
        output_lo, output_up = mip_utils.propagate_bounds(
            layer, input_lo, input_up)
        np.testing.assert_allclose(output_lo, np.array([0, 2., 0.]))
        np.testing.assert_allclose(output_up, np.array([3., 6., 0.]))

    def test_leaky_relu1(self):
        # negative_slope > 0
        layer = torch.nn.LeakyReLU(0.1)
        input_lo = torch.tensor([-2., 2., -3.], dtype=torch.float64)
        input_up = torch.tensor([3., 6., -1.], dtype=torch.float64)
        output_lo, output_up = mip_utils.propagate_bounds(
            layer, input_lo, input_up)
        np.testing.assert_allclose(output_lo, np.array([-0.2, 2., -0.3]))
        np.testing.assert_allclose(output_up, np.array([3., 6., -0.1]))

    def test_leaky_relu2(self):
        # negative_slope < 0
        layer = torch.nn.LeakyReLU(-0.1)
        input_lo = torch.tensor([-10, 2, -8, -1], dtype=torch.float64)
        input_up = torch.tensor([0.5, 8, -2, 4], dtype=torch.float64)
        output_lo, output_up = mip_utils.propagate_bounds(
            layer, input_lo, input_up)
        np.testing.assert_allclose(output_lo, np.array([0, 2, 0.2, 0]))
        np.testing.assert_allclose(output_up, np.array([1., 8., 0.8, 4]))

    def test_linear_layer_no_bias(self):
        layer = torch.nn.Linear(3, 2, bias=False)
        layer.weight.data = torch.tensor([[-1, 0, 2], [3, -2, -1]],
                                         dtype=torch.float64)
        input_lo = torch.tensor([-2, 2, -3], dtype=torch.float64)
        input_up = torch.tensor([3, 6, -1], dtype=torch.float64)
        output_lo, output_up = mip_utils.propagate_bounds(
            layer, input_lo, input_up)
        np.testing.assert_allclose(output_lo.detach().numpy(),
                                   np.array([-9., -17]))
        np.testing.assert_allclose(output_up.detach().numpy(),
                                   np.array([0., 8.]))
        # Now check we can take the gradient.
        loss = output_lo.sum() + output_up.sum()
        loss.backward()

        for s in np.linspace(0, 1, 11):
            output = layer(s * input_lo + (1 - s) * input_up)
            np.testing.assert_array_less(output.detach().numpy(),
                                         output_up.detach().numpy() + 1E-10)
            np.testing.assert_array_less(output_lo.detach().numpy() - 1E-10,
                                         output.detach().numpy())

    def test_linear_layer_bias(self):
        layer = torch.nn.Linear(3, 2, bias=True)
        layer.weight.data = torch.tensor([[-1, 0, 2], [3, -2, -1]],
                                         dtype=torch.float64)
        layer.bias.data = torch.tensor([2., -1.], dtype=torch.float64)
        input_lo = torch.tensor([-2, 2, -3], dtype=torch.float64)
        input_up = torch.tensor([3, 6, -1], dtype=torch.float64)
        output_lo, output_up = mip_utils.propagate_bounds(
            layer, input_lo, input_up)
        np.testing.assert_allclose(output_lo.detach().numpy(),
                                   np.array([-7., -18.]))
        np.testing.assert_allclose(output_up.detach().numpy(),
                                   np.array([2., 7.]))
        # Now check we can take the gradient.
        loss = output_lo.sum() + output_up.sum()
        loss.backward()

        for s in np.linspace(0, 1, 11):
            output = layer(s * input_lo + (1 - s) * input_up)
            np.testing.assert_array_less(output.detach().numpy(),
                                         output_up.detach().numpy() + 1E-10)
            np.testing.assert_array_less(output_lo.detach().numpy() - 1E-10,
                                         output.detach().numpy())

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

max_issues_repo_name: StanfordASL/neural-network-lyapunov

max_forks_repo_name: StanfordASL/neural-network-lyapunov


text:
import cv2
import numpy as np

def isset(v):
    try:
        type (eval(v))
    except:
        return 0
    else:
        return 1

# create video capture
cap = cv2.VideoCapture(0)

while(1):

    # Read the frames frome a camera
    _,frame = cap.read()
    frame = cv2.blur(frame,(3,3))

    # Or get it from a JPEG
    #frame = cv2.imread('frame0010.jpg', 1)

    # Convert the image to hsv space and find range of colors
    hsv = cv2.cvtColor(frame, cv2.COLOR_BGR2HSV)

    # This is for RED tag
    #thresh = cv2.inRange(hsv,np.array((120, 80, 80)), np.array((180, 255, 255)))

    # This is for GREEN tag
    thresh = cv2.inRange(hsv,np.array((50, 80, 80)), np.array((120, 255, 255)))

    thresh2 = thresh.copy()

    # find contours in the threshold image
    contours,hierarchy = cv2.findContours(thresh, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)

    # finding contour with maximum area and store it as best_cnt
    max_area = 0
    for cnt in contours:
        area = cv2.contourArea(cnt)
        if area > max_area:
            max_area = area
            best_cnt = cnt

    # finding centroids of best_cnt and draw a circle there
    if isset('best_cnt'):
        M = cv2.moments(best_cnt)
        cx,cy = int(M['m10']/M['m00']), int(M['m01']/M['m00'])
        cv2.circle(frame,(cx,cy),5,255,-1)
        print "Central pos: (%d, %d)" % (cx,cy)
    else:
        print "[Warning]Tag lost..."

    # Show the original and processed image
    cv2.imshow('frame', frame)
    cv2.imshow('thresh', thresh2)

    # if key pressed is 'Esc' then exit the loop
    if cv2.waitKey(33)== 27:
        break

# Clean up and exit the program
cv2.destroyAllWindows()
cap.release()

max_stars_repo_name: akilawickey/irobot

max_issues_repo_name: akilawickey/irobot

max_forks_repo_name: akilawickey/irobot


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

# See: http://pandas.pydata.org/pandas-docs/stable/10min.html#min

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

# Create a "Series"
s = pd.Series([1, 3, 5, np.nan, 6, 8])
print(s)
print()

# Create a "DataFrame"
dates = pd.date_range('20130101', periods=6)
print(dates)
print()

df = pd.DataFrame(np.random.randn(6,4), index=dates, columns=list('ABCD'))
print(df)
print()

print(df.dtypes)
print()

max_stars_repo_name: jeremiedecock/snippets

max_issues_repo_name: jeremiedecock/snippets

max_forks_repo_name: jeremiedecock/snippets