text:
from shapenet.layer import HomogeneousShapeLayer
import torch
import numpy as np
import pytest
import warnings

@pytest.mark.parametrize("shapes,n_dims,use_cpp,params,target",
                         [
                             (np.zeros((45, 128, 2)), 2, False,
                              torch.ones(10, 44 + 1 + 1 + 2, 1, 1),
                              torch.ones(10, 128, 2)),
                             (np.zeros((45, 128, 2)), 2, True,
                              torch.ones(10, 44 + 1 + 1 + 2, 1, 1),
                              torch.ones(10, 128, 2)),
                             (np.zeros((45, 128, 3)), 3, False,
                              torch.ones(10, 44 + 3 + 3 + 3, 1, 1),
                              torch.ones(10, 128, 3)),
                             (np.zeros((45, 128, 3)), 3, True,
                              torch.ones(10, 44 + 3 + 3 + 3, 1, 1),
                              torch.ones(10, 128, 3))
                         ]
                         )
def test_homogeneous_shape_layer(shapes, n_dims, use_cpp, params, target):
    layer = HomogeneousShapeLayer(shapes, n_dims, use_cpp)
    params.requires_grad_(True)

    assert (layer(params.float()) == target.float()).all()
    try:
        result = layer(params.float())
        result.backward(torch.ones_like(result))
    except:
        assert False, "Backward not successful"


text:
# Load libraries
import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets import load_digits
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import validation_curve

# Load data
digits = load_digits()

# Create feature matrix and target vector
X, y = digits.data, digits.target

# Create range of values for parameter
param_range = np.arange(1, 250, 2)

# Calculate accuracy on training and test set using range of parameter values
train_scores, test_scores = validation_curve(RandomForestClassifier(),
                                             X,
                                             y,
                                             param_name="n_estimators",
                                             param_range=param_range,
                                             cv=3,
                                             scoring="accuracy",
                                             n_jobs=-1)

# Calculate mean and standard deviation for training set scores
train_mean = np.mean(train_scores, axis=1)
train_std = np.std(train_scores, axis=1)

# Calculate mean and standard deviation for test set scores
test_mean = np.mean(test_scores, axis=1)
test_std = np.std(test_scores, axis=1)

# Plot mean accuracy scores for training and test sets
plt.plot(param_range, train_mean, label="Training score", color="black")
plt.plot(param_range, test_mean, label="Cross-validation score", color="dimgrey")

# Plot accurancy bands for training and test sets
plt.fill_between(param_range, train_mean - train_std, train_mean + train_std, color="gray")
plt.fill_between(param_range, test_mean - test_std, test_mean + test_std, color="gainsboro")

# Create plot
plt.title("Validation Curve With Random Forest")
plt.xlabel("Number Of Trees")
plt.ylabel("Accuracy Score")
plt.tight_layout()
plt.legend(loc="best")
plt.show()

ext: py

max_stars_repo_path: Model Evaluation/validation-curve.py

max_stars_repo_name: WyckliffeAluga/data-chronicles

max_issues_repo_path: Model Evaluation/validation-curve.py

max_issues_repo_name: WyckliffeAluga/data-chronicles

max_forks_repo_path: Model Evaluation/validation-curve.py

max_forks_repo_name: WyckliffeAluga/data-chronicles

reason: import numpy


text:
# -*- coding: utf-8 -*-
"""
Created on Sun Sep 23 17:18:21 2018

@author: Florian
"""

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import glob
import sys
import traceback
from os import listdir

import fabio
import h5py
import matplotlib.pyplot as plt
import numpy
from scipy import optimize

from .fio_reader import read

class Scan(object):

    def __init__(self):
        self.filename = None

        self.scan_motor_names = []
        self.scan_motors = []

        self.scan_cmd = None
        self.header_info = None

        self.is_lambda = False
        self.is_p100k = False
        self.is_p300k = False
        self.is_eiger = False
        self.is_mythen = False
        self.is_pe = False

        self.detectors = []

        self.image_data = {}

        self.remove_idxes = None

    def load_scan(self, filename=None, auto_load_images=True, auto_remove_double=False):

        self.scan_motor_names = []
        self.scan_motors = []

        self.is_lambda = False
        self.is_p100k = False
        self.is_p300k = False
        self.is_eiger = False
        self.is_mythen = False
        self.is_pe = False

        self.doubles_removed = False

        self.detectors = []
        self.image_data = {}

        if not filename is None:
            self.filename = filename

        if self.filename is None:
            raise Exception("filename is not defined")

        self.remove_idxes = None

        motor_positions, column_names, data, header_info = read(self.filename)

        self.scan_name = self.filename.split('/')[-1]
        self.scan_name = self.scan_name.rsplit('.', 1)[0]
        self.scan_dir = self.filename.rsplit('/', 1)[0]

        self.scan_cmd = header_info["scan_cmd"]
        self.motor_positions = motor_positions
        self.header_info = header_info
        self.data = data

        # extract rois from header_info
        if "rois" in self.header_info:
            self.rois = self.header_info["rois"]
        else:
            self.rois = None
        # extract UB matrix from header_info
        if "ubmatrix" in self.header_info:
            self.ubmatrix = self.header_info["ubmatrix"]
        else:
            self.ubmatrix = None

        # extract motor names:
        if not self.scan_cmd.find('lscan') == -1 and self.scan_cmd.find('hklscan') == -1:
            self.scan_motor_names = ['l']
            self.scan_motors = [numpy.array(data["kozhue6cctrl_l"])]
        elif not self.scan_cmd.find('kscan') == -1:
            self.scan_motor_names = ['k']
            self.scan_motors = [numpy.array(data["kozhue6cctrl_k"])]
        elif not self.scan_cmd.find('hscan') == -1:
            self.scan_motor_names = ['h']
            self.scan_motors = [numpy.array(data["kozhue6cctrl_h"])]
        elif not self.scan_cmd.find('hklscan') == -1:
            self.scan_motor_names = ['h', 'k', 'l']
            self.scan_motors = [numpy.array(data["kozhue6cctrl_h"]), numpy.array(data["kozhue6cctrl_k"]),
                                numpy.array(data["kozhue6cctrl_l"])]
        elif not self.scan_cmd.find('ascan') == -1 or not self.scan_cmd.find('dscan') == -1:
            mot1_name = self.scan_cmd.split()[1]
            self.scan_motor_names = [mot1_name]
            self.scan_motors = [numpy.array(data[mot1_name])]
        elif not self.scan_cmd.find('a2scan') == -1 or not self.scan_cmd.find('d2scan') == -1:
            mot1_name = self.scan_cmd.split()[1]
            mot2_name = self.scan_cmd.split()[4]
            self.scan_motor_names = [mot1_name, mot2_name]
            self.scan_motors = [numpy.array(data[mot1_name]), numpy.array(data[mot2_name])]
        elif not self.scan_cmd.find('a3scan') == -1 or not self.scan_cmd.find('d3scan') == -1:
            mot1_name = self.scan_cmd.split()[1]
            mot2_name = self.scan_cmd.split()[4]
            mot3_name = self.scan_cmd.split()[7]
            self.scan_motor_names = [mot1_name, mot2_name, mot3_name]
            self.scan_motors = [numpy.array(data[mot1_name]), numpy.array(data[mot2_name]),
                                numpy.array(data[mot3_name])]

        # check which 1D/2D detectors were used during the scan
        self.detectors = []

        for col_name in column_names:
            if not col_name.find('lambda') == -1:
                self.is_lambda = True
            elif not col_name.find('p100k') == -1:
                self.is_p100k = True
            elif not col_name.find('p300k') == -1:
                self.is_p300k = True
            elif not col_name.find('eiger') == -1:
                self.is_eiger = True
            elif not (col_name.find('mythint') == -1 and col_name.find('mythmax') == -1 and col_name.find(
                    'mythroi1') == -1 and col_name.find('mythroi2') == -1):
                self.is_mythen = True
            elif not col_name.find('pe_roi') == -1 or not col_name.find('pe_trigger') == -1:
                self.is_pe = True

        if self.is_lambda:
            self.detectors.append('lambda')
        if self.is_p100k:
            self.detectors.append('p100k')
        if self.is_p300k:
            self.detectors.append('p300k')
        if self.is_eiger:
            self.detectors.append('eiger')
        if self.is_mythen:
            self.detectors.append('mythen')
        if self.is_pe:
            self.detectors.append('pe')

        # load images if autoload is active
        if auto_load_images:
            self.load_image_stack()

        # remove doubled data points if autoremove is active
        if auto_remove_double:
            self.remove_double()

    def load_image_stack(self):
        '''
        loads the detector files for the detectors used in the scan.
        data is stored in the image_data dictionary
        '''

        if self.is_lambda:
            lambda_data = self._load_lambda()
            self.image_data["lambda"] = self._remove_double_fromstack(lambda_data)
        if self.is_p100k:
            p100k_data = self._load_pilatus("p100k")
            self.image_data["p100k"] = self._remove_double_fromstack(p100k_data)
        if self.is_eiger:
            eiger_data = self._load_eiger()
            self.image_data["eiger"] = self._remove_double_fromstack(eiger_data)
        if self.is_mythen:
            mythen_data = self._load_mythen("mythen")
            self.image_data["mythen"] = mythen_data

    def _remove_double_fromstack(self, stack):
        '''
        remove doubled data points from an image stack if remove_idxes is filled.
        i.e. if remove_double has already been used.
        This should only be used internaly if image files are loaded AFTER remove_double
        has been executed
        '''
        if not self.remove_idxes is None:
            return numpy.delete(stack, self.remove_idxes, axis=0)
        else:
            return stack

    def _load_eiger(self):
        '''
        load h5 file from a eiger1m detector und return a numpy 3D array with the data.

        If an index column for the detector is used, only files with corresponding index are load.
        i.e. data index and image index are matching

        Currently, only one h5 file per scan is supported!
        '''

        eiger_folder = "eiger1m"

        print("ScanName: %s" % self.scan_name)

        if len(self.scan_name.split('.')) == 1:
            folder_name = self.scan_name
        else:
            folder_part1 = self.scan_name.split('.')[0]
            folder_part2 = self.scan_name.split('.')[1].rsplit('_', 1)[1]
            folder_name = "%s_%s" % (folder_part1, folder_part2)

        eiger_file_path = "%s/%s/%s" % (self.scan_dir, folder_name, eiger_folder)

        self.eiger_file = "%s/%s_master.h5" % (eiger_file_path, folder_name)

        print("eiger_file: %s" % self.eiger_file)

        try:
            nxsfile = h5py.File(self.eiger_file, 'r')
        except:
            print("could not load hdf5 file. search for master file")

            files = listdir(eiger_file_path)

            for filename in files:
                if "_master.h5" in filename:
                    self.eiger_file = "%s/%s" % (eiger_file_path, filename)
                    nxsfile = h5py.File(self.eiger_file, 'r')
                    break

        nxsdata = nxsfile.get('/entry/data/data_000001')

        # print( nxsdata.shape )

        eiger_idx = None
        try:
            eiger_idx = self.data["eiger_index"]
        except:
            print("no index column for eiger available")

        if not eiger_idx is None:
            nr_images_in_stack = nxsdata.shape[0]
            remove_idxes = []
            for idx in range(nr_images_in_stack):
                if not idx in eiger_idx:
                    remove_idxes.append(idx)

            nxsdata = numpy.delete(nxsdata, remove_idxes, axis=0)
            # print( nxsdata.shape )

        return nxsdata

    def _load_lambda(self):
        '''
        load nxs file from a lambda detector und return a numpy 3D array with the data.

        If an index column for the detector is used, only files with corresponding index are load.
        i.e. data index and image index are matching

        Currently, only one nxs file per scan is supported!
        '''

        self.lambda_file = "%s/%s/lambda/%s_00000.nxs" % (self.scan_dir, self.scan_name, self.scan_name)
        # print ( "lambda_file: %s" % self.lambda_file)

        nxsfile = h5py.File(self.lambda_file)
        nxsdata = nxsfile.get('/entry/instrument/detector/data')

        # print( nxsdata.shape )

        lambda_idx = None
        try:
            lambda_idx = self.data["lambda_index"]
        except:
            print("no index column for lambda available")

        if not lambda_idx is None:
            nr_images_in_stack = nxsdata.shape[0]
            remove_idxes = []
            for idx in range(nr_images_in_stack):
                if not idx in lambda_idx:
                    remove_idxes.append(idx)

            nxsdata = numpy.delete(nxsdata, remove_idxes, axis=0)
            # print( nxsdata.shape )

        return nxsdata

    def _load_pilatus(self, detector):
        '''
        load tif or cbf files from a pilatus detector und return a numpy 3D array with the data.

        If an index column for the detector is used, only files with corresponding index are load.
        i.e. data index and image index are matching
        '''

        self.pilatus_dir = "%s/%s/%s/" % (self.scan_dir, self.scan_name, detector)

        file_list = glob.glob("%s%s*.tif" % (self.pilatus_dir, self.scan_name))
        file_list.extend(glob.glob("%s%s*.cbf" % (self.pilatus_dir, self.scan_name)))
        file_list = sorted(file_list)

        # print (file_list)
        nr_files = len(file_list)

        pilatus_idx = None
        try:
            pilatus_idx = self.data["%s_index" % detector]
            nr_files = len(pilatus_idx)
        except:
            print("no index column for %s available" % detector)

        firstdata = fabio.open(file_list[0]).data

        pilatus_data = numpy.zeros([nr_files, firstdata.shape[0], firstdata.shape[1]], dtype=firstdata.dtype)

        if not pilatus_idx is None:
            cnt = 0
            for idx in pilatus_idx:
                filename = "%s%s_%05d" % (self.pilatus_dir, self.scan_name, idx)
                if (filename + ".tif") in file_list:
                    pilatus_data[cnt, :, :] = fabio.open(filename + ".tif").data
                    cnt += 1
                elif (filename + ".cbf") in file_list:
                    pilatus_data[cnt, :, :] = fabio.open(filename + ".tif").data
                    cnt += 1

        # print( pilatus_data.shape )
        # print( len(pilatus_idx) )

        return pilatus_data

    def _load_mythen(self, detector):
        '''
        load .raw files from mythen detector
        '''

        self.mythen_dir = "%s/%s/%s/" % (self.scan_dir, self.scan_name, detector)

        file_list = glob.glob("%s%s*.raw" % (self.mythen_dir, self.scan_name))
        file_list = sorted(file_list)

        first_file = self._load_1D(file_list[0])

        # print (file_list)
        nr_files = len(file_list)

        image = numpy.zeros([nr_files, len(first_file)])

        cnt = 0
        for filename in sorted(file_list):
            image[cnt, :] = self._load_1D(filename)
            cnt += 1

        return image

    def _load_1D(self, filename, columns=2):

        data = numpy.loadtxt(filename)

        if columns == 2:
            return data[:, 1]

    def remove_double(self):

        if self.doubles_removed:
            return

        pos = self.scan_motors[0]

        remove_idxes = []

        last_pos = numpy.nan

        for idx in range(len(pos)):

            if idx > 0:
                if numpy.abs(pos[idx] - last_pos) == 0:
                    remove_idxes.append(idx - 1)

            last_pos = pos[idx]

        # print( "remove_idxes: %s" % repr(remove_idxes) )

        for key in self.data:
            self.data[key] = numpy.delete(self.data[key], remove_idxes)
            # print("removing from: %s (now: %d)" % (key, len(self.data[key]) ))

        for idx in range(len(self.scan_motors)):
            self.scan_motors[idx] = numpy.delete(self.scan_motors[idx], remove_idxes)

        self.remove_idxes = remove_idxes

        for key in self.image_data:
            self.image_data[key] = numpy.delete(self.image_data[key], remove_idxes, axis=0)
            print("stack size: %s" % repr(self.image_data[key].shape))

        self.doubles_removed = True

    def calc_q(self):

        energy = self.motor_positions["energyfmb"]

        wavelength = 12380 / energy

        k = 2 * numpy.pi / wavelength

        if 'tt' in self.scan_motor_names and 'om' in self.scan_motor_names:
            print('This is a theta2theta scan')

            for idx in range(len(self.scan_motor_names)):
                if self.scan_motor_names[idx] == 'tt':
                    break

            pos = self.scan_motors[idx]

            q = 2 * k * numpy.sin(pos / 2 * numpy.pi / 180)

            print("k = %f" % k)

            self.scan_motor_names.append('q')
            self.scan_motors.append(q)

    def __str__(self):

        s = "Scan Object" + "\n"
        s += "filename:    %s " % self.filename + "\n"
        s += "scan cmd:    %s " % self.scan_cmd + "\n"
        s += "scan motors: %s " % self.scan_motor_names + "\n"
        # s += "motor_positions: %s " % self.motor_positions + "\n"
        s += "header_info: %s " % self.header_info + "\n"
        s += "detectors: %s " % self.detectors + "\n"

        return s

class ScanAnalyzer(object):

    def __init__(self):
        pass

    @staticmethod
    def export2ascii(scan, detector, roi, filename, xcolName=None, bck=None):

        xcol_idx = 0

        if not xcolName is None:
            for idx in range(len(scan.scan_motor_names)):
                if scan.scan_motor_names[idx] == xcolName:
                    xcol_idx = idx
                    break

        result = ScanAnalyzer.extract_rois(scan, {"fit_roi": roi})
        roi_intensity = numpy.array(result[detector]["fit_roi"])

        if not bck is None:
            result = ScanAnalyzer.extract_rois(scan, {"bck_int": bck})
            bck_int = numpy.array(result[detector]["bck_int"])

            roi_intensity = roi_intensity - bck_int

        xcol = scan.scan_motors[xcol_idx]
        xcol_name = scan.scan_motor_names[xcol_idx]

        data = open(filename, 'w')
        for pos, inten in zip(xcol, roi_intensity):
            # print("%6.4f, %6.4f" % (pos, inten-bck_inten))
            data.write("%6.4f %f\n" % (pos, inten))
        data.close()

    @staticmethod
    def extract_rois(scan, rois):

        image_data = scan.image_data

        if "atten_position" in scan.data:
            atten = scan.data["atten_position"]
        else:
            atten = 1

        roi_intensities = {}

        for detector in image_data:

            img_stack = image_data[detector]

            if len(img_stack.shape) == 3:
                print("2D detector image stack")

                img_nr = img_stack.shape[0]

                print("imgs in stack: %d" % img_nr)

                roi_intensity = {}
                roi_intensity["all"] = []

                for key in rois:
                    roi_intensity[key] = []

                for idx in range(img_nr):

                    img = img_stack[idx]

                    key = "all"
                    roi_intensity[key].append(numpy.sum(img, axis=None))

                    for key in rois:
                        roi_img = img[rois[key][1]:rois[key][3], rois[key][0]:rois[key][2]]
                        roi_sum = numpy.sum(roi_img, axis=None)
                        roi_intensity[key].append(roi_sum)

                for key in roi_intensity:
                    roi_intensity[key] = numpy.array(roi_intensity[key]) * numpy.array(atten)
                roi_intensities[detector] = roi_intensity

            elif len(img_stack.shape) == 2:

                print("line detector image")

                img_nr = img_stack.shape[0]

                print("imgs in stack: %d" % img_nr)

                roi_intensity = {}
                roi_intensity["all"] = []

                for key in rois:
                    roi_intensity[key] = []

                for idx in range(img_nr):

                    img = img_stack[idx]

                    key = "all"
                    roi_intensity[key].append(numpy.sum(img, axis=None))

                    for key in rois:
                        roi_img = img[rois[key][0]:rois[key][1]]
                        roi_sum = numpy.sum(roi_img, axis=None)
                        roi_intensity[key].append(roi_sum)

                for key in roi_intensity:
                    roi_intensity[key] = numpy.array(roi_intensity[key]) * numpy.array(atten)
                roi_intensities[detector] = roi_intensity

        return roi_intensities

    @staticmethod
    def combine_scan_image_regions(scan_list, comb_rois):
        for scan in scan_list:
            print(scan)
            scan.remove_double()
        print(comb_rois)

        if len(scan_list[0].detectors) == 0:
            raise Exception("no 2D/1D detectors used in scan")

        detector = scan_list[0].detectors[0]

        detector_images = scan_list[0].image_data[detector]

        master_img = numpy.zeros(detector_images.shape)
        for idx in range(len(scan_list)):

            if "atten_position" in scan_list[idx].data:
                atten = numpy.array(scan_list[idx].data["atten_position"])
            else:
                atten = 1

            detector_images = scan_list[idx].image_data[detector]

            for idx2 in range(detector_images.shape[0]):
                detector_images[idx2] = atten[idx2] * detector_images[idx2]

                # atten_image = numpy.tensordot(atten,detector_images[:,comb_rois[idx][1]:comb_rois[idx][3], comb_rois[idx][0]:comb_rois[idx][2]], axes=([0],[0]))

            atten_image = detector_images[:, comb_rois[idx][1]:comb_rois[idx][3], comb_rois[idx][0]:comb_rois[idx][2]]

            master_img[:, comb_rois[idx][1]:comb_rois[idx][3], comb_rois[idx][0]:comb_rois[idx][2]] = atten_image

        print(master_img[0])

        combined_scan = Scan()

        combined_scan.scan_motors = scan_list[0].scan_motors
        combined_scan.scan_motor_names = scan_list[0].scan_motor_names
        combined_scan.detectors = [detector]
        combined_scan.image_data = {detector: master_img}
        combined_scan.doubles_removed = True

        combined_scan.data = {}

        return combined_scan

    @staticmethod
    def writeHDF(filename, scan):

        image_stack = scan.image_data[scan.detectors[0]]

        with h5py.File(filename, 'w') as f:
            dset = f.create_dataset("data", (image_stack.shape), dtype='float')

            for idx in range(image_stack.shape[0]):
                # print ("min: %f   max: %f" % (numpy.min(image_stack[idx,:,:]), numpy.max(image_stack[idx,:,:]) ) )
                dset[idx, :, :] = image_stack[idx, :, :]

        print("hdf writing complete")

class ScanFitter(object):

    def __init__(self):
        pass

    @staticmethod
    def simpleFit(scan, detector, roi, peaks):

        result = ScanAnalyzer.extract_rois(scan, {"fit_roi": roi})

        roi_intensity = numpy.array(result[detector]["fit_roi"])

        xcol = scan.scan_motors[0]
        xcol_name = scan.scan_motor_names[0]

        params = []

        peak_list = []

        for key in peaks:

            if peaks[key]["type"] == "gauss":
                params += peaks[key]["params"][0:3]

                peak_list.append(key)

        print(params)
        params += [peaks[key]["params"][3]]

        print(params)

        # fitfunc = lambda p, x: ScanFitter._gauss(x, p) # Target function
        # errfunc = lambda p, x, y: fitfunc(p, x) - y # Distance to the target function

        fitfunc = lambda p, x: ScanFitter._multi_gauss(x, p)  # Target function
        errfunc = lambda p, x, y: numpy.abs(
            numpy.log10(fitfunc(p, x)) - numpy.log10(y))  # Distance to the target function

        p1, success = optimize.leastsq(errfunc, params[:], args=(xcol, roi_intensity))

        print("p1 : " + repr(p1))

        fitfunc(params, xcol)

        # ydata = {'data' : roi_intensity, 'before_fit' : ScanFitter._multi_gauss(xcol, params), 'after_fit' : ScanFitter._multi_gauss(xcol, p1)  }
        # ydata = {'data' : roi_intensity, 'before_fit' : ScanFitter._multi_gauss(xcol, params) }

        ydata = {'data': roi_intensity, 'fit': ScanFitter._multi_gauss(xcol, p1)}

        for idx in range(len(peak_list)):
            name = peak_list[idx]

            print("% d: %s" % (idx, name))

            gauss_params = p1[idx * 3:idx * 3 + 3]

            gauss_params = numpy.append(gauss_params, p1[-1])

            print("gauss_params: " + repr(gauss_params))

            result = ScanFitter._gauss(xcol, gauss_params)

            ydata[name] = result

        ScanPlot.plotCurves(xcol, ydata)

    @staticmethod
    def gaussLorentzianFit(scan, detector, roi, peaks, wichtung_params=None):

        result = ScanAnalyzer.extract_rois(scan, {"fit_roi": roi})

        roi_intensity = numpy.array(result[detector]["fit_roi"])

        xcol = scan.scan_motors[0]
        xcol_name = scan.scan_motor_names[0]

        params = [0, 0, 0, 0, 0, 0, 0]

        peak_list = []

        if not len(peaks) == 2:
            raise Exception("online one gauss and one lorentian are alowed for gausLorentzianFit")

        gauss_loaded = False
        lorentzian_loaded = False

        for key in peaks:

            if peaks[key]["type"] == "gauss":

                if gauss_loaded:
                    raise Exception("online one gauss and one lorentian are alowed for gausLorentzianFit")

                params[3:6] = peaks[key]["params"][0:3]

                peak_list.append(key)

                gauss_loaded = True

            if peaks[key]["type"] == "lorentzian":

                if lorentzian_loaded:
                    raise Exception("online one gauss and one lorentian are alowed for gausLorentzianFit")

                params[0:3] = peaks[key]["params"][0:3]

                params[6] = peaks[key]["params"][3]

                peak_list.append(key)

                lorentzian_loaded = False

        print(params)

        if wichtung_params is None:
            wichtung_params = [0, 1, 0]

        wichtung = lambda p, x: p[2] * numpy.heaviside(x - p[0], 1) * numpy.heaviside(-x + p[1], 1) + 1

        fitfunc = lambda p, x: ScanFitter._gauss_lorentzian(x, p)  # Target function
        # errfunc = lambda p, x, y: numpy.abs(numpy.log10(fitfunc(p, x)) - numpy.log10(y)) # Distance to the target function

        errfunc = lambda p, x, y: wichtung(wichtung_params, x) * numpy.abs(
            numpy.log10(fitfunc(p, x) + 1) - numpy.log10(y + 1))  # Distance to the target function
        # errfunc = lambda p, x, y: numpy.abs( (fitfunc(p, x)+1) - (y+1) ) # Distance to the target function

        p1, success = optimize.leastsq(errfunc, params[:], args=(xcol, roi_intensity))

        # p1, success = optimize.curve_fit(fitfunc, xcol, roi_intensity)

        print("p1 : " + repr(p1))

        fitfunc(params, xcol)

        ydata = {'data': roi_intensity, 'fit': ScanFitter._gauss_lorentzian(xcol, p1),
                 'start': ScanFitter._gauss_lorentzian(xcol, params)}
        # ydata = {'data' : roi_intensity, 'start' : ScanFitter._gauss_lorentzian(xcol, params)  }
        # ydata = {'data' : roi_intensity, 'fit' : ScanFitter._gauss_lorentzian(xcol, p1)  }

        for idx in range(len(peak_list)):

            name = peak_list[idx]

            print("% d: %s" % (idx, name))

            if peaks[name]["type"] == "gauss":
                params = p1[3:6]
                params = numpy.append(params, p1[-1])
                print("params: " + repr(params))
                result = ScanFitter._gauss(xcol, params)
            elif peaks[name]["type"] == "lorentzian":
                params = p1[0:3]
                params = numpy.append(params, p1[-1])
                print("params: " + repr(params))
                result = ScanFitter._lorentzian(xcol, params)

            ydata[name] = result

        ScanPlot.plotCurves(xcol, ydata)

    @staticmethod
    def multiVoigtFit(scan, detector, roi, peaks):

        result = ScanAnalyzer.extract_rois(scan, {"fit_roi": roi})

        roi_intensity = numpy.array(result[detector]["fit_roi"])

        xcol = scan.scan_motors[0]
        xcol_name = scan.scan_motor_names[0]

        params = []

        bounds_low = []
        bounds_high = []

        '''
            0 : position
            1 : gauss_intensity
            2 : gauss_width
            3 : lorentzian_intensity
            4 : lorentzian_width
        '''

        peak_list = []

        for key in peaks:

            if peaks[key]["type"] == "voigt":
                params += peaks[key]["params"][0:5]

                peak_list.append(key)

                bounds_low += [10, 0, 0, 0, 0]
                bounds_high += [30, 1e10, 20, 1e10, 20]

        print(params)
        params += [peaks[key]["params"][5]]
        bounds_low += [0]
        bounds_high += [1e5]

        print(params)

        # fitfunc = lambda p, x: ScanFitter._gauss(x, p) # Target function
        # errfunc = lambda p, x, y: fitfunc(p, x) - y # Distance to the target function

        fitfunc = lambda p, x: ScanFitter._multi_voigt(x, p)  # Target function
        # errfunc = lambda p, x, y: numpy.abs(numpy.log10(fitfunc(p, x)+.1) - numpy.log10(y+.1)) # Distance to the target function

        # errfunc = lambda p: numpy.sum(numpy.abs(numpy.log10(fitfunc(p, xcol)+.1) - numpy.log10(roi_intensity+.1)))
        # errfunc = lambda p: ScanFitter.diff_func( fitfunc(p, xcol)+.1, roi_intensity+.1)

        # p1, success = optimize.least_squares(errfunc, params[:], args=(xcol, roi_intensity), bounds=(bounds_low,bounds_high))
        # p1, success = optimize.least_squares(errfunc, params[:], args=(xcol, roi_intensity))
        # p1, success = optimize.leastsq(errfunc, params[:], args=(xcol, roi_intensity))

        # print (errfunc(params[:]))

        # print("errorfunc:" +  str(errfunc(params[:])))

        # p1, success = optimize.least_squares(errfunc, params[:], args=(xcol, roi_intensity))
        print(len(params))

        try:
            # p1, success = optimize.least_squares(errfunc, params[:])

            param_tuple = tuple(params)
            print("param_tuple: " + str(param_tuple))

            # p1, success = optimize.curve_fit(ScanFitter._multi_voigt, xcol, roi_intensity, p0 = param_tuple)

            fit_fun2 = lambda xcol, *p: numpy.log10(ScanFitter._multi_voigt(xcol, *p) + 1)

            # print ("fit_fun2 : " + repr(fit_fun2(xcol, *param_tuple) ) )

            # p2, p2cov = optimize.curve_fit(fit_fun2, xcol, numpy.log10(roi_intensity+1), p0 = param_tuple, bounds=(bounds_low,bounds_high), method='dogbox')
            p2, p2cov = optimize.curve_fit(fit_fun2, xcol, numpy.log10(roi_intensity + .1), p0=param_tuple)

            print("p2cov: " + str(p2cov))

        except:
            etype, value, etraceback = sys.exc_info()
            print(etype)
            print(value)
            print("traceback: " + str(etraceback))

            traceback.print_exc()
            return

        # print ("p1 : " + repr(p1) )

        # ydata = {'data' : roi_intensity, 'before_fit' : ScanFitter._multi_voigt(xcol, params), 'after_fit' : ScanFitter._multi_voigt(xcol, p1)  }
        # ydata = {'data' : roi_intensity, 'before_fit' : ScanFitter._multi_voigt(xcol, tuple(params) ), 'after_fit' : ScanFitter._multi_voigt(xcol, tuple(p1))  }
        # ydata = {'data' : roi_intensity, 'before_fit' : ScanFitter._multi_voigt(xcol, *param_tuple ), 'curve_fit' : ScanFitter._multi_voigt(xcol, *p1), 'log_fit' : ScanFitter._multi_voigt(xcol, *p2)  }
        ydata = {'data': roi_intensity, 'before_fit': ScanFitter._multi_voigt(xcol, *param_tuple),
                 'log_fit': ScanFitter._multi_voigt(xcol, *p2)}
        # ydata = {'data' : roi_intensity, 'before_fit' : ScanFitter._multi_gauss(xcol, params) }

        # ydata = {'data' : roi_intensity, 'fit' : ScanFitter._multi_voigt(xcol, p1)  }

        for idx in range(len(peak_list)):
            name = peak_list[idx]

            print("% d: %s" % (idx, name))

            voigt_params = p2[idx * 5:idx * 5 + 5]

            voigt_params = numpy.append(voigt_params, p2[-1])

            print("voigt_params: " + repr(voigt_params) + "type: " + str(type(voigt_params)))

            result = ScanFitter._pseudo_voigt(xcol, voigt_params)

            ydata[name] = result

        ScanPlot.plotCurves(xcol, ydata)

    @staticmethod
    def diff_func(array1, array2):

        # print(" %6.4f %6.4f " % (sum(numpy.log10(array1)),sum(numpy.log10(array2))) )

        diff = numpy.sum(numpy.abs(numpy.log10(array1) - numpy.log10(array2)))

        # print (diff)

        return diff

    @staticmethod
    def multiFit(scan, detector, roi, peaks):

        result = ScanAnalyzer.extract_rois(scan, {"fit_roi": roi})

        roi_intensity = numpy.array(result[detector]["fit_roi"])

        xcol = scan.scan_motors[0]
        xcol_name = scan.scan_motor_names[0]

        params = []

        peak_list = []

        for key in peaks:

            if peaks[key]["type"] == "gauss":
                params += peaks[key]["params"][0:3]

                peak_list.append(key)

        print(params)
        params += [peaks[key]["params"][3]]

        print(params)

        fitfunc = lambda p, x: ScanFitter._multi_gauss(x, p)  # Target function
        errfunc = lambda p, x, y: numpy.sum(
            numpy.abs(numpy.log10(fitfunc(p, x)) - numpy.log10(y)))  # Distance to the target function

        # p1, success = optimize.leastsq(errfunc, params[:], args=(xcol, roi_intensity))

        p1, success = optimize.curve_fit(fitfunc, xcol, roi_intensity)

        print("p1 : " + repr(p1))

        print(fitfunc(params, xcol))

        ydata = {'data': roi_intensity, 'fit': ScanFitter._multi_gauss(xcol, p1)}

        for idx in range(len(peak_list)):
            name = peak_list[idx]

            print("% d: %s" % (idx, name))

            gauss_params = p1[idx * 3:idx * 3 + 3]

            gauss_params = numpy.append(gauss_params, p1[-1])

            print("gauss_params: " + repr(gauss_params))

            result = ScanFitter._gauss(xcol, gauss_params)

            ydata[name] = result

        ScanPlot.plotCurves(xcol, ydata)

    @staticmethod
    def _multi_gauss(x, params, nr_peaks=None):

        nr_peaks = int(len(params) / 3)

        result = numpy.zeros([len(x), ])

        for idx in range(nr_peaks):
            gauss_params = params[idx * 3:idx * 3 + 3]

            # print("gauss_params: " + repr(gauss_params))

            # gauss_params.append(0.)
            gauss_params = numpy.append(gauss_params, 0.0)
            result += ScanFitter._gauss(x, gauss_params)

        result += params[-1]

        return result

    @staticmethod
    def _multi_voigt(x, *params):

        nr_peaks = int(len(params) / 5)

        # print ("peaks: %d" % nr_peaks)
        # print ("params: " + str(params))

        result = numpy.zeros([len(x), ])

        for idx in range(nr_peaks):
            voigt_params = params[idx * 5:idx * 5 + 5]

            # print("gauss_params: " + repr(gauss_params))

            # gauss_params.append(0.)
            voigt_params = numpy.append(voigt_params, 0.0)
            result += ScanFitter._pseudo_voigt(x, voigt_params)

        result += params[-1]

        # print("multi voigt done")

        return result

    @staticmethod
    def _gauss_lorentzian(x, params):

        if not len(params) == 7:
            raise Exception("ERROR on _gauss_lorentzian: len(params) == 7 required (%d given) " % len(params))

        result = numpy.zeros([len(x), ])

        lorentzian_params = params[0:3]
        lorentzian_params = numpy.append(lorentzian_params, 0.0)

        gauss_params = params[3:6]
        gauss_params = numpy.append(gauss_params, 0.0)

        result += ScanFitter._lorentzian(x, lorentzian_params)
        result += ScanFitter._gauss(x, gauss_params)

        result += params[-1]

        return result

    @staticmethod
    def _gauss(x, params):
        '''
        copy of the MATLAB gauss function:

        function [ int ] = gauss( x, params )
        %GAUSS function
        %   x:      x-values
        %   params: 1 - sigma
        %           2 - shift
        %           3 - intensity
        %           4 - background

        gauss = @(x,params) params(3).*1./sqrt(2*pi*params(1).^2) * exp( -0.5.*((x-params(2))/params(1)).^2 ) + params(4);
        int = gauss(x,params);
        '''

        intensity = params[2] / (numpy.sqrt(2 * numpy.pi * params[0] ** 2)) * numpy.exp(
            -0.5 * ((numpy.array(x) - params[1]) / params[0]) ** 2) + params[3]

        return intensity

    @staticmethod
    def _lorentzian(x, params):
        '''

        lorentzian:

        x:  x-values
        params 0: omega_0   (maximum position)
        params 1: gamma     (width)
        params 2: intensity
        params 4: background

        f(omega) = 1/( (omega**2 - omega_0**2 + gamma**2*omega_0**2) )

        '''

        intensity = params[2] * 1 / ((x ** 2 - params[0] ** 2) ** 2 + params[1] ** 2 * params[0] ** 2) + params[3]

        return intensity

    @staticmethod
    def _pseudo_voigt(x, params):

        '''
        params:
            0 : position
            1 : gauss_intensity
            2 : gauss_width
            3 : lorentzian_intensity
            4 : lorentzian_width
            5 : background
        '''

        gauss_params = [params[2], params[0], params[1], 0]
        lorentzian_params = [params[0], params[4], params[3], 0]

        return ScanFitter._gauss(x, gauss_params) + ScanFitter._lorentzian(x, lorentzian_params) + params[5]

class ScanPlot(object):

    def __init__(self):
        pass

    @staticmethod
    def plotRois(scan, detector, rois, xcolName=None):
        '''
        uses matplotlib to plot roi intensities of scan for given detector

        '''
        xcol_idx = 0

        if not xcolName is None:
            for idx in range(len(scan.scan_motor_names)):
                if scan.scan_motor_names[idx] == xcolName:
                    xcol_idx = idx
                    break

        xcol = scan.scan_motors[xcol_idx]

        xcol_name = scan.scan_motor_names[xcol_idx]

        result = ScanAnalyzer.extract_rois(scan, rois)

        roi_intensity = result[detector]

        nr_plots = len(roi_intensity)

        img_size = scan.image_data[detector].shape[0]

        if "background" in roi_intensity:

            background_area = (rois["background"][2] - rois["background"][0]) * (
                    rois["background"][3] - rois["background"][1])

            print(background_area)

            background = numpy.array(roi_intensity["background"])
            nr_plots -= 1
        else:
            background = 0
            background_area = 1

        sqrt_nr = int(numpy.sqrt(nr_plots) + 0.5)
        cnt = 1
        for key in sorted(roi_intensity):

            if not key == "background":
                print(key)
                plt.subplot(sqrt_nr, sqrt_nr, cnt)

                if not key == "all":
                    roi_area = (rois[key][2] - rois[key][0]) * (rois[key][3] - rois[key][1])
                    print(roi_area)
                else:
                    roi_area = numpy.prod(img_size)

                # plt.semilogy(xcol, roi_intensity[key] )
                plt.semilogy(xcol, numpy.array(roi_intensity[key]) / roi_area - background / background_area)
                plt.ylabel(key)
                plt.xlabel(xcol_name)

                cnt += 1

        plt.show()

    @staticmethod
    def plotCurves(xcol, ydata):

        # print (ydata)

        fig, ax = plt.subplots()

        for key in ydata:
            if key == "substrate":
                ax.semilogy(xcol, numpy.array(ydata[key]), label=key, linestyle=':')
            elif key == "film":
                ax.semilogy(xcol, numpy.array(ydata[key]), label=key, linestyle=':')
            else:
                ax.semilogy(xcol, numpy.array(ydata[key]), label=key, linestyle='-')

        legend = ax.legend(loc='upper left')

        plt.show()

if __name__ == '__main__':

    # filename = '/asap3/petra3/gpfs/p08/2019/data/11006183/raw/cfosto190311_00183.fio'
    # filename = '/asap3/petra3/gpfs/p08/2019/data/11006183/raw/cfosto190311_00187.fio'
    # filename = '/asap3/petra3/gpfs/p08/2019/data/11006183/raw/nfosto190211_00297.fio'
    # rois = {"roi1" : [754, 249, 775, 265], "roi2" : [754, 249, 1500, 265], "roi3" : [800, 249, 1500, 265] }

    filename = '/asap3/petra3/gpfs/p08/2019/commissioning/c20190208_000_start19/raw/ptysz02032019_00423.fio'

    rois = {"roi1": [163, 71, 174, 76], "roi2": [164, 69, 174, 79], "roi3": [0, 0, 487, 194]}
    scan1 = Scan()

    # scan1.load_scan(filename)
    # print (scan1)

    # scan1.remove_double()
    # ScanPlot.plotRois(scan1, "p100k", rois)

    # filenames = ['/asap3/petra3/gpfs/p08/2019/data/11006183/raw/nfomao190322_01039.fio', '/asap3/petra3/gpfs/p08/2019/data/11006183/raw/nfomao190322_01040.fio']
    # rois = {"roi1" : [765, 279, 830, 290], "roi2" : [802, 279, 1500, 290], "roi3" : [820, 279, 1500, 290], "background" : [1000, 279, 1340, 290] }

    filenames = ['/asap3/petra3/gpfs/p08/2019/data/11006183/raw/nfosto190211_00286.fio',
                 '/asap3/petra3/gpfs/p08/2019/data/11006183/raw/nfosto190211_00287.fio']
    rois = {"roi1": [761, 279, 812, 290], "roi2": [761, 279, 1220, 290], "roi3": [940, 279, 1220, 290],
            "background": [1220, 279, 1340, 290]}

    scans = []
    for name in filenames:
        scan = Scan()
        scan.load_scan(name)
        scans.append(scan)
    combined_scan = ScanAnalyzer.combine_scan_image_regions(scans, [[0, 0, 880, 516], [880, 0, 1340, 516]])

    ScanAnalyzer.writeHDF("test.h5", combined_scan)

    print(combined_scan)

    try:
        print("daten: " + combined_scan.data)
    except:
        pass
    try:
        print("filename: " + combined_scan.filename)
    except:
        pass

    # ScanPlot.plotRois(combined_scan, "lambda", rois)

    # ScanFitter.multiFit(combined_scan, "lambda", rois["roi2"], {"film" : {"type" : "gauss", "params" : [0.015,0.95,2e4,4e4]}, "bck" : {"type" : "gauss", "params" : [0.15,0.95,1e3,4e4]}, "substrate" :  {"type" : "gauss", "params" : [0.003,1.01,2e5,1.2e5]} })

    ScanFitter.gaussLorentzianFit(combined_scan, "lambda", rois["roi2"],
                                  {"film": {"type": "gauss", "params": [0.015, 0.95, 2e4, 4e4]},
                                   "substrate": {"type": "lorentzian", "params": [1.01, 0.003, 1e3, 1.2e5]}})

    # filename = '/asap3/petra3/gpfs/p08/2019/data/11005886/raw/M71801_FG966_probeC_00183.fio'
    # rois = {"roi1" : [163, 71, 174, 76], "roi2" : [164, 69, 174, 79], "roi3" : [0, 0, 487, 194] }
    # scan1 = Scan()
    # scan1.load_scan(filename)
    # print(scan1)

max_stars_repo_name: schreiber-lab/mlreflect

max_issues_repo_name: schreiber-lab/mlreflect

max_forks_repo_name: schreiber-lab/mlreflect

reason: import numpy,from scipy


text:
import math
import numpy as np

import torch
from torch import nn
import torch.nn.functional as F

from utils.nar_beam_search import Beamsearch
from utils.functions import get_best

class NARModel(nn.Module):

    def __init__(self,
                 problem,
                 embedding_dim,
                 encoder_class,
                 n_encode_layers,
                 aggregation="sum",
                 normalization="layer",
                 learn_norm=True,
                 track_norm=False,
                 gated=True,
                 n_heads=8,
                 mask_graph=False,
                 *args, **kwargs):
        """
        Models with a GNN/Transformer/MLP encoder and the Non-autoregressive decoder using attention mechanism

        Args:
            problem: TSP/TSPSL, to identify the learning paradigm
            embedding_dim: Hidden dimension for encoder/decoder
            encoder_class: GNN/Transformer/MLP encoder
            n_encode_layers: Number of layers for encoder
            aggregation: Aggregation function for GNN encoder
            normalization: Normalization scheme ('batch'/'layer'/'none')
            learn_norm: Flag for enabling learnt affine transformation during normalization
            track_norm: Flag to enable tracking training dataset stats instead of using batch stats during normalization
            gated: Flag to enbale anisotropic GNN aggregation
            n_heads: Number of attention heads for Transformer encoder/MHA in decoder
            mask_graph: Flag to use graph mask during decoding

        References:
            - A. Nowak, S. Villar, A. S. Bandeira, and J. Bruna. A note on learning algorithms for quadratic assignment with graph neural networks. arXiv preprint arXiv:1706.07450, 2017.
            - C. K. Joshi, T. Laurent, and X. Bresson. An efficient graph convolutional network technique for the travelling salesman problem. arXiv preprint arXiv:1906.01227, 2019.
        """

        assert problem.NAME in ['tsp', 'tspsl'], "NAR Attention Decoder only supports TSP and TSPSL."

        super(NARModel, self).__init__()

        self.problem = problem
        self.embedding_dim = embedding_dim
        self.encoder_class = encoder_class
        self.n_encode_layers = n_encode_layers
        self.aggregation = aggregation
        self.normalization = normalization
        self.learn_norm = learn_norm
        self.track_norm = track_norm
        self.gated = gated
        self.n_heads = n_heads
        self.mask_graph = mask_graph

        # Input embedding layer
        self.init_embed = nn.Linear(2, embedding_dim, bias=True)

        # Encoder model
        self.embedder = self.encoder_class(n_layers=n_encode_layers,
                                           n_heads=n_heads,
                                           hidden_dim=embedding_dim,
                                           aggregation=aggregation,
                                           norm=normalization,
                                           learn_norm=learn_norm,
                                           track_norm=track_norm,
                                           gated=gated)

        # Edge prediction layer
        self.project_node_emb = nn.Linear(embedding_dim, embedding_dim, bias=True)
        self.project_graph_emb = nn.Linear(embedding_dim, embedding_dim, bias=True)
        self.edge_pred = nn.Linear(embedding_dim, 2, bias=True)

    def set_decode_type(self, decode_type, **kwargs):
        self.decode_type = decode_type

    def forward(self, nodes, graph, supervised=False, targets=None, class_weights=None, return_pi=False):
        """
        Args:
            nodes: Input graph nodes (B x V x 2)
            graph: Graph as **NEGATIVE** adjacency matrices (B x V x V)
            supervised: Toggles SL training, teacher forcing and BCE loss computation
            targets: Targets for teacher forcing and BCE loss
            return_pi: Toggles returning the output sequences
                       Not compatible with DataParallel as the results
                       may be of different lengths on different GPUs
        """

        # Supervised learning
        if self.problem.NAME == 'tspsl' and supervised:
            assert targets is not None, "Pass targets during training in supervised mode"
            assert class_weights is not None, "Pass class weights during training in supervised mode for NAR models"

            if return_pi:
                # Perform greedy search during training if we want to log cost, pi during training
                logits, _, pi, cost = self.beam_search(nodes, graph, beam_size=3)
                loss = nn.CrossEntropyLoss(weight=class_weights, reduction='mean')(
                    logits.permute(0, 3, 1, 2), targets)
                return cost, loss, pi

            else:
                # Only perform _inner function for faster training
                logits, _ = self._inner(nodes, graph)
                loss = nn.CrossEntropyLoss(weight=class_weights, reduction='mean')(
                    logits.permute(0, 3, 1, 2), targets)
                return torch.zeros(nodes.shape[0]), loss

        # Reinforcement learning or inference
        else:
            _, log_p, pi, cost = self.greedy_search(nodes, graph)
            ll = self._calc_log_likelihood(log_p[:, :, :, 1], pi)

            if return_pi:
                return cost, ll, pi
            return cost, ll

    def _init_embed(self, nodes):
        return self.init_embed(nodes)

    def _inner(self, nodes, graph):
        """
        Returns:
            logits: Unnormalized logits over graph edges (B x V x V x 2)
            log_p: Log-Softmax over final dimension of `logits` (B x V x V)
        """
        batch_size, num_nodes, _ = nodes.shape

        # Compute node embeddings
        embeddings = self.embedder(self._init_embed(nodes), graph)

        # Compute edge embeddings (B x V x V x H)
        Ux = self.project_node_emb(embeddings)
        Gx = self.project_graph_emb(embeddings.mean(dim=1))
        edge_embeddings = F.relu(Ux[:, :, None, :] + Ux[:, None, :, :] + Gx[:, None, None, :])

        if self.mask_graph:
            edge_embeddings[graph[:, :, :, None].expand_as(edge_embeddings)] = 0

        # Compute logits
        logits = self.edge_pred(edge_embeddings)  # B x V x V x 2
        log_p = F.log_softmax(logits, dim=3)

        return logits, log_p

    def _calc_log_likelihood(self, _log_p, a):
        # Get log_p corresponding to selected actions
        log_p = _log_p.gather(2, a.unsqueeze(-1)).squeeze(-1)
        assert (log_p > -1000).data.all(), "Logprobs should not be -inf, check sampling procedure!"

        # Calculate log_likelihood
        return log_p.sum(1)

    def beam_search(self, nodes, graph, beam_size):
        """Method to perform graph search (beam search or greedy search)

        Args:
            nodes: Input graph nodes (B x V x 2)
            graph: Graph as **NEGATIVE** adjacency matrices (B x V x V)
            beam_size: Beam search width (=1 to enable greedy search)

        Returns:
            logits, log_p: Outputs of inner function
            pi: Tour sequence for shortest out of `beam_size` candidate tours
            cost: Tour length for shortest out of `beam_size` candidate tours
        """
        batch_size, num_nodes, _ = nodes.shape

        # Compute logits
        logits, log_p = self._inner(nodes, graph)

        # Perform beamsearch
        with torch.no_grad():
            _log_p = log_p.clone().detach()[:, :, :, 1]
            _log_p[_log_p == 0] = -1e-10  # Set 0s (i.e. log(1)s) to very small negative number

            beamsearch = Beamsearch(beam_size, batch_size, num_nodes, device=_log_p.device, decode_type=self.decode_type)
            trans_probs = _log_p.gather(1, beamsearch.get_current_state())
            for step in range(num_nodes - 1):
                beamsearch.advance(trans_probs)
                trans_probs = _log_p.gather(1, beamsearch.get_current_state())

            # Find TSP tour with highest probability among beam_size candidates
            ends = torch.zeros(batch_size, 1, device=_log_p.device)
            pi = beamsearch.get_hypothesis(ends)

            if beam_size == 1:
                # Compute tour costs
                cost, _ = self.problem.get_costs(nodes, pi)

            elif beam_size > 1:
                # Beam search
                sequences = []
                costs = []
                ids = []

                # Iterate over all positions in beam
                for pos in range(0, beam_size):
                    ends = pos * torch.ones(batch_size, 1, device=_log_p.device)  # New positions

                    try:
                        pi_temp = beamsearch.get_hypothesis(ends)
                        cost_temp, _ = self.problem.get_costs(nodes, pi_temp)
                        cost_temp, pi_temp = cost_temp.cpu().numpy(), pi_temp.cpu().numpy()

                        sequences.append(pi_temp)
                        costs.append(cost_temp)
                        ids.append(list(range(batch_size)))

                    except AssertionError:
                        # Handles error if the temporary solution is an invalid tour
                        continue

                sequences = np.array(sequences)
                costs = np.array(costs)
                ids = np.array(ids)

                # Reshape/permute sequences, costs, ids
                valid_beam, batch_size, num_nodes = sequences.shape
                sequences = sequences.reshape(valid_beam * batch_size, num_nodes)
                s_idx = []
                for i in range(batch_size):
                    for j in range(valid_beam):
                        s_idx.append(i + j * batch_size)

                # Get sequences and costs of shortest tours
                pi, cost = get_best(sequences[s_idx], costs.T.flatten(), ids.T.flatten())

        return logits, log_p, pi, cost

    def greedy_search(self, nodes, graph):
        return self.beam_search(nodes, graph, beam_size=1)

    def sample_many(self, nodes, graph, batch_rep=1, iter_rep=1):
        assert batch_rep == 1 and iter_rep == 1, "Sampling solutions is not supported. Use beam search instead."
        return self.greedy_search(nodes, graph)

max_stars_repo_name: netspractice/learning-tsp

max_issues_repo_name: netspractice/learning-tsp

max_forks_repo_name: netspractice/learning-tsp


text:
# Copyright 2017 The TensorFlow Authors. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ==============================================================================
"""Tests for window_ops."""

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import functools

import numpy as np

from tensorflow.contrib.signal.python.kernel_tests import test_util
from tensorflow.contrib.signal.python.ops import window_ops
from tensorflow.python.framework import dtypes
from tensorflow.python.framework import ops
from tensorflow.python.platform import test

def _scipy_raised_cosine(length, symmetric=True, a=0.5, b=0.5):
  """A simple implementation of a raised cosine window that matches SciPy.

  https://en.wikipedia.org/wiki/Window_function#Hann_window
  https://github.com/scipy/scipy/blob/v0.14.0/scipy/signal/windows.py#L615

  Args:
    length: The window length.
    symmetric: Whether to create a symmetric window.
    a: The alpha parameter of the raised cosine window.
    b: The beta parameter of the raised cosine window.

  Returns:
    A raised cosine window of length `length`.
  """
  if length == 1:
    return np.ones(1)
  odd = length % 2
  if not symmetric and not odd:
    length += 1
  window = a - b * np.cos(2.0 * np.pi * np.arange(length) / (length - 1))
  if not symmetric and not odd:
    window = window[:-1]
  return window

class WindowOpsTest(test.TestCase):

  def setUp(self):
    self._window_lengths = [1, 2, 3, 4, 5, 31, 64, 128]
    self._dtypes = [(dtypes.float16, 1e-2),
                    (dtypes.float32, 1e-6),
                    (dtypes.float64, 1e-9)]

  def _compare_window_fns(self, np_window_fn, tf_window_fn):
    with self.session(use_gpu=True):
      for window_length in self._window_lengths:
        for periodic in [False, True]:
          for tf_dtype, tol in self._dtypes:
            np_dtype = tf_dtype.as_numpy_dtype
            expected = np_window_fn(window_length,
                                    symmetric=not periodic).astype(np_dtype)
            actual = tf_window_fn(window_length, periodic=periodic,
                                  dtype=tf_dtype).eval()
            self.assertAllClose(expected, actual, tol, tol)

  def test_hann_window(self):
    """Check that hann_window matches scipy.signal.hann behavior."""
    # The Hann window is a raised cosine window with parameters alpha=0.5 and
    # beta=0.5.
    # https://en.wikipedia.org/wiki/Window_function#Hann_window
    self._compare_window_fns(
        functools.partial(_scipy_raised_cosine, a=0.5, b=0.5),
        window_ops.hann_window)

  def test_hamming_window(self):
    """Check that hamming_window matches scipy.signal.hamming's behavior."""
    # The Hamming window is a raised cosine window with parameters alpha=0.54
    # and beta=0.46.
    # https://en.wikipedia.org/wiki/Window_function#Hamming_window
    self._compare_window_fns(
        functools.partial(_scipy_raised_cosine, a=0.54, b=0.46),
        window_ops.hamming_window)

  def test_constant_folding(self):
    """Window functions should be constant foldable for constant inputs."""
    for window_fn in (window_ops.hann_window, window_ops.hamming_window):
      for dtype, _ in self._dtypes:
        for periodic in [False, True]:
          g = ops.Graph()
          with g.as_default():
            window = window_fn(100, periodic=periodic, dtype=dtype)
            rewritten_graph = test_util.grappler_optimize(g, [window])
            self.assertEqual(1, len(rewritten_graph.node))

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

max_stars_repo_name: fraudies/tensorflow

max_issues_repo_name: fraudies/tensorflow

max_forks_repo_name: fraudies/tensorflow


text:
# -*- coding: utf-8 -*-
"""
Created on Wed Jul 12 11:27:14 2017

@author: Alex
"""

import pandas as pd
import numpy as np
from datetime import timedelta
from flarepy import utils
import warnings

def compare_to_primary(data_ref, data, windows=[timedelta(minutes=1), timedelta(minutes=2), timedelta(minutes=3)]):
    """
    Basic function to compare a secondary dataset to a primary reference.

    Parameters
    ----------
    data_ref : ~`pandas.Series` or ~`pandas.Series`
        The primary/reference list of flares indexed by peaktime.
        The actual data valus arn't relivent, it compares based on the index.
        Note: currently expects values in this index to be unique.

    data : ~`pandas.Series`
        The secondary flare list, that will be compared to the primary (reference).

    windows : `list` of `datetime.timedelta`
        The list of windows to search in, these get progressively larger.
        Ideally you want the increase per window to be small enough to only add
        a single match.

    Returns
    -------
    ser_matched : ~`pandas.Series`
        The list of matches between the primary (reference) and the secondary list.
        The index is the index from the secondary, the values are the associated
        index from the primary reference list.

    ser_not_in_primary : ~`pandas.Series`
        The list of flares in the secondary and not in the primary (reference).
        The index is the index from the secondary, the values are all NaN.

    ser_not_in_secondary : ~`pandas.Series`
        The list of flares in the primary (reference) and not in the secondary.
        The index is the index from the secondary, the values are all NaN.
    """
    # Sanatise inputs
    ser_ref = data_ref
    if isinstance(ser_ref, pd.DataFrame):
        ser_ref = pd.Series(data=np.nan, index=ser_ref.index)
    elif isinstance(ser_ref, np.ndarray):
        ser_ref = pd.Series(data=np.nan, index=pd.to_datetime(ser_ref))
    ser_data = data
    if isinstance(ser_data, pd.DataFrame):
        ser_data = pd.Series(data=np.nan, index=ser_data.index)
    elif isinstance(ser_ref, np.ndarray):
        ser_data = pd.Series(data=np.nan, index=pd.to_datetime(ser_data))

    # Check for duplicates, remove if detected
    if not (np.all(np.invert(ser_ref.index.duplicated()))):
        warnings.warn("HEK data has flares with dupicate peak times. First entries have been kept.")
        ser_ref = ser_ref.groupby(ser_ref.index).first()
    if not (np.all(np.invert(ser_data.index.duplicated()))):
        warnings.warn("HEK data has flares with dupicate peak times. First entries have been kept.")
        ser_data = ser_data.groupby(ser_data.index).first()

    # For each of the values in the primary (reference) data
    ser_matched, ser_not_in_primary = compare_ser_results(ser_ref, ser_data, windows=windows)

    # For each of the values in the secondary
    ser_matched_alt, ser_not_in_secondary = compare_ser_results(ser_data, ser_ref, windows=windows)

    # Check if the matched series are equivilent. I'm not sure if this will be true.
    if not ser_matched.equals(ser_matched_alt):
        print('Note: the two matched lists don\'t match.')

    # Return the results
    return ser_matched, ser_not_in_primary, ser_not_in_secondary

def get_col_multiindex(windows, metrics, names=['Window','Metric']):
    """
    A function to create the correct column multiindex for a given set of windows
    and metrics.

    Parameters
    ----------

    data : ~`pandas.Series`
        The secondary flare list, that will be compared to the primary (reference).

    windows : `list` of `datetime.timedelta`
        The list of windows to search in, these get progressively larger.
        Ideally you want the increase per window to be small enough to only add
        a single match.

    names : `list` of `str` (['Window','Metric'])
        The list of names for each level of the multiindex.
        Generally left blank to use the default, otherwise the resulting pandas
        DataFrame may not work with other FlarePy tools, but added to give more
        versatility to the function.

    Returns
    -------
    output : ~`pandas.core.indexes.multi.MultiIndex`
        The multiindex correlating to the given windows and metrics.
    """
    # Add the 0s window if not included.
    if not windows[0].seconds is 0:
        windows.insert(0, timedelta(minutes=0))

    # Create list for each level
    lis_level_0 = []
    lis_level_1 = []
    for window in windows:
        # Window string
        str_win = str(window.seconds) + 's'

        # Add to the two lists
        for i in range(0, len(metrics)):
            lis_level_0.append(str_win)
            lis_level_1.append(metrics[i])

    # return the resulting multiindex
    return pd.MultiIndex.from_tuples(list(zip(*[lis_level_0, lis_level_1])), names=names)

def get_varied_window_stats(ser_data_ref, dic_ser_data, percentage=False, windows=[timedelta(minutes=1), timedelta(minutes=2), timedelta(minutes=3)], keydelimiter=';'):
    """
    Take all the results, perform comparisons and output a multiindex pandas
    DataFrame of them.

    Parameters
    ----------

    ser_data_ref : ~`pandas.Series`
        The primary/reference flare list, that will be used for comparison.

    dic_ser_data : `dict`
        The dictionary of all the results datasets to be compared.
        Note, for good measure it's generally worth having the reference in here.

    percentage : (False)
        ####### Needs to be implmented

    windows : `list` of `datetime.timedelta`
        The list of windows to search in, these get progressively larger.
        Each window will be checked and results added seperately so you can
        compare the results between windows for each method.

    keydelimiter : (';')
        The keys in the dictionary are created as a merger of the pre-processing
        method and method parameters details.
        This is the delimiter used to seperate out these parts.

    Returns
    -------
    ser_matched : ~`pandas.core.indexes.multi.MultiIndex`
        The multiindex correlating to the given windows and metrics.
    """
    #df_cwt_compared_to_hek, ser_in_cwt_not_in_hek, ser_in_hek_not_in_cwt = utils.compare_to_HEK(df_hek, df_peaks_cwt, windows=lis_windows)

    # Create the data and row index array
    lis_data = []
    lis_row_index = [[],[],[]]

    for key, df_detections in dic_ser_data.items():
        # Get the row index parameters and add to the row index array
        lis_index = key.split(keydelimiter)
        lis_row_index[0].append(lis_index[0])
        lis_row_index[1].append(lis_index[1])
        lis_row_index[2].append(lis_index[2])

        # Get the stats data values
        lis_data_line = []
        for i in range(0,len(windows)+1):#td_windows in windows:
            #print('windows[0:i]: ' + str(windows[0:i]))
            ser_matched, ser_not_in_primary, ser_not_in_secondary = compare_to_primary(ser_data_ref, df_detections, windows=windows[0:i])

            # Add the values
            lis_data_line.append(len(ser_not_in_primary))
            lis_data_line.append(len(ser_matched))
            lis_data_line.append(len(ser_not_in_secondary))

        # Add this line to the data
        lis_data.append(lis_data_line)

    # Create the column and row multiindices
    col_index = get_col_multiindex(windows, ['FA', 'M', 'FR'])
    #[np.array(['0s','0s','0s','2s','1s','1s','2s','2s','2s','3s','3s','3s']),np.array(['FA','M','FR','False Acceptance','Matched','False Rejection','False Acceptance','Matched','False Rejection','False Acceptance','Matched','False Rejection'])]
    row_index = pd.MultiIndex.from_tuples(list(zip(*lis_row_index)), names=['Pre-Processing','Method','Parameters'])

    """
    lis_data=[['bin=60s','HEK - Reference','',0,int_hek,0,0,int_hek,0,0,int_hek,0,0,int_hek,0],
          ['bin=60s boxcart=5','4min','N=4',3,12,8,3,12,8,3,12,8,3,12,8],
          ['bin=60s boxcart=5','CWT','W=1…49',4,16,4,4,16,4,4,16,4,4,16,4],['bin=60s boxcart=5','CWT','W=1…100',3,12,8,3,12,8,3,12,8,3,12,8]]
    lis_data=[[2,15,5,3,16,3,4,17,1,5,18,0],[3,12,8,4,13,6,5,14,4,6,15,2],[4,16,4,5,17,2,6,18,0,7,19,0],[3,12,8,4,13,6,5,14,4,6,15,2]]
    col_index = pd.MultiIndex.from_tuples(list(zip(*lis_col_index)), names=['Window','Metric'])
    lis_row_index = [np.array(['bin=12s boxcart=5','bin=60s boxcart=5','bin=60s boxcart=5','bin=60s boxcart=5']),np.array(['4min','4min','CWT','CWT']),np.array(['N=20','N=4','W=1…49','W=1…100'])]
    row_index = pd.MultiIndex.from_tuples(list(zip(*lis_row_index)), names=['Pre-Processing','Method','Parameters'])

    df_temp = pd.DataFrame(data=np.zeros([4,12]), columns=col_index, index=row_index)
    """
    return pd.DataFrame(data=lis_data, columns=col_index, index=row_index)

def compare_ser_results(ser_pri, ser_sec, windows=[timedelta(minutes=1), timedelta(minutes=2), timedelta(minutes=3)]):
    """
    The underlaying method used to compare one pandas.Series to another.
    The algorithm starts by checking for an exact match, if no match is found
    then this checks progressively wider windows (arround the primary/reference)
    as defined by the windows argument.
    Windows should increase in such a way as to reduce the chances of multiple
    matches, if 2 matches are found in the lowest matching window then we return
    the first.
    If we get more then 2 matches we assume it's too big a window and don't
    return matches.

    Parameters
    ----------
    ser_pri : ~`pandas.Series`
        The primary/reference list of flares indexed by peaktime.
        The actual data valus arn't relivent, it compares based on the index.
        Note: currently expects values in this index to be unique.

    ser_sec : ~`pandas.Series`
        The secondary flare list, that will be compared to the primary (reference).

    windows : `list`
        The list of windows to search in, these get progressively larger.
        Ideally you want the increase per window to be small enough to only add
        a single match.

    Returns
    -------
    ser_matched : ~`pandas.Series`
        The list of matches between the primary (reference) and the secondary list.
        The index is the index from the secondary, the values are the associated
        index from the primary (reference) list.

    ser_not_in_primary : ~`pandas.Series`
        The list of flares in the secondary and not in the primary (reference).
        The index is the index from the secondary, the values are all NaN.
    """
    ####print('ser_pri: '+str(len(ser_pri)))
    ####print(ser_pri)
    # Holders for matched/unmatched flares
    ser_matched = pd.Series()
    ser_not_in_primary = pd.Series(data=np.nan, index=ser_sec.index)

    #print('compare_ser_results: windows: '+str(windows))

    # For each of the values in the primary (reference) dataset
    ####print('\n')
    ####print('str(len(ser_sec)): '+str(len(ser_sec)))
    ####print(ser_sec)
    ####print('\n\n')
    for i in range(0, len(ser_sec)):
        dt_det = ser_sec.index[i]
        #print('dt_det: '+str(dt_det)+'\n')

        # Check if the flare detection is the same as a HEK entry
        if dt_det in ser_pri:
            #print('here0')
            ser_matched[dt_det] = dt_det
            ser_not_in_primary = ser_not_in_primary.drop(dt_det)
        else:
            ####print('here1')
            # If it's not exactly the same, check within the given windows
            for td_window in windows:
                ####print('here2')
                dt_window_min = dt_det - td_window
                dt_window_max = dt_det + td_window
                ####print('dt_window_min: '+str(dt_window_min))
                ####print('dt_window_max: '+str(dt_window_max))

                # truncate to this window and see if any HEK results are present
                ser_pri_trun = ser_pri.truncate(dt_window_min, dt_window_max)
                ####print('ser_pri_trun: '+str(ser_pri_trun))

                if len(ser_pri_trun) == 1:
                    #print('here3')
                    # We have a unique match
                    ser_matched[dt_det] = ser_pri_trun.index[0]
                    ser_not_in_primary = ser_not_in_primary.drop(dt_det)
                    break
                elif len(ser_pri_trun) == 2:
                    ####print('here4')
                    # We have 2 matches, if diff between windows is small enough
                    # these are equidistant from the detection. Choose the first.
                    ser_matched[dt_det] = ser_pri_trun.index[0]
                    ser_not_in_primary = ser_not_in_primary.drop(dt_det)
                    break
                elif len(ser_pri_trun) > 2:
                    ####print('here5')
                    # Too many matches, windows are 2 large.
                    break

    return ser_matched, ser_not_in_primary

def compare_to_HEK(df_all_hek, df_all_detections, windows=[timedelta(minutes=1), timedelta(minutes=2), timedelta(minutes=3)]):
    """

    Parameters
    ----------
    df_all_hek : ~`pandas.DataFrame`
        The dataframe containing HEK flare results for the same timerange and
        the detections, indexed by peaktimes.

        Note: HEK sometimes has multiple flares with the same peaktime, generally
        one of these will end up being ill-formed, for example have 0 duration.

    df_all_detections : ~`pandas.DataFrame`
        The pandas DataFrame of flare detections, that will be compared to the
        HEK (reference), should be indexed by flare peaktime and can optionally
        have start/end times.

    windows : `list`
        The list of windows to search in, these get progressively larger.
        Ideally you want the increase per window to be small enough to only add
        a single match.

    Returns
    -------
    df_out  : ~`pandas.DataFrame`
        The list of matches between the primary (reference) and the secondary list.
        The index is the index from the secondary, the values are the associated
        index from the primary (reference) list.

    ser_not_in_primary : ~`pandas.Series`
        The list of detected flares not in the HEK (reference) data.
        The index is the index from the detection data, the values are all NaN.

    ser_not_in_secondary : ~`pandas.Series`
        The list of flares in the HEK (reference) data and not in the detections.
        The index is the index from the HEK data, the values are all NaN.
    """
    # Does HEK data have any duplicates
    if not (np.all(np.invert(df_all_hek.index.duplicated()))):
        warnings.warn("HEK data has flares with dupicate peak times. First entries have been kept.")
        df_all_hek = df_all_hek.groupby(df_all_hek.index).first()

    # Find matches
    ser_matched, ser_not_in_primary, ser_not_in_secondary = compare_to_primary(df_all_hek['event_peaktime'], df_all_detections['event_peaktime'], windows=[timedelta(minutes=1), timedelta(minutes=2), timedelta(minutes=3)])

    # For matches, find values of event_peaktime, fl_peakflux, event_starttime, event_endtime, ... relative to HEK
    ####print('df_all_hek.index: \n'+str(df_all_hek.index))
    ####print('ser_matched.values: \n'+str(ser_matched.values))
    ####print('type(ser_matched.values): \n'+str(type(ser_matched.values)))
    #### The fllowing doesn't work, I'm doing a hack to get arround it
    #df_matched_hek = df_all_hek[ser_matched.values]
    #df_matched_det = df_all_detections[ser_matched.index]
    lis_indices = []
    for dt_index in ser_matched.values:
        lis_indices.append(df_all_hek.index.get_loc(dt_index))
    df_matched_hek = pd.DataFrame(df_all_hek.iloc[lis_indices])
    #print('\n\n')
    #print(df_matched_hek)
    #print('\n\n')
    lis_indices = []
    for dt_index in ser_matched.index:
        lis_indices.append(ser_matched.index.get_loc(dt_index))
    df_matched_det = pd.DataFrame(df_all_detections.iloc[lis_indices])
    arr_rel_event_peaktime = pd.to_datetime(ser_matched.index) - pd.to_datetime(ser_matched.values)
    ####print('arr_rel_event_peaktime:\n'+str(arr_rel_event_peaktime))

    # If we haven't calculated the HEK peak flux from the flare class then do so
    if not 'fl_peakflux' in df_all_hek.columns:
        df_matched_hek['fl_peakflux'] = utils.flare_class_to_intensity(df_matched_hek['fl_goescls'].values)
        #print('HEK fl_peakflux:\n'+str(df_matched_hek['fl_peakflux']))
    arr_rel_fl_peakflux = df_matched_det['fl_peakflux'].values / df_matched_hek['fl_peakflux'].values

    # Note: if the detected flares don't have start/end times determined, then return zero values. (would rather NaN's)
    if 'event_starttime' in df_all_detections.columns:
        arr_rel_event_starttime = pd.to_datetime(df_matched_det['event_starttime'].values) - pd.to_datetime(df_matched_hek['event_starttime'].values)
    else:
        arr_rel_event_starttime = np.zeros([len(df_matched_hek)])
    if 'event_endtime' in df_all_detections.columns:
        arr_rel_event_endtime = pd.to_datetime(df_matched_det['event_endtime'].values) - pd.to_datetime(df_matched_hek['event_endtime'].values)
    else:
        arr_rel_event_endtime = np.zeros([len(df_matched_hek)])

    # Make the dataframe for the relative values
    dic_matched_rel_data = { 'det_event_peaktime': ser_matched.index,
                             'hek_event_peaktime': ser_matched.values,
                             'fl_peakflux_rel_hek': arr_rel_fl_peakflux,
                             'event_peaktime_rel_hek': arr_rel_event_peaktime,
                             'event_peaktime_rel_hek_(s)': (arr_rel_event_peaktime/1000000000.0).astype(int),
                             'event_peaktime_rel_hek_(min)': (arr_rel_event_peaktime/60000000000.0).astype(int),
                             'event_starttime_rel_hek': arr_rel_event_starttime,
                             'event_starttime_rel_hek_(s)': (arr_rel_event_starttime/1000000000.0).astype(int),
                             'event_starttime_rel_hek_(min)': (arr_rel_event_starttime/60000000000.0).astype(int),
                             'event_endtime_rel_hek': arr_rel_event_endtime,
                             'event_endtime_rel_hek_(s)': (arr_rel_event_endtime/1000000000.0).astype(int),
                             'event_endtime_rel_hek_(min)': (arr_rel_event_endtime/60000000000.0).astype(int)
                             }
    df_out = pd.DataFrame(data=dic_matched_rel_data, index=ser_matched.index)

    return df_out, ser_not_in_primary, ser_not_in_secondary

def compare_many_to_primary(data_ref, dic_datasets, windows=[timedelta(minutes=1), timedelta(minutes=2), timedelta(minutes=3)]):
    """
    Basic function to compare multiple datasets to one primary (reference).

    ####This isn't complete yet!
    """
    # Sanatise inputs
    ser_ref = data_ref
    if isinstance(ser_ref, pd.DataFrame):
        ser_ref = pd.Series(data=np.nan, index=ser_ref.index)
    elif isinstance(ser_ref, np.ndarray):
        ser_ref = pd.Series(data=np.nan, index=pd.to_datetime(ser_ref))

    # Check all datasets
    lis_ser_data = []
    for dataset in datasets:
        ser_dataset = dataset
        if isinstance(ser_dataset, pd.DataFrame):
            ser_dataset = pd.Series(data=np.nan, index=ser_data.index)
        elif isinstance(ser_dataset, np.ndarray):
            ser_dataset = pd.Series(data=np.nan, index=pd.to_datetime(ser_data))
        lis_ser_data.append(ser_dataset)

    # Holders for matched/unmatched flares
    dic_temp = {}

    ser_matched = pd.Series()
    ser_unmatched = pd.Series(data=np.nan, index=ser_data.index)

    # For each dataset.
    for dataset in datasets:
        # For each of the values in the primary (reference) data
        for i in range(0, len(ser_data)):
            dt_det = ser_data.index[i]

            # Check if the flare detection is the same as a HEK entry
            if dt_det in ser_ref.index:
                ser_matched[dt_det] = dt_det
                ser_unmatched = ser_unmatched.drop(dt_det)
            else:
                # If not exactly the same, check within the given windows
                for td_window in windows:
                    dt_window_min = dt_det - td_window
                    dt_window_max = dt_det + td_window

                    # truncate to this window and see if any HEK results are present
                    ser_ref_trun = ser_ref.truncate(dt_window_min, dt_window_max)
                    if len(ser_ref_trun) == 1:
                        # We have a unique match
                        ser_matched[dt_det] = ser_ref_trun.index[0]
                        ser_unmatched = ser_unmatched.drop(dt_det)
                        break
                    elif len(ser_ref_trun) == 2:
                        # We have 2 matches, if diff between windows is small enough
                        # these are equidistant from the detection. Choose the first.
                        ser_matched[dt_det] = ser_ref_trun.index[0]
                        ser_unmatched = ser_unmatched.drop(dt_det)
                        break
                    elif len(ser_ref_trun) > 2:
                        # Too many matches, windows are 2 large.
                        break
    # Return the results
    return df_matched, ser_unmatched

max_stars_repo_name: Alex-Ian-Hamilton/flarepy

max_issues_repo_name: Alex-Ian-Hamilton/flarepy

max_forks_repo_name: Alex-Ian-Hamilton/flarepy


text:
"""Utils to check the samplers and compatibility with scikit-learn"""

# Adapated from scikit-learn
# Authors: Guillaume Lemaitre <g.lemaitre58@gmail.com>
# License: MIT

from __future__ import division

import sys
import traceback

from collections import Counter

import pytest

import numpy as np
from scipy import sparse
from pytest import raises

from sklearn.datasets import make_classification
from sklearn.cluster import KMeans
from sklearn.preprocessing import label_binarize
from sklearn.utils.estimator_checks import check_estimator \
    as sklearn_check_estimator, check_parameters_default_constructible
from sklearn.exceptions import NotFittedError
from sklearn.utils.testing import assert_allclose
from sklearn.utils.testing import set_random_state
from sklearn.utils.multiclass import type_of_target

from imblearn.over_sampling.base import BaseOverSampler
from imblearn.under_sampling.base import BaseCleaningSampler, BaseUnderSampler
from imblearn.ensemble.base import BaseEnsembleSampler
from imblearn.over_sampling import SMOTE
from imblearn.under_sampling import NearMiss, ClusterCentroids

from imblearn.utils.testing import warns

def _yield_sampler_checks(name, Estimator):
    yield check_target_type
    yield check_samplers_one_label
    yield check_samplers_no_fit_error
    yield check_samplers_X_consistancy_sample
    yield check_samplers_fit
    yield check_samplers_fit_sample
    yield check_samplers_ratio_fit_sample
    yield check_samplers_sparse
    yield check_samplers_pandas
    yield check_samplers_multiclass_ova

def _yield_all_checks(name, estimator):
    # trigger our checks if this is a SamplerMixin
    if hasattr(estimator, 'sample'):
        for check in _yield_sampler_checks(name, estimator):
            yield check

def check_estimator(Estimator):
    """Check if estimator adheres to scikit-learn conventions and
    imbalanced-learn

    This estimator will run an extensive test-suite for input validation,
    shapes, etc.
    Additional tests samplers if the Estimator inherits from the corresponding
    mixin from imblearn.base

    Parameters
    ----------
    Estimator : class
        Class to check. Estimator is a class object (not an instance).
    """
    name = Estimator.__name__
    # test scikit-learn compatibility
    sklearn_check_estimator(Estimator)
    check_parameters_default_constructible(name, Estimator)
    for check in _yield_all_checks(name, Estimator):
        check(name, Estimator)

def check_target_type(name, Estimator):
    X = np.random.random((20, 2))
    y = np.linspace(0, 1, 20)
    estimator = Estimator()
    set_random_state(estimator)
    with warns(UserWarning, match='should be of types'):
        estimator.fit(X, y)

def check_samplers_one_label(name, Sampler):
    error_string_fit = "Sampler can't balance when only one class is present."
    sampler = Sampler()
    X = np.random.random((20, 2))
    y = np.zeros(20)
    try:
        sampler.fit(X, y)
    except ValueError as e:
        if 'class' not in repr(e):
            print(error_string_fit, Sampler, e)
            traceback.print_exc(file=sys.stdout)
            raise e
        else:
            return
    except Exception as exc:
        print(error_string_fit, traceback, exc)
        traceback.print_exc(file=sys.stdout)
        raise exc

def check_samplers_no_fit_error(name, Sampler):
    sampler = Sampler()
    X = np.random.random((20, 2))
    y = np.array([1] * 5 + [0] * 15)
    with raises(NotFittedError, match="instance is not fitted yet."):
        sampler.sample(X, y)

def check_samplers_X_consistancy_sample(name, Sampler):
    sampler = Sampler()
    X = np.random.random((30, 2))
    y = np.array([1] * 20 + [0] * 10)
    sampler.fit(X, y)
    X_different = np.random.random((40, 2))
    y_different = y = np.array([1] * 25 + [0] * 15)
    with raises(RuntimeError, match="X and y need to be same array earlier"):
        sampler.sample(X_different, y_different)

def check_samplers_fit(name, Sampler):
    sampler = Sampler()
    X = np.random.random((30, 2))
    y = np.array([1] * 20 + [0] * 10)
    sampler.fit(X, y)
    assert hasattr(sampler, 'ratio_')

def check_samplers_fit_sample(name, Sampler):
    sampler = Sampler()
    X, y = make_classification(n_samples=1000, n_classes=3,
                               n_informative=4, weights=[0.2, 0.3, 0.5],
                               random_state=0)
    target_stats = Counter(y)
    X_res, y_res = sampler.fit_sample(X, y)
    if isinstance(sampler, BaseOverSampler):
        target_stats_res = Counter(y_res)
        n_samples = max(target_stats.values())
        assert all(value >= n_samples for value in Counter(y_res).values())
    elif isinstance(sampler, BaseUnderSampler):
        n_samples = min(target_stats.values())
        assert all(value == n_samples for value in Counter(y_res).values())
    elif isinstance(sampler, BaseCleaningSampler):
        target_stats_res = Counter(y_res)
        class_minority = min(target_stats, key=target_stats.get)
        assert all(target_stats[class_sample] > target_stats_res[class_sample]
                   for class_sample in target_stats.keys()
                   if class_sample != class_minority)
    elif isinstance(sampler, BaseEnsembleSampler):
        y_ensemble = y_res[0]
        n_samples = min(target_stats.values())
        assert all(value == n_samples
                   for value in Counter(y_ensemble).values())

def check_samplers_ratio_fit_sample(name, Sampler):
    # in this test we will force all samplers to not change the class 1
    X, y = make_classification(n_samples=1000, n_classes=3,
                               n_informative=4, weights=[0.2, 0.3, 0.5],
                               random_state=0)
    sampler = Sampler()
    expected_stat = Counter(y)[1]
    if isinstance(sampler, BaseOverSampler):
        ratio = {2: 498, 0: 498}
        sampler.set_params(ratio=ratio)
        X_res, y_res = sampler.fit_sample(X, y)
        assert Counter(y_res)[1] == expected_stat
    elif isinstance(sampler, BaseUnderSampler):
        ratio = {2: 201, 0: 201}
        sampler.set_params(ratio=ratio)
        X_res, y_res = sampler.fit_sample(X, y)
        assert Counter(y_res)[1] == expected_stat
    elif isinstance(sampler, BaseCleaningSampler):
        ratio = {2: 201, 0: 201}
        sampler.set_params(ratio=ratio)
        X_res, y_res = sampler.fit_sample(X, y)
        assert Counter(y_res)[1] == expected_stat
    elif isinstance(sampler, BaseEnsembleSampler):
        ratio = {2: 201, 0: 201}
        sampler.set_params(ratio=ratio)
        X_res, y_res = sampler.fit_sample(X, y)
        y_ensemble = y_res[0]
        assert Counter(y_ensemble)[1] == expected_stat

def check_samplers_sparse(name, Sampler):
    # check that sparse matrices can be passed through the sampler leading to
    # the same results than dense
    X, y = make_classification(n_samples=1000, n_classes=3,
                               n_informative=4, weights=[0.2, 0.3, 0.5],
                               random_state=0)
    X_sparse = sparse.csr_matrix(X)
    if isinstance(Sampler(), SMOTE):
        samplers = [Sampler(random_state=0, kind=kind)
                    for kind in ('regular', 'borderline1',
                                 'borderline2', 'svm')]
    elif isinstance(Sampler(), NearMiss):
        samplers = [Sampler(version=version)
                    for version in (1, 2, 3)]
    elif isinstance(Sampler(), ClusterCentroids):
        # set KMeans to full since it support sparse and dense
        samplers = [Sampler(random_state=0,
                            voting='soft',
                            estimator=KMeans(random_state=1,
                                             algorithm='full'))]
    else:
        samplers = [Sampler()]

    for sampler in samplers:
        set_random_state(sampler)
        X_res_sparse, y_res_sparse = sampler.fit_sample(X_sparse, y)
        X_res, y_res = sampler.fit_sample(X, y)
        if not isinstance(sampler, BaseEnsembleSampler):
                assert sparse.issparse(X_res_sparse)
                assert_allclose(X_res_sparse.A, X_res)
                assert_allclose(y_res_sparse, y_res)
        else:
            for x_sp, x, y_sp, y in zip(X_res_sparse, X_res,
                                        y_res_sparse, y_res):
                assert sparse.issparse(x_sp)
                assert_allclose(x_sp.A, x)
                assert_allclose(y_sp, y)

def check_samplers_pandas(name, Sampler):
    pd = pytest.importorskip("pandas")
    # Check that the samplers handle pandas dataframe and pandas series
    X, y = make_classification(n_samples=1000, n_classes=3,
                               n_informative=4, weights=[0.2, 0.3, 0.5],
                               random_state=0)
    X_pd, y_pd = pd.DataFrame(X), pd.Series(y)
    sampler = Sampler()
    if isinstance(Sampler(), SMOTE):
        samplers = [Sampler(random_state=0, kind=kind)
                    for kind in ('regular', 'borderline1',
                                 'borderline2', 'svm')]

    elif isinstance(Sampler(), NearMiss):
        samplers = [Sampler(version=version)
                    for version in (1, 2, 3)]

    else:
        samplers = [Sampler()]

    for sampler in samplers:
        set_random_state(sampler)
        X_res_pd, y_res_pd = sampler.fit_sample(X_pd, y_pd)
        X_res, y_res = sampler.fit_sample(X, y)
        assert_allclose(X_res_pd, X_res)
        assert_allclose(y_res_pd, y_res)

def check_samplers_multiclass_ova(name, Sampler):
    # Check that multiclass target lead to the same results than OVA encoding
    X, y = make_classification(n_samples=1000, n_classes=3,
                               n_informative=4, weights=[0.2, 0.3, 0.5],
                               random_state=0)
    y_ova = label_binarize(y, np.unique(y))
    sampler = Sampler()
    set_random_state(sampler)
    X_res, y_res = sampler.fit_sample(X, y)
    X_res_ova, y_res_ova = sampler.fit_sample(X, y_ova)
    assert_allclose(X_res, X_res_ova)
    if issubclass(Sampler, BaseEnsembleSampler):
        for batch_y, batch_y_ova in zip(y_res, y_res_ova):
            assert type_of_target(batch_y_ova) == type_of_target(y_ova)
            assert_allclose(batch_y, batch_y_ova.argmax(axis=1))
    else:
        assert type_of_target(y_res_ova) == type_of_target(y_ova)
        assert_allclose(y_res, y_res_ova.argmax(axis=1))

max_stars_repo_name: fengtongtong/imbalanced-learning

max_issues_repo_name: AnnYN/imbalanced-learn

max_forks_repo_name: AnnYN/imbalanced-learn


text:
import os
import sys
import torch
import numpy as np
import uuid

from torchvision import transforms, utils, datasets
from torch import nn, optim
from torch.autograd import Variable
from domain.model import VAE

setattr(sys.modules['__main__'], 'VAE', VAE)
sx = 128
imgsize = sx*sx*3
nusers = 100

def sigmoid(x):
    x = x.data.numpy()
    return 1 / (1 + np.exp(-x))

def generateRatings(z):
    if type(z) is list:
        z = Variable(torch.from_numpy(np.array(z)).float())
    mu_usrs = torch.load('domain/useremb.pt', map_location='cpu')
    pred = torch.matmul(mu_usrs, z)
    pred = sigmoid(pred)
    return pred

def decode(z):
    vae = torch.load('domain/vae_ratings_bestmodel.pt', map_location='cpu')
    if type(z) is list:
        z = Variable(torch.from_numpy(np.array(z)).float())
    fixed_z = z.view(1, 256)
    reconst_z = vae.decoder(fixed_z)
    sx = 128
    reconst_images = reconst_z.view(reconst_z.size(0), 3, sx, sx)

    # build path to serve image
    img_name = str(uuid.uuid4())[:8]
    img_url_prefix = 'output/' + img_name
    local_path_prefix = './' + img_url_prefix

    utils.save_image(reconst_images, local_path_prefix + '_128.png')
    return img_url_prefix + '_128.png'

max_stars_repo_name: crslab/cloud-vision-explorer

max_issues_repo_name: crslab/cloud-vision-explorer

max_forks_repo_name: crslab/cloud-vision-explorer


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Feb 28 15:05:15 2019

@author: dberke

Test script for conversions.py and the associated functions.
"""

import numpy as np
import pytest

import unyt as u

from varconlib.conversions import (air_indexEdlen53, vac2airESO, air2vacESO,
                                   vac2airMorton00, air2vacMortonIAU)

@pytest.fixture(scope='module')
def vacuum_array():
    return u.unyt_array(range(3800, 6900, 10), units=u.angstrom)

@pytest.fixture(scope='module')
def air_array_2D():
    return u.unyt_array(range(5000, 6000, 10), units=u.angstrom).reshape(10,
                                                                         10)

@pytest.fixture(scope='module')
def air_array():
    return u.unyt_array(range(3800, 6900, 10), units=u.angstrom)

@pytest.fixture(scope='module')
def unitless_array():
    return np.array([5000, 5001, 5002, 5003, 5004])

class TestEdlen1953(object):

    def testAirIndex(self):
        assert air_indexEdlen53(5000) == pytest.approx(1.0002789636500335)

    def testAirUnits(self, vacuum_array):
        assert vac2airESO(vacuum_array).units == u.angstrom

    def test2DVacuumUnits(self, air_array_2D):
        assert air2vacESO(air_array_2D).units == u.angstrom

    def testReshapeArrays(self, air_array_2D, verbose=True):
        assert air_array_2D.ndim == 2
        assert np.shape(air2vacESO(air_array_2D)) == (10, 10)

    def testVacuumUnits(self, air_array):
        assert air2vacESO(air_array, verbose=True).units == u.angstrom

    def testNoUnits(self, unitless_array):

        assert vac2airESO(unitless_array) == pytest.approx([4998.60557074,
                                                            4999.60530505,
                                                            5000.60503934,
                                                            5001.60477364,
                                                            5002.60450793])

class TestMorton2000(object):

    def testVac2Air(self, unitless_array):
        assert vac2airMorton00(unitless_array) == pytest.approx([4998.605522,
                                                                 4999.605256,
                                                                 5000.604990,
                                                                 5001.604724,
                                                                 5002.604459])

    def testAir2Vac(self, unitless_array):
        assert air2vacMortonIAU(unitless_array) == pytest.approx([5001.394848,
                                                                  5002.395114,
                                                                  5003.395380,
                                                                  5004.395646,
                                                                  5005.395911])

max_stars_repo_name: DBerke/varconlib

max_issues_repo_name: DBerke/varconlib

max_forks_repo_name: DBerke/varconlib


text:
import os, cv2
import numpy as np
import hdf5storage as h5io

from heatmap_overlay import *

def visual_color_img(RootDir, MethodNames, with_fix=0):

	imgsDir = RootDir + 'images/'
	fixsDir = RootDir + 'fixations/maps/'
	salsDir = RootDir + 'Results/Saliency/'

	img_ext = '.jpg'
	sal_ext = '.png'

	for idx_m in range(len(MethodNames)):
		print("---" + str(idx_m + 1) + "/" + str(len(MethodNames)) + "---: " + MethodNames[idx_m])

		salmap_dir = salsDir + MethodNames[idx_m] + '/'
		out_path = salmap_dir + 'Visual_color/'
		if not os.path.exists(out_path):
			os.makedirs(out_path)

		sal_names = [f for f in os.listdir(salmap_dir) if f.endswith(sal_ext)]
		sal_names.sort()

		for idx_n in range(len(sal_names)):
			print(str(idx_n + 1) + "/" + str(len(sal_names)) + ": " + sal_names[idx_n])

			file_name = sal_names[idx_n][:-4]
			outname = out_path + file_name + sal_ext
			if os.path.exists(outname):
				continue

			img    = cv2.imread(imgsDir + file_name + img_ext,-1)
			salmap = cv2.imread(salmap_dir + file_name + sal_ext,-1)

			fixname = fixsDir + file_name + '.mat'
			if with_fix and os.path.exists(fixname):
				fixmap = h5io.loadmat(fixname)["I"]

			overmap = heatmap_overlay(img,salmap)
			if with_fix and os.path.exists(fixname):
				fixpts_dilate = cv2.dilate(fixmap, np.ones((5, 5), np.uint8))
				fixpts_dilate = np.repeat(np.expand_dims(fixpts_dilate, axis=2), 3, axis=2)
				overmap[fixpts_dilate > 0.5] = 1

			overmap = overmap / np.max(overmap) *255
			cv2.imwrite(out_path + file_name + sal_ext, im2uint8(overmap))

if __name__ == "__main__":

	DataSet = 'salicon'
	RootDir = 'E:/Code/IIP_Saliency_Video/DataSet/salicon/val/'
	MethodNames = ['zk-twos-st']

	WITH_FIX = 0
	visual_color_img(RootDir, MethodNames, with_fix=WITH_FIX)