text:
# For data manipulation
import csv
import numpy as np
import pandas as pd
# For creating SVR model
from sklearn.svm import SVR
# For creating SVM model
from sklearn import svm, preprocessing
from sklearn.model_selection import cross_validate
from sklearn import metrics # For evaluation
from sklearn.model_selection import train_test_split
# For plotting data
import matplotlib.pyplot as plt
from matplotlib.pylab import rcParams
# Import warnings filter
from warnings import simplefilter
# For writing up in excel
import xlsxwriter
import openpyxl
# For saving the model
import pickle
# etc
import datetime
from datetime import date, timedelta
import os, os.path
import re
import time

print('import success')

## FOR IGNORING FUTURE WARNINGS
simplefilter(action='ignore', category=FutureWarning)

### SETUP THE FIGURE SIZE FOR PLOTTING
rcParams['figure.figsize'] = 10,6 # width 8, height 6

### DEFINE TODAY'S DATE
today_date = date.today().strftime('%Y-%m-%d')
print('SCRIPT RUN ON ' + today_date)
print('================================\n')

########################################## CONFIGURATION REQUIRED ##########################################

### PERCENTAGE OF TRAINING AND TESTING DATA
train_percent = 0.8
test_percent = 0.2

### CSV FILEPATH
time_period_directory = '1_year/'
stock_name = 'Tesco'
csv_stock_name = 'TSCO.L.csv'
csv_stock_path = '../../datasets/UK Stock Data/' + time_period_directory

### SAVE TIME PERIOD IN NEW VARIABLE
x = re.match("(.+)_(.+)/", str(time_period_directory))
time_period = str(x.group(1)) + ' ' + str(x.group(2))

### DEFINE THE FEATURE
feature = 'Date'

### MODEL SPECIFICATIONS
model_name = 'SVR'
model_summary_directory = '../.././'
model_summary_filename = 'Model Summary.xlsx'

########################################################
########################################################

### PREPROCESS DATA FIRST | return -> dates, prices and the original dates
def preprocess_data():
    ################################################### DATA PREPROCESSING

    ## Read data
    df = pd.read_csv(csv_stock_path + csv_stock_name) # Tesco Data

    ## Remove null values
    df = df.dropna()
    df = df.reset_index()
    df = df.drop('index', axis=1)

    ## Sort by Date
    df = df.sort_values('Date')

    ## Set the index to Date
    df.reset_index(inplace=True)

    ## Get date values
    dates = []
    # for i in range(len(df)):
    #     dates.append(int(df['Date'][i].split('-')[2]))

    df['Date'] = pd.to_datetime(df.Date, format='%Y-%m-%d')
    original_dates = df['Date'] # Save the original dates for plotting
    df.set_index("Date", inplace=True)
    df = df.drop('index', axis=1)

    ## Put dataframe to new variable called 'data'
    data = df.copy()

    ## Preview the data
    print('---- PREVIEW INITIAL DATA ----')
    print('Head data:\n' + str(data.head()))
    print('Tail data:\n' + str(data.tail()))

    ## Initialize the data X and Y variable (x: date, y: close price)
    prices = data['Close']
    for i in range(len(data)):
        dates.append(i)

    print(dates)

    ## Plotting initial data
    plt.plot(prices)
    plt.title('{} Close Price Preview'.format(stock_name))
    plt.xlabel('Date')
    plt.ylabel('Close Price')
    plt.show()

    return dates, prices, original_dates

def split_data(data_1, data_2, data_3):

    ## Convert to a list
    dates_data = data_1 # Already in a list
    price_data = data_2.tolist()
    original_dates_data = data_3.tolist()

    ## Split the data into train data and test data
    split = int(len(price_data)*train_percent)
    x_train_data, x_test_data = dates_data[0:split], dates_data[split:]
    y_train_data, y_test_data = price_data[0:split], price_data[split:]
    dates_train_data, dates_test_data = original_dates_data[0:split], original_dates_data[split:]

    ## Convert to numpy date format
    dates_train_data = np.array(dates_train_data, dtype='datetime64[D]')
    dates_test_data = np.array(dates_test_data, dtype='datetime64[D]')

    return x_train_data, x_test_data, y_train_data, y_test_data, dates_train_data, dates_test_data

### FOR MAKING PREDICTION | return -> predicted price on specific date on linear, polynomial and RBF model
def svr_predict(dates, prices, original_dates):

    x_train, x_test, y_train, y_test, dates_train, dates_test = split_data(dates, prices, original_dates)

    # Preview the first five X Train, Y Train, X Test, Y Test data
    print('== PREVIEW TRAIN AND TEST DATA')
    print('---------------------------------')
    print('X TRAIN:\n' + str([x_train[x] for x in range(0,5)]))
    print('Y TRAIN:\n' + str([y_train[x] for x in range(0,5)]))
    print('---------------------------------')
    print('X TEST:\n' + str([x_test[x] for x in range(0,5)]))
    print('Y TEST:\n' + str([y_test[x] for x in range(0,5)]))

    ## Convert to 1xn dimension
    x_train = np.reshape(x_train, (len(x_train), 1))
    x_test = np.reshape(x_test, (len(x_test), 1))

    # x : Training vector {array-like, sparse matrix}, shape (n_samples, n_features)
    # y : Target vector relative to x array-like, shape (n_samples,)

    ################################################### MODEL DEVELOPMENT

    svm_linear_param = {
        'kernel': 'linear',
        'C': 1e3
    }

    svm_polynomial_param = {
        'kernel': 'poly',
        'C': 1e3,
        'degree': 2
    }

    svm_rbf_param = {
        'kernel': 'rbf',
        'C': 1e3,
        'gamma': 0.1
    }

    ## Define hyperparameters for the SVR model
    svr_linear = SVR(kernel=svm_linear_param['kernel'], C=svm_linear_param['C']) # For linear classification, Epsilon = 0.1
    svr_polynomial = SVR(kernel=svm_polynomial_param['kernel'], C=svm_polynomial_param['C'], degree=svm_polynomial_param['degree']) # Polynomial order = 2, Epsilon = 0.1
    svr_rbf = SVR(kernel=svm_rbf_param['kernel'], C=svm_rbf_param['C'], gamma=svm_rbf_param['gamma']) # For non-linear classification, Epsilon = 0.1, Gamma = 0.1, Degree = 3

    ## Fit training data (x = dates, y = prices)
    print('----------------------------')
    print('Training SVR linear model..')
    svr_linear.fit(x_train, y_train)
    print('Training SVR polynomial model..')
    svr_polynomial.fit(x_train, y_train)
    print('Training SVR RBF model..')
    svr_rbf.fit(x_train, y_train)

    y_pred_linear = svr_linear.predict(x_test)
    y_pred_polynomial = svr_polynomial.predict(x_test)
    y_pred_rbf = svr_rbf.predict(x_test)

    pkl_filename1 = 'svr linear'
    pkl_filename2 = 'svr polynomial'
    pkl_filename3 = 'svr rbf'
    with open(pkl_filename1, 'wb') as file:
        pickle.dump(svr_linear, file)
    with open(pkl_filename2, 'wb') as file:
        pickle.dump(svr_polynomial, file)
    with open(pkl_filename3, 'wb') as file:
        pickle.dump(svr_rbf, file)

    ################################################### MODEL EVALUATION

    ## Evaluate the model using several metrics (Mean Absolute Error, Mean Squared Error, Root Mean Squared Error)
    mae_linear = metrics.mean_absolute_error(y_test, y_pred_linear)
    mae_polynomial = metrics.mean_absolute_error(y_test, y_pred_polynomial)
    mae_rbf = metrics.mean_absolute_error(y_test, y_pred_rbf)
    rmse_linear = np.sqrt(np.mean(np.power((np.array(y_test)-np.array(y_pred_linear)), 2)))
    rmse_polynomial = np.sqrt(np.mean(np.power((np.array(y_test)-np.array(y_pred_polynomial)), 2)))
    rmse_rbf = np.sqrt(np.mean(np.power((np.array(y_test)-np.array(y_pred_rbf)), 2)))

    print('-------------------------------------------')
    print('=== SVR Linear ===')
    print('MAE: ' + str(mae_linear))
    print('RMSE: ' + str(rmse_linear))
    print('=== SVR Polynomial ===')
    print('MAE: ' + str(mae_polynomial))
    print('RMSE: ' + str(rmse_polynomial))
    print('=== SVR Radial Basis Function ===')
    print('MAE: ' + str(mae_rbf))
    print('RMSE: ' + str(rmse_rbf))
    print('===========================================')

    ################################################### MODEL VISUALIZATION

    ## Compare the results between models for further observation
    result_linear = pd.DataFrame({'Actual': y_test, 'Predicted': y_pred_linear})
    result_polynomial = pd.DataFrame({'Actual': y_test, 'Predicted': y_pred_polynomial})
    result_rbf = pd.DataFrame({'Actual': y_test, 'Predicted': y_pred_rbf})
    print('== RESULTS COMPARISON ON TESTING DATA ==')
    print(result_linear)
    print('--------------------------')
    print(result_polynomial)
    print('--------------------------')
    print(result_rbf)

    ## Plotting results made by SVR
    print('== VISUALIZE RESULTS ==')
    plt.scatter(dates_test, y_test, color='orange', label='Actual Values')
    plt.plot(dates_test, y_pred_linear, color='green', label='Linear model')
    plt.plot(dates_test, y_pred_polynomial, color='blue', label='Polynomial model')
    plt.plot(dates_test, y_pred_rbf, color='red', label='RBF model')
    plt.xlabel('Date')
    plt.ylabel('Price')
    plt.title('Support Vector Regression | Predicted Values vs. Real Values')
    plt.legend()
    plt.show()

    ## Plotting the overall data
    plt.plot(dates_train, y_train, color='steelblue', label='Training Data')
    plt.plot(dates_test, y_test, color='cyan', label='Test Data')
    plt.plot(dates_test, y_pred_linear, c='green', label='SVR Linear result')
    plt.plot(dates_test, y_pred_polynomial, c='blue', label='SVR Polynomial result')
    plt.plot(dates_test, y_pred_rbf, c='red', label='SVR RBF result')
    plt.title('Support Vector Regression | {} Stock Prices Prediction'.format(stock_name))
    plt.xlabel('Dates')
    plt.ylabel('Prices')
    plt.legend()
    plt.show()

    save_model_summary(stock_name, time_period, feature, svm_linear_param, mae_linear, rmse_linear)
    save_model_summary(stock_name, time_period, feature, svm_polynomial_param, mae_polynomial, rmse_polynomial)
    save_model_summary(stock_name, time_period, feature, svm_rbf_param, mae_rbf, rmse_rbf)

## SAVING THE MODEL SUMMARY FOR MODEL COMPARISON
def save_model_summary(stock_name, time_period, feature, hyperparameter, mae, rmse):
    model_summary = {
        'A': stock_name, # Dataset
        'B': time_period, # Period
        'C': feature, # Features
        'D': hyperparameter, # Hyperparameters
        'E': '%.3f' % (mae), # MAE
        'F': '%.3f' % (rmse) # RMSE
    }

    ## Load workbook
    wb = openpyxl.load_workbook(model_summary_directory + model_summary_filename)

    ## Choose the sheet based on the model name
    sheet = wb.get_sheet_by_name(model_name)

    ## Start to write in a new line
    row_num = sheet.max_row + 1

    ## Write the model summary on specific row
    for key, value in model_summary.items():
        if isinstance(value, float):
            sheet[str(key) + str(row_num)].value = value
        else:
            sheet[str(key) + str(row_num)].value = str(value)

    wb.save(filename=model_summary_directory+model_summary_filename)
    print('Model summary successfully saved')

if __name__ == '__main__':
    DATE_VALUES, PRICE_VALUES, ORIGINAL_DATES = preprocess_data()
    svr_predict(DATE_VALUES, PRICE_VALUES, ORIGINAL_DATES)
    print('=================================================')


text:
'''
Created on Dec 15, 2016

@author: husensofteng
'''
import sys, os
import requests
import urllib
import numpy as np
import json

def get_cellnames_from_cellinfodict(cellinfodict_inputfile, cell_names_start_with="#"):
    cellinfo_lines = []
    with open(cellinfodict_inputfile, 'r') as cellinfodict_infile:
        cellinfo_lines = cellinfodict_infile.readlines()

    cell_lines = []
    cell_name = ""
    for line in cellinfo_lines:
        if line.startswith('//'):#skip lines start with // researved for comments in the config file
            continue
        elif line.startswith('***'):
            break
        elif line.startswith(cell_names_start_with):
            cell_name=line.strip().strip(cell_names_start_with)
            if cell_name not in cell_lines:
                cell_lines.append(cell_name)
    return cell_lines

def get_data_API(biosamples_out_dir="./", biosample_term_names_to_get=[]):
    HEADERS = {'accept': 'application/json'}
    biosample_term_names_to_get_str = '&biosample_term_name='+ '&biosample_term_name='.join(biosample_term_names_to_get)
    #Changed RNA-seq to total+RNA-seq in the URL
    URL = "https://www.encodeproject.org/search/?type=Experiment&assay_slims=Transcription&assay_title=total+RNA-seq&status=released&assembly=hg19&replicates.library.biosample.donor.organism.scientific_name=Homo+sapiens&files.file_type=tsv&files.analysis_step_version.analysis_step.pipelines.title=RNA-seq+of+long+RNAs+%28paired-end%2C+stranded%29&frame=object{}".format(biosample_term_names_to_get_str)
    response_json_dict = requests.get(URL, headers=HEADERS).json()
    #print json.dumps(response_json_dict, indent=4, separators=(',', ': '))
    bio_sample_files = {}
    '''
    iterate through each file in @graph list, for each file 'files' key send a request to get the file's info in json
    check file_type, assembly,... and the get link to the file in 'href'
    '''
    for cell in range(0, len(response_json_dict['@graph'])):#it contains a list, one element for each experiment
        bio_sample = response_json_dict['@graph'][cell]['biosample_term_name']
        if 'hg19' not in response_json_dict['@graph'][cell]['assembly'] or response_json_dict['@graph'][cell]['status']!='released':
            continue
        if bio_sample not in bio_sample_files.keys():
            bio_sample_files[bio_sample] = []
        for f in response_json_dict['@graph'][cell]['files']:
            f_info = requests.get("https://www.encodeproject.org/{}/".format(f), headers=HEADERS).json()
            if 'assembly' in f_info.keys() and f_info['file_type']=='tsv' and f_info['output_type']=='gene quantifications':
                if f_info['assembly'] == 'hg19':
                    if not os.path.exists(biosamples_out_dir+bio_sample):
                        os.mkdir(biosamples_out_dir+bio_sample)
                    print biosamples_out_dir+bio_sample + ':' + f_info['href']
                    bio_sample_files[bio_sample].append(f_info['href'])
                    if not os.path.exists(biosamples_out_dir+bio_sample+'/'+f_info['href'].split('/')[-1]):
                        urllib.urlretrieve('https://www.encodeproject.org/{}'.format(f_info['href']), biosamples_out_dir+bio_sample+'/'+f_info['href'].split('/')[-1])

    print "Finished downloading files for: ", bio_sample_files.keys()
    #now the content of bio_sample_files can be further processed to combine files of each bio_sample

def process_RNA_seq_datafolder(input_dir_path, num_header_lines=1, gene_id_index=0, gene_value_index=5,sep='\t'):
    bio_samples = {}
    bio_samples_output_dict = "bio_samples_dict.json"
    if os.path.exists(bio_samples_output_dict):
        with open(bio_samples_output_dict, 'r') as bio_samples_infile_dict:
            bio_samples = json.load(bio_samples_infile_dict)
        return bio_samples
    for bio_sample_input_dir in os.listdir(input_dir_path):
        if os.path.exists(input_dir_path+'/'+bio_sample_input_dir+'/'+bio_sample_input_dir+".bed"):
            continue
        bio_samples[bio_sample_input_dir] = {}
        for bio_sample_input_file in os.listdir(input_dir_path+'/'+bio_sample_input_dir):
            if bio_sample_input_file.split('.')[-1] != "tsv":
                continue
            with open(input_dir_path+'/'+bio_sample_input_dir+'/'+bio_sample_input_file, 'r') as bio_sample_infile:
                for i in range(0, num_header_lines):
                    bio_sample_infile.readline()#skip the header lines
                genes = bio_sample_infile.readlines()
                for g in genes:
                    if g=="" or g.startswith('//'):
                        continue
                    gene_id = g.strip().split(sep)[gene_id_index]
                    gene_value = float(g.strip().split(sep)[gene_value_index])
                    if gene_id not in bio_samples[bio_sample_input_dir].keys():
                        bio_samples[bio_sample_input_dir][gene_id] = []
                    bio_samples[bio_sample_input_dir][gene_id].append(gene_value)
    with open(bio_samples_output_dict, 'w') as bio_samples_output_dict_outfile:
        json.dump(bio_samples, bio_samples_output_dict_outfile)
    return bio_samples

def get_expr_per_bio_sample(bio_samples, gencode_id_info_dict, output_dir_path):

    genes_with_gencode_info_dict = {}
    for bio_sample in bio_samples:
        genes_with_gencode_info_dict[bio_sample] = {}

        if os.path.exists(output_dir_path+'/'+bio_sample + '/' + bio_sample+".bed"):
            return genes_with_gencode_info_dict

        with open(output_dir_path+'/'+bio_sample + '/' + bio_sample+".bed", 'w') as bio_sample_outfile:
            for gene_id in bio_samples[bio_sample]:
                if gene_id in gencode_id_info_dict.keys():
                    genes_with_gencode_info_dict[bio_sample][gene_id] = np.mean(bio_samples[bio_sample][gene_id])
                    bio_sample_outfile.write('\t'.join(gencode_id_info_dict[gene_id]) + '\t' + str(np.mean(bio_samples[bio_sample][gene_id])) + '\n')
    return genes_with_gencode_info_dict
    #combine results of all bio samples to one file (one col per bio sample) --- extend the the GTEx file

def get_gene_names_and_ids_from_genecode(genecode_genes_input_file, genecode_genes_only_genes_bed_output_file, num_header_lines=0, sep='\t'):
    gencode_id_info_dict = {}
    gene_id_name_dict = {}
    gencode_id_info_dict_savefile = genecode_genes_only_genes_bed_output_file+"saveddict"
    if os.path.exists(gencode_id_info_dict_savefile):
        with open(gencode_id_info_dict_savefile, 'r') as gencode_id_info_dict_file:
            return json.load(gencode_id_info_dict_file)
    with open(genecode_genes_input_file, 'r') as genecode_genes_infile, open(genecode_genes_only_genes_bed_output_file, 'w') as genecode_genes_only_genes_bed_outfile:
        for i in range(0, num_header_lines):
            genecode_genes_infile.readline()

        l = genecode_genes_infile.readline()
        while l:
            if l.startswith('#') or l=="":
                l = genecode_genes_infile.readline()
                continue
            sl = l.strip().split(sep)
            if sl[2]=="gene":
                sl_info_dict = dict([x.split('=') for x in sl[8].split(';')])
                gencode_id_info_dict[sl_info_dict['gene_id']] = [sl[0], sl[3], sl[4], sl[5], sl[6], sl_info_dict['gene_id'],
                                                              sl_info_dict['gene_name'], sl_info_dict['gene_status'], sl_info_dict['gene_type']]
                gene_id_name_dict[sl_info_dict['gene_id']] = sl_info_dict['gene_name']
                genecode_genes_only_genes_bed_outfile.write(sep.join(gencode_id_info_dict[sl_info_dict['gene_id']]) + '\n')
            l = genecode_genes_infile.readline()

    with open(gencode_id_info_dict_savefile, 'w') as gencode_id_info_dict_savefile_outfile:
        json.dump(gencode_id_info_dict, gencode_id_info_dict_savefile_outfile)
    return gencode_id_info_dict#, gene_id_name_dict

if __name__ == '__main__':

    biosamples_dir_path=sys.argv[1]
    if not os.path.exists(biosamples_dir_path):
        os.mkdir(biosamples_dir_path)
    cellinfodict_inputfile = sys.argv[2]
    biosample_term_names_to_get = get_cellnames_from_cellinfodict(cellinfodict_inputfile)
    #biosample_term_names_to_get = ['A549','GM12878','Gastric','HCT116','HEK293','HeLa-S3','HepG2','IMR-90','Ishikawa','K562','Kidney_curated','LNCaP clone FGC','MCF 10A','MCF-7','Ovary','Panc1','Pancreas_curated','Prostate_curated','SK-N-SH','T47D','U2OS','astrocyte','epithelial cell of esophagus','keratinocyte','osteoblast','urothelium cell line']
    genecode_genes_input_file=sys.argv[3]
    genecode_genes_only_genes_bed_output_file=genecode_genes_input_file+"_onlygenes.bed"

    get_data_API(biosamples_dir_path+'/', biosample_term_names_to_get)
    bio_samples = process_RNA_seq_datafolder(biosamples_dir_path)#give a dir that contains all bio_samples that have previously been downloaded
    gencode_id_info_dict = get_gene_names_and_ids_from_genecode(genecode_genes_input_file=genecode_genes_input_file, genecode_genes_only_genes_bed_output_file=genecode_genes_only_genes_bed_output_file)
    genes_with_gencode_info_dict = get_expr_per_bio_sample(bio_samples, gencode_id_info_dict, output_dir_path=biosamples_dir_path)

ext: py

reason: import numpy


text:
# Under MIT License, see LICENSE.txt

from math import ceil, sqrt
import numpy

from RULEngine.Util.constant import *
from RULEngine.Util.Position import Position

from ai.Algorithm.IntelligentModule import IntelligentModule

__author__ = 'RoboCupULaval'

class InfluenceMap(IntelligentModule):
    """
    Une classe représentant une influence map d'un terrain de robosoccer. C'est une matrice de cellule représentant le
    terrain auquelle on ajoute des points avec des influences qui se propagent sur le reste du terrain.

    La matrice et le code respecte le format suivant:
        - axe X ou axe 0 est l'axe représentant les rangées / l'axe 0 de la matrice /
          la largeur (axe Y) du terrain physique
        - axe Y ou axe 1 est l'axe représentant les colonnes / l'axe 1 de la matrice /
          la longueur (axe X) du terrain physique
        - Le point d'origine (0, 0) se situe au coin supérieur gauche du terrain physique.
        - Ce qui veut dire que c'est tableaux sont row-major.

    L'algorithm de propagation de l'influence est:
    (force appliqué à la case) = (force du point d'origine de l'influence) * ((facteur de réduction) ** (distance))
    transfomé en int arrondie vers 0.
    """

    def __init__(self, game_state, resolution=100, strength_decay=0.90, strength_peak=100, effect_radius=40,
                 have_static=False, have_it_executed=False):
        """
            Constructeur de la classe InfluenceMap

            :param game_state:  référence vers le game state
            :param resolution:    résolution des cases (défaut = 100)
            :param strength_decay: facteur de réduction de l'influence par la distance (défaut = 0.8)
            :param strength_peak:  maximum de la force appliquable par un point (est aussi le min) (défaut = 100)
            :param effect_radius:  distance qui borne la propagation de l'influence autour de l'origine (défaut = 40)
        """
        assert isinstance(resolution, int), "Creation InfluenceMap avec param resolution autre que int"
        assert isinstance(strength_decay, float), "Creation InfluenceMap avec param strength_decay autre que int"
        assert isinstance(effect_radius, int), "Creation InfluenceMap avec param effect_radius autre que int"
        assert isinstance(strength_peak, int), "Creation InfluenceMap avec param strength_peak autre que int"
        assert isinstance(have_static, bool), "Creation InfluenceMap avec param have_static autre que bool"
        assert isinstance(have_it_executed, bool), "Creation InfluenceMap avec param have_it_executed autre que bool"
        assert 0 < resolution, "Creation InfluenceMap avec param resolution <= 0"
        assert 0 < strength_decay < 1, "Creation InfluenceMap avec param strength_decay pas dans intervalle ]0, 1["
        assert 0 < strength_peak, "Creation InfluenceMap avec param strength_decay <= à 0"
        assert 0 < effect_radius, "Creation InfluenceMap avec param effect_radius <= 0"

        super().__init__(game_state)

        # board parameters
        self._resolution = resolution
        self._strength_decay = strength_decay
        self._strength_peak = strength_peak
        self._effect_radius = effect_radius
        self._border_strength = - strength_peak * 0.03  # TODO change this variable for something not out of thin air!

        # things on the baord parameters
        self._ball_position_on_board = ()

        self._number_of_rows, self._number_of_columns = self._calculate_rows_and_columns()

        if self.state is not None:
            self.have_it_executed = have_it_executed
            self._last_updated = self.state.get_timestamp()

        self._adjust_effect_radius()

        # different tableau pour enregistrer les différentes représentation
        self._static_boards = None
        self._starterboard = self._create_standard_influence_board()
        self._board = self._create_standard_influence_board()

        # those board are optionnal
        self._borders_board = self._create_standard_influence_board()
        self._goals_board = self._create_standard_influence_board()

        # todo determine how to choose if you want different kinds of static boards (ex: borders, goals, to determine..)
        if have_static:
            self._create_static_board()

    def update(self):
        if self.have_it_executed:
            if self.state.get_timestamp() - self._last_updated > 1:
                # purge the board with a new one (static or not)
                if self._static_boards is not None:
                    self._board = numpy.copy(self._static_boards)
                else:
                    self._board = self._create_standard_influence_board()

                self._update_and_draw_robot_position()
                self._update_ball_position()
                self.state.debug_manager.add_influence_map(self.export_board())
                self._last_updated = self.state.timestamp

    def export_board(self):
        return self._board.tolist()

    def find_points_over_strength_square(self, top_left_position, bottom_right_position, strength):
        """
        Retourne les points qui se trouve au-dessus ou égale à la force demandé dans le tableau principal(_board)

        :param Positon top_left_position: la rangé supérieure
        :param Position bottom_right_position: la rangé inférieure
        :param int strength: le force à trouvé qui est égale ou au dessus.

        :return: un liste de point qui se trouve au dessus ou égale à la force demandé
        """
        assert isinstance(top_left_position, Position), "Cette méthode requiert un object Position"
        assert isinstance(bottom_right_position, Position), "Cette méthode requiert un object Position"

        top_row, left_column = self._transform_field_to_board_position(top_left_position)
        bot_row, right_column = self._transform_field_to_board_position(bottom_right_position)
        ind_x, ind_y = numpy.where(self._board[top_row:bot_row, left_column:right_column] >= strength)
        ind_x = ind_x.tolist()
        ind_x = [x+top_row for x in ind_x]
        ind_y = ind_y.tolist()
        ind_y = [x+left_column for x in ind_y]
        indices = zip(ind_x, ind_y)
        return list(indices)

    def find_points_under_strength_square(self, top_left_position, bottom_right_position, strength):
        """
        Retourne les points qui se trouve au-dessous ou égale à la force demandé dans le tableau principal(_board)

        :param Positon top_left_position: la rangé supérieure
        :param Position bottom_right_position: la rangé inférieure
        :param int strength: le force à trouvé qui est égale ou au dessous.

        :return: un liste de point qui se trouve au-dessous ou égale à la force demandé
        :rtype: une liste de tuple rangée * colonne (int * int) list
        """
        assert isinstance(top_left_position, Position), "Cette méthode requiert un object Position"
        assert isinstance(bottom_right_position, Position), "Cette méthode requiert un object Position"

        top_row, left_column = self._transform_field_to_board_position(top_left_position)
        bot_row, right_column = self._transform_field_to_board_position(bottom_right_position)
        ind_x, ind_y = numpy.where(self._board[top_row:bot_row, left_column:right_column] <= strength)
        ind_x = ind_x.tolist()
        ind_x = [x+top_row for x in ind_x]
        ind_y = ind_y.tolist()
        ind_y = [x+left_column for x in ind_y]
        indices = zip(ind_x, ind_y)
        return list(indices)

    def find_max_value_in_board(self):
        """
        Permet de trouver les points du tableau qui ont la plus grande valeur du tableau

        :return: la valeur maximale du tableau et la liste de point (rangée * colonne) des point qui ont la valeur max
        :rtype: tuple (int * (int * int) list)
        """
        uniques = numpy.unique(self._board)
        max_in_board = uniques[-1]
        x, y = numpy.where(self._board >= max_in_board)
        x = x.tolist()
        y = y.tolist()
        indices = zip(x, y)
        return max_in_board, indices

    def find_min_value_in_board(self):
        """
        Permet de trouver les points du tableau qui ont la plus petite valeur du tableau

        :return: la valeur minimale du tableau et la liste de point (rangée * colonne) des point qui ont la valeur min
        :rtype: tuple (int * (int * int) list)
        """
        uniques = numpy.unique(self._board)
        min_in_board = uniques[0]
        x, y = numpy.where(self._board <= min_in_board)
        x = x.tolist()
        y = y.tolist()
        indices = zip(x, y)
        return min_in_board, indices

    def find_max_value_in_circle(self, center, radius):
        pass

    def get_influence_at_position(self, position):
        assert isinstance(position, Position), "accessing this function require a Position object"

        row_column_in_board = self._transform_field_to_board_position(position)
        return self._board[row_column_in_board[0], row_column_in_board[1]]

    def get_ball_influence(self):
        return self._board[self._ball_position_on_board[0], self._ball_position_on_board[1]]

    def get_number_of_cells(self):
        return self._number_of_rows * self._number_of_columns

    def print_board_to_file(self):
        """
        Create a file in the same running directory with a representation of the current board.
        """

        numpy.savetxt("IMBoard", self._board, fmt='%4i')
        # todo remove this line while not in debug mode
        print(self._starterboard.shape[0], " x ", self._starterboard.shape[1])

    def str(self):
        # todo comment and make sure this is right! Embelish?
        return "Influence Map - ", str(self._number_of_rows), " x ", str(self._number_of_columns)

    def _update_and_draw_robot_position(self):
        """
        Fetch la position des robots dans le gamestate et les applique sur le tableau principal
        """
        robots_position = []
        for i in range(self.state.get_count_player()):
            try:
                friend_position = self._transform_field_to_board_position(self.state.get_player_position(i))
                robots_position.append(friend_position)
                self._add_point_and_propagate_influence(friend_position[0], friend_position[1], self._board, 100)
            except:
                pass

            try:
                enemy_position = self._transform_field_to_board_position(self.state.get_enemy_position(i))
                robots_position.append(enemy_position)
                self._add_point_and_propagate_influence(enemy_position[0], enemy_position[1], self._board, -100)
            except:
                pass
        self._clamp_board(self._board)

    def _update_ball_position(self):
        self._ball_position_on_board = self._transform_field_to_board_position(self.state.get_ball_position())

    def _calculate_rows_and_columns(self):
        """
        Utilise la resolution pour calculer le nombre de rangée(rows) et de colonne(columns) pour le terrain

        :return: le nombre de rangées et le nombre de colonnes
        :rtype: tuple (int * int)
        """
        numberofrow = int(ceil((abs(FIELD_Y_BOTTOM) + FIELD_Y_TOP) / self._resolution))
        numberofcolumn = int(ceil((abs(FIELD_X_LEFT) + FIELD_X_RIGHT) / self._resolution))

        return numberofrow, numberofcolumn

    def _create_standard_influence_board(self):
        """
        Crée un objet numpy.ndarray, une liste à 2 dimenson de self._number_of_rows par self._number_of_columns d'int16.

        :return: Un numpy.ndarray d'int16 de self._number_of_rows par self._number_of_columns.
        :rtype: numpy.ndarray dtype=numpy.int16
        """
        return numpy.zeros((self._number_of_rows, self._number_of_columns), numpy.int16)

    def _adjust_effect_radius(self):
        """
        Permet d'ajuster le rayon d'effet de l'influence (self._effect_radius) si celui-ci éxède la propagation normale.

        Si la propagation de l'influence se retrouve à zéro avant d'avoir atteint le rayon d'effet défini à la création
        de la classe, cette méthode l'ajuste pour quelle donne exactement le rayon qui permet d'afficher tous les points
        qui ne donnent pas zéro.

        influence dans les cases:
        100 91 86 ... 6 5 4 3 2 1 0 0 0 0 0 0 0  sans effect_radius ajusté (si effect_radius est trop gros)
                v
        100 91 86 ... 6 5 4 3 2 1                avec effect_radius ajusté

        Coupe le nombre d'itérations des boucles de mise en place de l'influence lorsque le calcul donnerait une
        influence à 0.
        """
        distance_from_center = 0
        strenght_value_at_current_distance = self._strength_peak
        while strenght_value_at_current_distance >= 1:
            distance_from_center += 1
            decay = self._compute_value_by_distance(self._strength_peak, distance_from_center)
            strenght_value_at_current_distance = decay

        if self._effect_radius > distance_from_center:
            self._effect_radius = distance_from_center

    def _initialize_borders_board(self):
        """
        Ajoute des bordures sur le board dédie _borders_board
        """
        temp_array = self._create_standard_influence_board()
        self._put_and_propagate_borders(temp_array)
        numpy.add(temp_array, self._borders_board, out=self._borders_board)

    def _put_and_propagate_borders(self, board_to_apply):
        """
        Mets des bordures sur un array numpy.ndarray vierge.

        :param board_to_apply: Un numpy.ndarray un array vierge pour y ajouter des bordures.
        :type board_to_apply: numpy.ndarray dtype=numpy.int16
        """
        assert(isinstance(board_to_apply, numpy.ndarray))
        assert(board_to_apply.shape[0] == self._number_of_rows)
        assert(board_to_apply.shape[1] == self._number_of_columns)
        assert(board_to_apply.dtype == numpy.int16)

        board_to_apply[0] = self._border_strength
        board_to_apply[:, 0] = self._border_strength

        # keep the effectradius low while making the border speed up greatly the initialization, also you don't need
        # so much border
        # TODO see if this variable is okay, maybe change the way it is determined
        border_variance = 2
        temp_effectradius = int(ceil(ROBOT_RADIUS / self._resolution)) + border_variance

        # TODO make sure this does what you think it does.
        # Top border
        for border in range(0, self._number_of_columns):
            # only for the rows affected by the change.
            for x in range(0, temp_effectradius):

                columnmin = max(0, (border - temp_effectradius - 1))
                columnmax = min(self._number_of_columns, (border + temp_effectradius + 1))

                # for every columns affected
                for y in range(columnmin, columnmax):
                    decay = self._compute_value_by_distance(self._border_strength, self._distance(0, border, x, y))
                    board_to_apply[x, y] += decay

        # left border
        for border in range(0, self._number_of_rows):

            for y in range(0, temp_effectradius):

                rowmin = max(0, (border - temp_effectradius - 1))
                rowmax = min(self._number_of_rows, border + temp_effectradius + 1)

                for x in range(rowmin, rowmax):
                    decay = self._compute_value_by_distance(self._border_strength, self._distance(border, 0, x, y))
                    board_to_apply[x, y] += decay

        # Prend l'image créer et la flip l-r et u-d puis additionne
        temp_inverse_board = numpy.copy(board_to_apply)
        temp_inverse_board = temp_inverse_board[::-1, ::-1]
        board_to_apply += temp_inverse_board

    def _initialize_goals_board(self, v_h_goal_offset):
        """
        Mets des buts sur le starterboard et le board
        """
        # todo fetch which side we are on.
        temp_array = self._create_standard_influence_board()

        self._put_goals_and_propagate(v_h_goal_offset, temp_array)
        numpy.add(temp_array, self._goals_board, out=self._goals_board)

    def _calculate_goal_vertical_offset(self):
        """
        Calcule la dimension vertical du but p/r à la résolution.
        :return: la valeur verticale du but ajusté à la résolution
        :rtype: int
        """
        return int(ceil(FIELD_GOAL_Y_TOP / self._resolution))

    def _calculate_goal_horizontal_offset(self):
        """
        Calcule la dimension horizontale du but p/r à la résolution.

        :return: la valeur horizontale du but ajusté à la résolution
        :rtype: int
        """
        return int(ceil(FIELD_GOAL_SEGMENT / self._resolution))

    def _put_goals_and_propagate(self, v_h_offset, board_to_apply):
        """
        Mets des buts sur le numpy.ndarray fourni.

        :param v_h_offset: la dimension verticale et horizontale des buts.
        :type v_h_offset: tuple (int * int)
        :param board_to_apply: Un numpy.ndarray un array vierge pour y ajouter des buts.
        :type board_to_apply: numpy.ndarray dtype=numpy.int16
        """

        # TODO take into account what team you are ie: orientation and strength adjustment
        for x in range(int(self._number_of_rows / 2 - v_h_offset[0]),
                       int(self._number_of_rows / 2 + v_h_offset[0]) + 1,
                       6):
            self._add_point_and_propagate_influence(x, 0, board_to_apply, 100)

        numpy.add((numpy.negative(board_to_apply[:, ::-1])), board_to_apply, out=board_to_apply)

    def _create_static_board(self):
        """
        Crée des tableaux pour les bordures et les buts et les enregistres dans les variables de classes déjà présente.

        Il reste à déterminer comment on s'occupe de choisir quel tableau static à créer.
        """
        # TODO see to make flags to control what goes into the static board
        self._static_boards = self._create_standard_influence_board()

        self._initialize_borders_board()
        numpy.add(self._borders_board, self._static_boards, out=self._static_boards)

        v_h_goal_offset = (self._calculate_goal_vertical_offset(), self._calculate_goal_horizontal_offset())
        self._initialize_goals_board(v_h_goal_offset)
        numpy.add(self._goals_board, self._static_boards, out=self._static_boards)
        self._clamp_board(self._static_boards)
        # numpy.savetxt("Debug", self._static_boards, fmt='%5i')

    def _compute_value_by_distance(self, strength, distance):
        """
        Calcule la valeur qu'une recoit d'un point de force strength à distance distance.

        :param strength: la force du point central
        :type strength: int
        :param distance: la distance entre le point central et le point où calculer la value
        :type distance: int or float

        :return: la valeur calculé
        :rtype: int

        """
        return int(strength * (self._strength_decay ** distance))

    def _clamp_board(self, board_to_clamp):
        """
        Arrondis toutes les cellules du tableau pour qu'ils soient dans [-self._strength_peak, self._strength_peak].

        :param board_to_clamp: Un numpy.ndarray à clampé
        :type board_to_clamp: numpy.ndarray dtype=numpy.int16
        """
        numpy.clip(board_to_clamp, -self._strength_peak, self._strength_peak, out=board_to_clamp)

    def _add_point_and_propagate_influence(self, row, column, board_to_apply, strength=0):
        """
        Pose un point et propage son influence sur le tableau donné.

        Il est important de clamper le board après avoir utiliser cette méthode!

        :param int row: la rangée du point d'origine de la force
        :param int column: int La colonne du point d'origine de la force
        :param board_to_apply: numpy.ndarray l'array sur lequel on ajoute un point et on le propage
        :type board_to_apply: numpy.ndarray dtype=numpy.int16
        :param int strength: la force du point à appliquer
        """
        assert (isinstance(row, int))
        assert (isinstance(column, int))
        assert (0 <= row < self._number_of_rows)
        assert (0 <= column < self._number_of_columns)
        assert (isinstance(board_to_apply, numpy.ndarray))
        assert (board_to_apply.shape[0] == self._number_of_rows)
        assert (board_to_apply.shape[1] == self._number_of_columns)
        assert (isinstance(strength, int))
        assert (-self._strength_peak <= strength <= self._strength_peak)

        rowmin = max(0, (row - self._effect_radius))
        rowmax = min((self._number_of_rows - 1), (row + self._effect_radius))

        columnmin = max(0, (column - self._effect_radius))
        columnmax = min((column + self._effect_radius), self._number_of_columns)

        cases_iterator = numpy.nditer(board_to_apply[rowmin:(rowmax + 1), columnmin:(columnmax + 1)],
                                      flags=['multi_index'], op_flags=['readwrite'])

        while not cases_iterator.finished:
            to_put = self._compute_value_by_distance(strength, self._distance(cases_iterator.multi_index[0],
                                                                              cases_iterator.multi_index[1],
                                                                              (row - rowmin),
                                                                              (column - columnmin)))
            influence_already_in_case = cases_iterator[0]
            influence_to_put_instead = to_put + influence_already_in_case
            cases_iterator[0] = influence_to_put_instead
            cases_iterator.iternext()

    def _transform_field_to_board_position(self, position):
        assert(isinstance(position, Position))
        assert(FIELD_X_LEFT <= position.x <= FIELD_X_RIGHT)
        assert(FIELD_Y_BOTTOM <= position.y <= FIELD_Y_TOP)
        # this should hold up

        xpos = -position.x + ((abs(FIELD_X_LEFT) + FIELD_X_RIGHT) / 2)
        ypos = position.y + ((abs(FIELD_Y_BOTTOM) + FIELD_Y_TOP) / 2)

        xpos = int(xpos / self._resolution)
        ypos = int(ypos / self._resolution)

        if ypos == self._number_of_rows:
            ypos -= 1
        if xpos == self._number_of_columns:
            xpos -= 1

        return ypos, xpos

    def _transform_board_to_field_position(self, row, column):
        assert(isinstance(row, int))
        assert(isinstance(column, int))
        assert(0 <= row <= self._number_of_rows)
        assert(0 <= column <= self._number_of_columns)

        # This should hold up
        ypos = row * self._resolution
        xpos = column * self._resolution

        ypos = (ypos - ((abs(FIELD_Y_BOTTOM) + FIELD_Y_TOP) / 2)) + (self._resolution / 2)
        xpos = (xpos - ((abs(FIELD_X_LEFT) + FIELD_X_RIGHT) / 2)) + (self._resolution / 2)

        tempposition = Position(xpos, ypos)
        return tempposition

    def _distance(self, x1, y1, x2, y2):
        assert (isinstance(x1, int or float))
        assert (isinstance(x2, int))
        assert (isinstance(y1, int))
        assert (isinstance(y2, int))

        return sqrt(((x2 - x1) ** 2) + ((y2 - y1) ** 2))


text:
# ==BEGIN LICENSE==
#
# MIT License
#
# Copyright (c) 2018 SRI Lab, ETH Zurich
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
#
# ==END LICENSE==

import numpy as np
import numpy.random as r

from dpfinder.algorithms.algorithms import Algorithm, precision_np, precision_tf
from dpfinder.utils.utils import n_hot

class Sum(Algorithm):

	def __init__(self, array_size):
		super().__init__()
		self.array_size = array_size
		self.epsilon = 0.1

	def get_random_input(self):
		a_init = np.random.uniform(-1, 1, size=self.array_size)
		d_init = np.random.uniform(-1, 1, size=1)
		a_init = a_init.astype(np.float64)
		d_init = d_init.astype(np.float64)
		return a_init, d_init

	def get_random_output(self, a, b):
		expected = np.sum(a)
		o_init = np.empty(2, dtype=np.float64)
		o_init[0] = expected - 3
		o_init[1] = expected + 3
		return o_init

max_stars_repo_name: barryZZJ/dp-finder

max_issues_repo_name: barryZZJ/dp-finder

max_forks_repo_name: barryZZJ/dp-finder


text:
import numpy as np
import pandas as pd

def test_1():
    sr = pd.Series([1, 2, 3])
    assert np.ndarray == type(sr.values) and [
        1, 2, 3] == sr.values.tolist(), 'series value error'

    assert np.ndarray == type(sr.index.values) and [
        0, 1, 2] == sr.index.values.tolist(), 'series index error'

def test_2():
    # copy 跟asarray类似, 都是只有参数是ndarray,并且不改变dtype的时候,才会不复制.
    # 1. 复制: 参数为array不会改
    a1 = [1, 2, 3]
    sr = pd.Series(a1, copy=False)
    sr.iloc[0] = 100
    assert [1, 2, 3] == a1 and [
        100, 2, 3] == sr.values.tolist(), 'series copy error'
    # 2. 不复制: 参数为ndarray会改
    a2 = np.array([1, 2, 3])
    sr = pd.Series(a2, copy=False)
    sr.iloc[0] = 100
    assert [100, 2, 3] == a2.tolist() and [
        100, 2, 3] == sr.values.tolist(), 'series copy 2 error'
    # 3. 复制: 参数为ndarray,但是修改dtype 也不会改
    a3 = np.array([1, 2, 3])
    assert np.dtype('i4') == a3.dtype, 'series copy 3 error'
    sr = pd.Series(a3, dtype=np.dtype('i8'), copy=False)
    sr.iloc[0] = 100
    assert [1, 2, 3] == a3.tolist() and [
        100, 2, 3] == sr.values.tolist(), 'series copy 3 error'

def test_3():
    # index
    a1 = pd.Series(['xiao', 'xinmiao'], index=['x', 'y'])
    assert 'xinmiao' == a1['y'], 'series index error'

    a2 = pd.Series({'x': 10, 'y': 24})
    assert 24 == a2['y'], 'series index 2 error'

max_stars_repo_name: relax-space/pandas-first

max_issues_repo_name: relax-space/pandas-first

max_forks_repo_name: relax-space/pandas-first


text:
import gym
import numpy as np
from models.q_table import Q_Table
from observation_wrapper import ObservationWrapper

# todo convert all of this to pytorch so we can use autograd :)
class QAgent():
    def __init__(self, obs_size: int, action_size: int, alpha: float, gamma: float, epsilon: float):
        self.obs_size = obs_size
        self.action_size = action_size

        self.model = Q_Table(self.action_size, alpha, gamma, epsilon)

    def step(self, env: gym.Env, state: np.array):
        action = self.model.get_action(state, env.uniform_random_action())

        while not env.valid_moves()[action]:
            action = env.uniform_random_action() # force random action if action is illegal
            if isinstance(action, np.ndarray): # break if action is an array( i.e. trading environment, there are no invalid actions)
                break

        next_state, reward, done, info = env.step(action)

        self.model.update_q(state, action, reward, next_state)

        return next_state, reward, done, info

max_stars_repo_name: sirmammingtonham/IMT-Agent

max_issues_repo_name: sirmammingtonham/IMT-Agent

max_forks_repo_name: sirmammingtonham/IMT-Agent


text:
"""
Objects containing methods for quantifying Sit to stand transitions based on accelerometer data

Lukas Adamowicz
July 2019
Pfizer
"""
from numpy import ceil, log2, arange, abs, sum, sqrt, diff
from numpy.linalg import norm
from scipy.fftpack import fft
from sit2standpy.utility import Transition

class TransitionQuantifier:
    """
    Quantification of a sit-to-stand transition.
    """
    def __init__(self):
        pass

    def quantify(self, times, fs, raw_acc=None, mag_acc_f=None, mag_acc_r=None, v_vel=None, v_pos=None):
        """
        Compute quantitative values from the provided signals

        Parameters
        ----------
        times : tuple
            Tuple of the start and end timestamps for the transition
        mag_acc_f : {None, numpy.ndarray}, optional
            Filtered acceleration magnitude during the transition.
        mag_acc_r : {None, numpy.ndarray}, optional
            Reconstructed acceleration magnitude during the transition.
        v_vel : {None, numpy.ndarray}, optional
            Vertical velocity during the transition.
        v_pos : {None, numpy.ndarray}, optional
            Vertical position during the transition.

        Returns
        -------
        transition : Transition
            Transition object containing metrics quantifying the transition.
        """
        if raw_acc is not None:
            acc = norm(raw_acc, axis=1)
            sparc, _, _ = TransitionQuantifier.sparc(acc, fs)
        else:
            sparc = None
        if mag_acc_f is not None:
            max_acc = mag_acc_f.max()
            min_acc = mag_acc_f.min()
        else:
            max_acc = None
            min_acc = None
        if v_vel is not None:
            max_v_vel = v_vel.max()
            min_v_vel = v_vel.min()
        else:
            max_v_vel = None
            min_v_vel = None
        if v_pos is not None:
            v_disp = v_pos[-1] - v_pos[0]
        else:
            v_disp = None

        self.transition_ = Transition(times=times, v_displacement=v_disp, max_v_velocity=max_v_vel,
                                      min_v_velocity=min_v_vel, max_acceleration=max_acc, min_acceleration=min_acc,
                                      sparc=sparc)

        return self.transition_

    @staticmethod
    def sparc(x, fs, padlevel=4, fc=10.0, amp_th=0.05):
        """
        SPectral ARC length metric for quantifying smoothness

        Parameters
        ----------
        x : numpy.ndarray
            Array containing the data to be analyzed for smoothness
        fs : float
            Sampling frequency
        padlevel : int, optional
            Indicates the amount of zero-padding to be done to the movement data for estimating the spectral arc length.
            Default is 4.
        fc : float, optional
            The max cutoff frequency for calculating the spectral arc length metric. Default is 10.0
        amp_th : float, optional
            The amplitude threshold to be used for determining the cut off frequency up to which the spectral arc length
            is to be estimated. Default is 0.05

        Returns
        -------
        sal : float
            The spectral arc length estimate of the given data's smoothness
        (f, Mf) : (numpy.ndarray, numpy.ndarray)
            The frequency and the magnitude spectrum of the input data. This spectral is from 0 to the Nyquist frequency
        (f_sel, Mf_sel) : (numpy.ndarray, numpy.ndarray)
            The portion of the spectrum that is selected for calculating the spectral arc length

        References
        ----------
        S. Balasubramanian, A. Melendez-Calderon, A. Roby-Brami, E. Burdet. "On the analysis of movement smoothness."
        Journal of NeuroEngineering and Rehabilitation. 2015.
        """
        # number of zeros to be padded
        nfft = int(pow(2, ceil(log2(len(x))) + padlevel))

        # frequency
        f = arange(0, fs, fs / nfft)
        # normalized magnitude spectrum
        Mf = abs(fft(x, nfft))
        Mf = Mf / Mf.max()

        # indices to choose only the spectrum withing the given cutoff frequency Fc
        # NOTE: this is a low pass filtering operation to get rid of high frequency noise from affecting the next step
        # (amplitude threshold based cutoff for arc length calculation
        fc_inx = ((f <= fc) * 1).nonzero()
        f_sel = f[fc_inx]
        Mf_sel = Mf[fc_inx]

        # choose the amplitude threshold based cutoff frequency. Index of the last point on the magnitude spectrum is
        # greater than or equal to the amplitude threshold
        inx = ((Mf_sel >= amp_th) * 1).nonzero()[0]
        fc_inx = range(inx[0], inx[-1] + 1)
        f_sel = f_sel[fc_inx]
        Mf_sel = Mf_sel[fc_inx]

        # calculate the arc length
        sal = -sum(sqrt(pow(diff(f_sel) / (f_sel[-1] - f_sel[0]), 2) + pow(diff(Mf_sel), 2)))

        return sal, (f, Mf), (f_sel, Mf_sel)

reason: from numpy,from scipy


text:
#Author: Shadab Alam, March 2015
#This program takes the pair count and normalization from
#PairCount's code output and compute 2d correlation and
#monopole and quadrupole moment of corrlation

import numpy as np
import os
import sys

#import pylab as pl
#import matplotlib
#matplotlib.use('Agg')
#import matplotlib.pyplot as pl

f_write=[]
f_plot=[]

def load_PairCount_output(root):
   #load sbins
   lines=open(root+'-DD.dat').readlines()
   sbins=np.array([np.float(x) for x in lines[0].split()])
   ns=sbins.size-1
   #load mu bins
   mubins=np.array([np.float(x) for x in lines[1].split()])
   nmu=mubins.size-1
   #load pair counts
   DD=np.loadtxt(root+'-DD.dat',skiprows=2)
   DR=np.loadtxt(root+'-DR.dat',skiprows=2)
   try:
      RR=np.loadtxt(root+'-RR.dat',skiprows=2)
   except:
      RR=np.loadtxt(root[:-4]+'0001-RR.dat',skiprows=2)
   #load wieghts
   lines=open(root+'-norm.dat').readlines()
   Wsum_data=np.float(lines[0].split()[-1])
   Wsum_rand=np.float(lines[1].split()[-1])

   Wrat=Wsum_data/Wsum_rand
   DR=DR*Wrat
   RR=RR*Wrat*Wrat
   #print ns,nmu,DD.shape, Wsum_data, Wsum_rand
   #print sbins, mubins

   #combined all output in a dictionary and return
   pcdict={'sbins': sbins, 'ns': ns, 'mubins':mubins, 'nmu':nmu,
	 'DD': DD, 'DR': DR, 'RR': RR ,
	 'Wsum_data': Wsum_data, 'Wsum_rand': Wsum_rand, 'Wrat':Wrat}
   return pcdict

def load_PairCount_output_cross(root):
   #load sbins
   lines=open(root+'-D1D2.dat').readlines()
   sbins=np.array([np.float(x) for x in lines[0].split()])
   ns=sbins.size-1
   #load mu bins
   mubins=np.array([np.float(x) for x in lines[1].split()])
   nmu=mubins.size-1
   #load pair counts
   D1D2=np.loadtxt(root+'-D1D2.dat',skiprows=2)
   D1R2=np.loadtxt(root+'-D1R2.dat',skiprows=2)
   R1D2=np.loadtxt(root+'-R1D2.dat',skiprows=2)
   R1R2=np.loadtxt(root+'-R1R2.dat',skiprows=2)

   #load wieghts
   lines=open(root+'-norm.dat').readlines()
   Wsum_data1=np.float(lines[0].split()[-1])
   Wsum_rand1=np.float(lines[1].split()[-1])
   Wsum_data2=np.float(lines[2].split()[-1])
   Wsum_rand2=np.float(lines[3].split()[-1])

   Wrat1=Wsum_data1/Wsum_rand1
   Wrat2=Wsum_data2/Wsum_rand2

   #Normalize the pair counts
   D1R2=D1R2*Wrat2
   R1D2=R1D2*Wrat1
   R1R2=R1R2*Wrat1*Wrat2

   #combined all output in a dictionary and return
   pcdict={'sbins': sbins, 'ns': ns, 'mubins':mubins, 'nmu':nmu,
	 'D1D2': D1D2, 'D1R2': D1R2, 'R1D2': R1D2,'R1R2': R1R2 }

   return pcdict

def combine_count(root1,root2):
  #load the count for each of the root
  sbins1, ns1, mubins1, nmu1, DD1,DR1,RR1, Wsum_data1,Wsum_rand1=load_PairCount_output(root1)
  sbins2, ns2, mubins2, nmu2, DD2,DR2,RR2, Wsum_data2,Wsum_rand2=load_PairCount_output(root2)

  Wrat=(Wsum_data1+Wsum_data2)/(Wsum_rand1+Wsum_rand2)
  #giving normalized pair counts
  DD=DD1+DD2
  DR=(DR1+DR2)*Wrat
  RR=(RR1+ RR2)*Wrat*Wrat

  if(np.sum(sbins1==sbins2)==sbins1.size and ns1==ns2
     and np.sum(mubins1==mubins2)==mubins1.size and nmu1==nmu2):
     return sbins1, ns1, mubins1, nmu1, DD,DR,RR
  else:
     return 'error'

def compute_xi02(root,pcdict,xitype='auto', nscomb=1,write=0):
   if(write):
      outfile=root+'-xi2D.dat'
      f_write.append(outfile)
      fout=open(outfile,'w')

   if(xitype=='auto'):
      pclist=['DD', 'DR','RR']
   elif(xitype=='cross'):
      pclist=['D1D2', 'D1R2','R1D2','R1R2']

   #for rebinning
   nsrebins=np.int(np.ceil(pcdict['ns']/np.float(nscomb)))
   srebins=np.zeros(nsrebins)
   nmurebins=pcdict['nmu']
   murebins=np.zeros(nmurebins)

   xi2drebin=np.zeros(nsrebins*nmurebins).reshape(nsrebins,nmurebins)

   rebin=-1
   for ii in range(0,pcdict['ns'],nscomb):
      ss=np.mean(pcdict['sbins'][ii:ii+nscomb+1])
      rebin=rebin+1
      srebins[rebin]=ss
      for jj in range(0,nmurebins):
         mu=np.mean(pcdict['mubins'][jj:jj+2])
         if(rebin==0):
            murebins[jj]=mu

	 pcbin={}
	 for pctype in pclist:
	    pcbin[pctype]=np.sum(pcdict[pctype][ii:ii+nscomb,jj])

         if(xitype=='auto'):
            xi2d=(pcbin['DD']-2*pcbin['DR']+pcbin['RR'])/pcbin['RR']
	 elif(xitype=='cross'):
            xi2d=(pcbin['D1D2']-pcbin['D1R2']-pcbin['R1D2']+pcbin['R1R2'])/pcbin['R1R2']

         xi2drebin[rebin,jj]=xi2d
         if(write):
            fout.write('%12.8lf %12.8lf %12.8lf' %(ss,mu,xi2d))
            for pctype in pclist:
	       fout.write('%12.8lf '%pcbin[pctype])
	    fout.write('\n')

   if(write):
      fout.close()

   #print srebins, murebins
   return srebins, murebins, xi2drebin

#Assume sampling in rp-rpi computes xi2D and wp
def compute_xi2D_rppi(xi2droot,wproot,pcdict,xitype='auto',nscomb=1,write=0):
   if(write):
      outfile=xi2droot+'-xi2D.dat'
      f_write.append(outfile)
      fout=open(outfile,'w')

      outfile=wproot+'-wp.dat'
      f_write.append(outfile)
      fwp=open(outfile,'w')

   if(xitype=='auto'):
      pclist=['DD', 'DR','RR']
   elif(xitype=='cross'):
      pclist=['D1D2', 'D1R2','R1D2','R1R2']

   wp=np.zeros(pcdict['nper'])
   #for rebinning
   nper_rebins=pcdict['nper']
   npar_rebins=np.int(np.ceil(pcdict['npar']/np.float(nscomb)))
   rper_rebins=np.zeros(nper_rebins)
   rpar_rebins=np.zeros(npar_rebins)

   xi2drebin=np.zeros(nper_rebins*npar_rebins).reshape(nper_rebins,npar_rebins)

   for ii in range(0,pcdict['nper'],1):
      ss=np.mean(pcdict['rper'][ii:ii+2])
      rper_rebins[ii]=ss
      rebin=-1
      for jj in range(0,pcdict['npar'],nscomb):
         mu=np.mean(pcdict['rpar'][jj:jj+nscomb+1])
         rebin=rebin+1
         rpar_rebins[rebin]=mu

	 pcbin={}
	 for pctype in pclist:
	    pcbin[pctype]=np.sum(pcdict[pctype][ii,jj:jj+nscomb])

         if(xitype=='auto'):
            xi2d=(pcbin['DD']-2*pcbin['DR']+pcbin['RR'])/pcbin['RR']
	 elif(xitype=='cross'):
            xi2d=(pcbin['D1D2']-pcbin['D1R2']-pcbin['R1D2']+pcbin['R1R2'])/pcbin['R1R2']

         xi2drebin[ii,rebin]=xi2d
         if(write):
            fout.write('%12.8lf %12.8lf %12.8lf' %(ss,mu,xi2d))
            for pctype in pclist:
	       fout.write('%12.8lf '%pcbin[pctype])
	    fout.write('\n')

      dr_par=np.mean(rpar_rebins[1:]-rpar_rebins[:-1])
      wp[ii]=np.sum(xi2drebin[ii,:])*dr_par
      if(write):
         fwp.write('%12.8lf %12.8lf\n'%(ss,wp[ii]))

   if(write):
      fout.close()
      fwp.close()

   return rper_rebins, rpar_rebins, xi2drebin,wp
#end of compute wp and xi2d with rprpi sampling

#a function to call for pair count to xi2d and xi02 in polar co-ordinate
def xi2d_xi02_froot(root,xi2droot,xi02root,xitype='auto',nscomb=1,samp='rmu',write=0):
   #compute the correlation functions
   if(xitype=='auto'):
      pcdict=load_PairCount_output(root)
   elif(xitype=='cross'):
      pcdict=load_PairCount_output_cross(root)

   if(samp=='logr-theta'):
      pcdict['sbins']=np.power(10,pcdict['sbins'])

   srebins, murebins, xi2drebin=compute_xi02(xi2droot,pcdict,xitype=xitype,nscomb=nscomb,write=write)
   xi02=Xi_Legendre(xi02root,srebins, murebins, xi2drebin,samp=samp,write=write)

   return srebins,murebins,xi2drebin,xi02

#a function to call for pair count to xi2d and wp in cartesian co-ordinate
def xi2d_wp_froot(root,xi2droot,wproot,xitype='auto',nscomb=1,samp='log',write=0):
   #compute the correlation functions
   #rper,nper,rpar,npar, DD,DR,RR, Wsum_data,Wsum_rand=load_PairCount_output(root)
   if(xitype=='auto'):
      pcdict=load_PairCount_output(root)
   elif(xitype=='cross'):
      pcdict=load_PairCount_output_cross(root)

   #change dictionary key for cartesian coordinate
   keyreplace={'rper': 'sbins','nper': 'ns','rpar':'mubins','npar':'nmu'}
   for key in keyreplace.keys():
      pcdict[key]=pcdict.pop(keyreplace[key])

   if(samp=='logrp-pi'):
      pcdict['rper']=np.power(10,pcdict['rper'])

   rper_rebins,rpar_rebins,xi2drebin,wp=compute_xi2D_rppi(xi2droot,wproot,
                   pcdict,xitype=xitype,nscomb=nscomb,write=write)

   return rper_rebins, rpar_rebins, xi2drebin,wp

def xi2d_wp_xi02_froot_jn(inroot,xi2droot,outroot,xitype='auto',nscomb=1,samp='logrp-pi',compgZ=0,NJN=1):
   polar_coord=['rmu','rtheta','logr-theta']
   cart_coord=['rp-pi','logrp-pi']
   #files to write the compiled xi2d and wp functions from jackife
   xi2dfile=xi2droot+'-'+samp+'-NJN-%d.txt'%NJN

   for ii in range(0,NJN+1):
      if(NJN==0):
         root=inroot+'-'+samp
      elif(ii<NJN): #individual regions
         root=inroot+'-'+samp+'-JNdir/jk-%d'%ii
      else: #All data
         root=inroot+'-All-'+samp
      if(samp in polar_coord):
         srebins,murebins,xi2drebin,xi02=xi2d_xi02_froot(root,
	                    '','',xitype=xitype,samp=samp,nscomb=nscomb,write=0)
	 if(compgZ==1):
	    gZ=compute_gZ('',srebins, murebins, xi2drebin,samp=samp,write=0)

      elif(samp in cart_coord):
         rper_rebins, rpar_rebins, xi2drebin,wp=xi2d_wp_froot(root,
                           '','',xitype=xitype,samp=samp,nscomb=nscomb,write=0)
      else:
	 print 'Invalid sampling: %s'%samp
	 sys.exit()

      if(ii==0):
         n2d=xi2drebin.size
         xi2dmat=np.zeros(n2d*(NJN+1)).reshape(n2d,NJN+1)
         if(samp in polar_coord):
	    nxi02=xi02.shape[0]
	    rxi02=xi02[:,0]
	    xi0mat=np.zeros(nxi02*(NJN+1)).reshape(nxi02,NJN+1)
	    xi2mat=np.zeros(nxi02*(NJN+1)).reshape(nxi02,NJN+1)
	    xi1mat=np.zeros(nxi02*(NJN+1)).reshape(nxi02,NJN+1) #dipole usefule only for rtheta
	    if(compgZ==1):
	       gZmat=np.zeros(nxi02*(NJN+1)).reshape(nxi02,NJN+1)
	 else:
            nwp=rper_rebins.size
            wpmat = np.zeros(nwp*(NJN+1)).reshape(nwp,NJN+1)

      xi2dmat[:,ii]=xi2drebin.reshape(n2d)
      if(samp in polar_coord):
	 xi0mat[:,ii]=xi02[:,1]
	 xi2mat[:,ii]=xi02[:,2]
	 xi1mat[:,ii]=xi02[:,3]
	 if(compgZ==1):
	    gZmat[:,ii]=gZ[:,1]
      else:
         wpmat[:,ii]=wp

   if(samp in polar_coord):
      xi2daxes=[srebins,murebins]
      col1=rxi02
      outtags=['xi0','xi2']
      outmats=[xi0mat,xi2mat]
      if(samp=='rmu'):
         xi2dtag=['r','mu']
      elif(samp=='rtheta' or samp=='logr-theta'):
         xi2dtag=['r','theta']
         outmats.append(xi1mat)
	 outtags.append('xi1')
      else:
	 print 'Invalid sampling: %s'%samp

      if(compgZ==1):
	 outtags.append('gZ')
         outmats.append(gZmat)

   else:
      xi2daxes=[rper_rebins,rpar_rebins]
      xi2dtag=['rper','rpar']
      outmats=[wpmat]
      col1=rper_rebins
      outtags=['wp']

   #compute the mean and std err for xi2d and write to file
   meanxi2d=np.mean(xi2dmat[:,:NJN],axis=1)
   errxi2d=np.sqrt(NJN-1)*np.std(xi2dmat[:,:NJN],axis=1)
   xi2dwrite=np.column_stack([meanxi2d,errxi2d,
                                   xi2dmat[:,NJN],xi2dmat[:,:NJN]])
   with file(xi2dfile,'w') as fxi2d:
      fxi2d.write('#xi2d with jacknife, root:%s NJN=%d\n'%(root,NJN))
      for ww,arr in enumerate(xi2daxes):
         w_str='#%s: '%xi2dtag[ww]
         for rr in arr:
            w_str="%s %10.5f"%(w_str,rr)
         fxi2d.write(w_str+'\n')

      if(NJN>0):
         fxi2d.write('#meanxi2d sigmaxi2d Allxi2d jacknifecolumns\n')
         np.savetxt(fxi2d,xi2dwrite,fmt='% 18.8lf')
      elif(NJN==0):
         fxi2d.write('#xi2d\n')
	 np.savetxt(fxi2d,xi2dmat[:,NJN],fmt='% 18.8lf')

   print 'written: ',xi2dfile

   for pole, outmat in enumerate(outmats):
     #compute the mean and std err for wp and write to file
      meanout=np.mean(outmat[:,:NJN],axis=1)
      errout=np.sqrt(NJN-1)*np.std(outmat[:,:NJN],axis=1)
      outwrite=np.column_stack([col1,meanout,errout,outmat[:,NJN],outmat[:,:NJN]])
      outfile=outroot+'-'+outtags[pole]+'-'+samp+'-NJN-%d.txt'%NJN
      with file(outfile,'w') as fout:
         fout.write('#%s with jacknife, root:%s NJN=%d\n'%(outtags[pole],root,NJN))
	 if(NJN>0):
            fout.write('#r(Mpc/h) mean sigma All jacknifecolumns\n')
            np.savetxt(fout,outwrite,fmt='% 18.8lf')
	 elif(NJN==0):
            fout.write('#r(Mpc/h) val\n')
	    np.savetxt(fout,np.column_stack([col1,outmat[:,NJN]]),fmt='% 18.8lf')
      print 'written:',outfile

   return 0

def compute_xi02_cross(root,sbins, ns, mubins, nmu,
        D1D2,D1R2,R1D2,R1R2, nscomb=1,write=0):
   if(write):
      outfile=root+'cross-xi2D.dat'
      f_write.append(outfile)
      fout=open(outfile,'w')

   #for rebinning
   nsrebins=np.int(np.ceil(ns/np.float(nscomb)))
   srebins=np.zeros(nsrebins)
   nmurebins=nmu
   murebins=np.zeros(nmu)

   xi2drebin=np.zeros(nsrebins*nmu).reshape(nsrebins,nmu)

   rebin=-1
   for ii in range(0,ns,nscomb):
      ss=np.mean(sbins[ii:ii+nscomb+1])
      rebin=rebin+1
      srebins[rebin]=ss
      for jj in range(0,nmu):
         mu=np.mean(mubins[jj:jj+2])
         if(rebin==0):
            murebins[jj]=mu
         D1D2bin=np.sum(D1D2[ii:ii+nscomb,jj])
         D1R2bin=np.sum(D1R2[ii:ii+nscomb,jj])
         R1D2bin=np.sum(R1D2[ii:ii+nscomb,jj])
         R1R2bin=np.sum(R1R2[ii:ii+nscomb,jj])

         xi2d=(D1D2bin-D1R2bin-R1D2bin+R1R2bin)/R1R2bin
         xi2drebin[rebin,jj]=xi2d
         if(write):
            fout.write('%12.8lf %12.8lf %12.8lf\n' %(ss,mu,xi2d))

   if(write):
      fout.close()

   #print srebins, murebins
   return srebins, murebins, xi2drebin
#End of compute xi02 cross correlation

def Xi_Legendre(root,srebins, murebins, xi2drebin,samp='rmu',write=1):
   if(samp=='rmu'):
      dmu=np.average(murebins[1:]-murebins[:-1])
      diffmu=1 #to account for the fact that integral is 0 to 1
   elif(samp=='rtheta' or samp=='logr-theta'): #in this case murebin is actually theta and not mu
      dtheta=np.average(murebins[1:]-murebins[:-1])
      dmu=np.sin(murebins)*dtheta
      murebins=np.cos(murebins) #this converts theta to mu for integration
      diffmu=2 #to account for the fact that itegral is -1 to 1

   mu2=np.power(murebins,2)
   P0 =np.ones(mu2.size)/2.0    #(2l+1)Pl(mu)*sqrt(1-mu^2) monopole term
   P2=2.5*(3*mu2-1)/2.0         #p2=(3mu^2-1)/2  quadrupole term
   P1=3*murebins/2.0

   ns=srebins.size
   xi02=np.zeros(ns*4).reshape(ns,4)

   if(write):
      outfile=root+'-xi02.dat'
      f_write.append(outfile)
      fout=open(outfile,'w')

   for ii in range(0,ns):
      xi02[ii,0]=srebins[ii]
      xi02[ii,1]=np.sum(xi2drebin[ii,:]*P0*dmu)*2.0/diffmu
      xi02[ii,2]=np.sum(xi2drebin[ii,:]*P2*dmu)*2.0/diffmu
      xi02[ii,3]=np.sum(xi2drebin[ii,:]*P1*dmu)*2.0/diffmu
      if(write==1 and smap=='rmu'):
         fout.write('%12.8lf %12.8lf %12.8lf\n' %(xi02[ii,0],xi02[ii,1],xi02[ii,2]))
      if(write==1 and smap=='rtheta'):
         fout.write('%12.8lf %12.8lf %12.8lf %12.8lf\n' %(xi02[ii,0],xi02[ii,1],xi02[ii,2],xi02[ii,3]))

   #plot_xi02(xi02)
   return xi02

def compute_gZ(root,srebins, murebins, xi2drebin,samp='rtheta',write=1):
   '''This function computes the gravitational redshift from 2d correlation'''

   H=100 #hubble constant

   if(samp=='rtheta' or samp=='logr-theta'): #in this case murebin is actually theta and not mu
      dtheta=np.average(murebins[1:]-murebins[:-1])
      dmu=np.sin(murebins)*dtheta
      murebins=-np.cos(murebins) #this converts theta to mu for integration
      diffmu=2 #to account for the fact that itegral is -1 to 1
   else:
      print 'invalid sampling for gravitaional redshift: %s'%samp
      return 0

   ns=srebins.size
   gZ=np.zeros(ns*2).reshape(ns,2)
   numer=np.zeros(ns)
   denom=np.zeros(ns)

   if(write):
      outfile=root+'-gz.dat'
      f_write.append(outfile)
      fout=open(outfile,'w')

   #print np.min(murebins),np.max(murebins)
   for ii in range(0,ns):
      gZ[ii,0]=srebins[ii]
      #applying a mu limit
      mulim=1.0 #1.0
      ind1=murebins<mulim; ind2=murebins>-mulim; ind12=ind1*ind2
      #remove infinities and Nan
      indNaN=np.isnan(xi2drebin[ii,:])
      indinf=np.isinf(xi2drebin[ii,:])
      indcheck=indNaN+indinf
      #remove the infinity and its symmetric points
      indsym=np.copy(indcheck)

      for muinf in murebins[indcheck]:
         indsym=indsym+(np.abs(murebins+muinf)<1e-4)
      indcheck=indsym
      #print indcheck
      indall=~indcheck*ind12
      if(np.sum(indcheck)>0):
	 print 'infinities:(r,mu,xi)',srebins[ii],murebins[indcheck],xi2drebin[ii,indcheck]
	 #print ~indcheck*ind12
	 #import pylab as pl
	 #print indcheck
	 #pl.plot(murebins[~indcheck],xi2drebin[ii,~indcheck],'s-')
	 #pl.show()

      #summing
      #denom[ii]=np.sum((1+xi2drebin[ii,indall])*gZ[ii,0]*dtheta)
      #numer[ii]=H*np.sum((1+xi2drebin[ii,indall])*gZ[ii,0]*gZ[ii,0]*murebins[indall]*dtheta)
      #simpson rule for integration
      #Interpolation integration library
      import scipy.integrate as ssI
      from scipy import interpolate
      denom[ii]=ssI.simps(1+xi2drebin[ii,indall],murebins[indall])
      numer[ii]=H*ssI.simps((1+xi2drebin[ii,indall])*gZ[ii,0]*murebins[indall],murebins[indall])

      gZ[ii,1]=numer[ii]/denom[ii]

      if(write==1):
         fout.write('%12.8lf %12.8lf\n' %(gZ[ii,0],gZ[ii,1]))

   if(write==1):
      fout.close()

   if(0):#To integrate along r
      Rmin=1.0
      nR=15
      Rmax=60.0
      dLR=(np.log(Rmax)-np.log(Rmin))/nR
      #interpolate numer and denomr
      Idenomr = interpolate.splrep(srebins,denom, s=0,k=1)
      Inumerr = interpolate.splrep(srebins,numer, s=0,k=1)
      rr2d_tmp=np.linspace(0,Rmax,2000)
      denom_rtmp=interpolate.splev(rr2d_tmp,Idenomr, der=0)
      numer_rtmp=interpolate.splev(rr2d_tmp,Inumerr, der=0)

      rr=np.zeros(nR)
      gZnew=np.zeros(ns*2).reshape(ns,2)
      for ii in range(0,nR):
         r1=np.exp(np.log(Rmin)+ii*dLR)
         r2=np.exp(np.log(Rmin)+(ii+1)*dLR)
         rr[ii]=0.5*(r1+r2)
         ind1=rr2d_tmp>r1
         ind2=rr2d_tmp<r2
         denominator =ssI.simps(denom_rtmp[ind1*ind2],rr2d_tmp[ind1*ind2])
         numerator   =ssI.simps(numer_rtmp[ind1*ind2],rr2d_tmp[ind1*ind2])

         gZnew[ii,0]=rr[ii]
         gZnew[ii,1]=numerator/denominator
         if(0):
            import pylab as pl
            pl.figure(10)
            pl.plot(gZnew[ii,0],denominator,'b*')
            pl.plot(rr2d_tmp[ind1*ind2],denom_rtmp[ind1*ind2],'b-')
            pl.figure(20)
            pl.plot(gZnew[ii,0],numerator,'r*')
            pl.plot(rr2d_tmp[ind1*ind2],numer_rtmp[ind1*ind2],'r-')
            pl.figure(30)
            pl.plot(gZnew[ii,0],gZnew[ii,1],'k*')
            pl.plot(rr2d_tmp[ind1*ind2],numer_rtmp[ind1*ind2]/denom_rtmp[ind1*ind2],'k-')
            pl.show()
      #return np.column_stack([gZ,gZnew])
      #pl.plot(gZnew[:,0],gZnew[:,1])
      #pl.show()
      return gZnew

   return gZ
#end compute_gZ

def plot_xi02(xi02,mark='o-',lab='',lw=1):
   r=xi02[:,0]
   r2=r*r
   pl.plot(r,r2*xi02[:,1],mark,linewidth=lw,label=lab)
   pl.plot(r,r2*xi02[:,2],mark,linewidth=lw)
   pl.xlabel('r (Mpc/h)',fontsize=32)
   pl.ylabel(r'$r^2 \xi_l$',fontsize=32)
   if(lab!=''):
      pl.legend()

   return 0

def mean_cov(froot='',fold='',sky='NS',Njack=28,Sjack=9):
   xiarr=0
   reg=0
   if(sky=='NS'):
      jack=Njack+Sjack
   elif(sky=='N'):
      jack=Njack
   else:
      jack=Sjack

   title=''
   pl.figure()
   for s in sky:
      if s=='N':
         ext='-xi02.dat'
         nloop=Njack
         mark='r-'
         title=title+' N=red '
      elif s=='S':
         ext='-xi02.dat'
         nloop=Sjack
         mark='b-'
         title=title+' S=blue '
      for ii in range(0,nloop):
         if(fold=='../XI02/'):
            fname=froot+str(ii)+'-'+s+ext
         else:
            fname=froot+str(reg)+'-'+s+str(ii)+ext
         print fname
         data=np.loadtxt(fold+fname)
         if(ii==0 and xiarr==0):
            ns=data.shape[0]
            rr=data[:,0]
            rr2=rr*rr
            xi02=np.zeros(ns*2*jack).reshape(2*ns,jack)
            xiarr=1
         xi02[:ns,reg]=data[:,1]
         xi02[ns:,reg]=data[:,2]
         pl.plot(rr,rr2*data[:,1],mark)
         pl.plot(rr,rr2*data[:,2],mark)
         reg=reg+1

   #mean and error bar
   mean_xi=np.mean(xi02,axis=1)
   cov=np.cov(xi02)*jack
   xierr=np.sqrt(np.diagonal(cov))
   pl.errorbar(rr,rr2*mean_xi[:ns],yerr=rr2*xierr[:ns],fmt='o',color='k',markersize=10,label='mean')
   pl.errorbar(rr,rr2*mean_xi[ns:],yerr=rr2*xierr[ns:],fmt='o',color='k',markersize=10)

   #plot full sample
   data=np.loadtxt(fold+'All-'+sky+'-xi02.dat')
   pl.plot(rr,rr2*data[:,1],'m*-',linewidth=4,markersize=10,label='full sample')
   pl.plot(rr,rr2*data[:,2],'m*-',linewidth=4,markersize=10)
   pl.xlabel(r'$r Mpc/h $',fontsize=32)
   pl.ylabel(r'$r^2 \xi_l$',fontsize=32)
   pl.title(title,fontsize=32)
   pl.legend()
   pl.tight_layout()
   outplot='plots/'+froot+'-'+sky+'.png'
   f_plot.append(outplot)
   pl.savefig(outplot)

   pl.figure()
   pl.pcolor(cov)
   pl.colorbar()

   corr=np.copy(cov)
   for ii in range(0,cov.shape[0]):
      for jj in range(0,cov.shape[1]):
         corr[ii,jj]=cov[ii,jj]/np.sqrt(cov[ii,ii]*cov[jj,jj])

   pl.figure()
   pl.pcolor(corr)
   pl.colorbar()
   outplot='plots/'+froot+'-corr.png'
   f_plot.append(outplot)
   pl.savefig(outplot)

   return 0

def Mocks_xi_comb(rootN='N',rootS='S',nmock1=0,nmock2=1000,outfold='xi2d',plots_dir='plots',nscomb=1,write=0):

   #compute north only and south only and north+south correlations
   pl.figure(1)  #for xi02 each region
   pl.figure(2)  #for xi02 combined for each region
   for ii in range(nmock1,nmock2):
      if(ii%50==0):
         print ii
      #compute the correlation functions for north
      root1=rootN+'-'+str(ii).zfill(4)
      outroot=outfold+'ngc/XI2D/'+root1.split('/')[-1]
      sbins, ns, mubins, nmu, DD,DR,RR, Wsum_data,Wsum_rand=load_PairCount_output(root1)
      Wrat=Wsum_data/Wsum_rand
      srebins, murebins, xi2drebin=compute_xi02(outroot,sbins, ns, mubins, nmu,
                                  DD,DR,RR, Wrat,nscomb=nscomb,write=write)
      xi02=Xi_Legendre(root1,srebins, murebins, xi2drebin,outfold=outfold+'ngc/XI02/',write=write)
      pl.figure(1)
      plot_xi02(xi02,mark='r-')

      #compute the correlation function for south
      root2=rootS+'-'+str(ii).zfill(4)
      outroot=outfold+'sgc/XI2D/'+root2.split('/')[-1]
      sbins, ns, mubins, nmu, DD,DR,RR, Wsum_data,Wsum_rand=load_PairCount_output(root2)
      Wrat=Wsum_data/Wsum_rand
      srebins, murebins, xi2drebin=compute_xi02(outroot,sbins, ns, mubins, nmu,
                                  DD,DR,RR, Wrat,nscomb=nscomb,write=write)
      xi02=Xi_Legendre(root2,srebins, murebins, xi2drebin,outfold=outfold+'sgc/XI02/',write=write)
      pl.figure(2)
      plot_xi02(xi02,mark='b-')

      #compute the combined correlation function
      sbins, ns, mubins, nmu, DD,DR,RR=combine_count(root1,root2)
      outroot=outfold+'combgc/XI2D/'+root1.split('/')[-1]
      srebins, murebins, xi2drebin=compute_xi02(outroot,sbins, ns, mubins, nmu,
                                  DD,DR,RR, 1.0,nscomb=nscomb,write=write)
      xi02=Xi_Legendre(root1,srebins, murebins, xi2drebin,outfold=outfold+'combgc/XI02/',write=write)
      pl.figure(3)
      plot_xi02(xi02,mark='k-')

   pl.figure(1)
   pl.title(rootN.split('/')[-1])
   pl.savefig(plots_dir+rootN.split('/')[-1]+'-ngc.png')

   pl.figure(2)
   pl.title(rootS.split('/')[-1])
   pl.savefig(plots_dir+rootS.split('/')[-1]+'-sgc.png')

   pl.figure(3)
   pl.title(rootN.split('/')[-1])
   pl.savefig(plots_dir+rootN.split('/')[-1]+'-comb.png')

   return 0
#End of Mocks

def AllJN_realize(nscomb=8):
   rootN='../CMASS-N-JN/'
   rootS='../CMASS-S-JN/'
   outfold1='../XI02/'
   outfoldcomb='../XI02-NS/'

   njack_N=28
   njack_S=9
   reg=0

   #compute north only and north + souht correlations
   root2=rootS+'All-S'
   pl.figure(1)  #for xi02 each region
   pl.figure(2)  #for xi02 combined for each region
   for ii in range(0,njack_N):
      root1=rootN+'JN'+str(ii)+'-N'
      #compute the correlation functions
      sbins, ns, mubins, nmu, DD,DR,RR, Wsum_data,Wsum_rand=load_PairCount_output(root1)
      Wrat=Wsum_data/Wsum_rand
      srebins, murebins, xi2drebin=compute_xi02(root1,sbins, ns, mubins, nmu,
                                  DD,DR,RR, Wrat,nscomb=nscomb,write=1)
      xi02=Xi_Legendre(root1,srebins, murebins, xi2drebin,outfold=outfold1,write=1)
      pl.figure(3)
      plot_xi02(xi02,mark='ro-',lab='N-'+str(ii))
      pl.figure(1)
      plot_xi02(xi02,mark='r-')

      #compute the combined correlation function
      sbins, ns, mubins, nmu, DD,DR,RR=combine_count(root1,root2)
      root=rootN+'COMB'+str(reg)+'-N'+str(ii)
      srebins, murebins, xi2drebin=compute_xi02(root,sbins, ns, mubins, nmu,
                                  DD,DR,RR, 1.0,nscomb=nscomb,write=1)
      xi02=Xi_Legendre(root,srebins, murebins, xi2drebin,outfold=outfoldcomb,write=1)
      pl.figure(3)
      plot_xi02(xi02,mark='bo-',lab='COMB')
      outplot='plots/Reg-N'+str(ii)+'.png'
      f_plot.append(outplot)
      pl.savefig(outplot)
      pl.close(3)
      pl.figure(2)
      plot_xi02(xi02,mark='r-')

      reg=reg+1

   #compute south only and north + south correlations
   root2=rootN+'All-N'
   for ii in range(0,njack_S):
      root1=rootS+'JN'+str(ii)+'-S'
      #compute the correlation functions
      sbins, ns, mubins, nmu, DD,DR,RR, Wsum_data,Wsum_rand=load_PairCount_output(root1)
      Wrat=Wsum_data/Wsum_rand
      srebins, murebins, xi2drebin=compute_xi02(root1,sbins, ns, mubins, nmu,
                                  DD,DR,RR, Wrat,nscomb=nscomb,write=1)
      xi02=Xi_Legendre(root1,srebins, murebins, xi2drebin,outfold=outfold1,write=1)
      pl.figure(3)
      plot_xi02(xi02,mark='ro-',lab='S-'+str(ii))
      pl.figure(1)
      plot_xi02(xi02,mark='b-')

      #compute the combined correlation function
      sbins, ns, mubins, nmu, DD,DR,RR=combine_count(root1,root2)
      root=rootN+'COMB'+str(reg)+'-S'+str(ii)
      srebins, murebins, xi2drebin=compute_xi02(root,sbins, ns, mubins, nmu,
                                  DD,DR,RR, 1.0,nscomb=nscomb,write=1)
      xi02=Xi_Legendre(root,srebins, murebins, xi2drebin,outfold=outfoldcomb,write=1)
      pl.figure(3)
      plot_xi02(xi02,mark='bo-',lab='COMB')
      outplot='plots/Reg-S'+str(ii)+'.png'
      f_plot.append(outplot)
      pl.savefig(outplot)
      pl.close(3)
      pl.figure(2)
      plot_xi02(xi02,mark='b-')

      reg=reg+1

   root1=rootN+'All-N'
   sbins, ns, mubins, nmu, DD,DR,RR, Wsum_data,Wsum_rand=load_PairCount_output(root1)
   Wrat=Wsum_data/Wsum_rand
   srebins, murebins, xi2drebin=compute_xi02(root1,sbins, ns, mubins, nmu,
                                  DD,DR,RR, Wrat,nscomb=nscomb,write=1)
   xi02=Xi_Legendre(root1,srebins, murebins, xi2drebin,outfold=outfold1,write=1)
   pl.figure(3)
   plot_xi02(xi02,mark='ro-',lab=root1)
   pl.figure(1)
   plot_xi02(xi02,mark='r--',lw=1,lab=root1)

   root2=rootS+'All-S'
   sbins, ns, mubins, nmu, DD,DR,RR, Wsum_data,Wsum_rand=load_PairCount_output(root1)
   Wrat=Wsum_data/Wsum_rand
   srebins, murebins, xi2drebin=compute_xi02(root1,sbins, ns, mubins, nmu,
                                  DD,DR,RR, Wrat,nscomb=nscomb,write=1)
   xi02=Xi_Legendre(root1,srebins, murebins, xi2drebin,outfold=outfold1,write=1)
   pl.figure(3)
   plot_xi02(xi02,mark='bo-',lab=root2)
   pl.figure(1)
   plot_xi02(xi02,mark='b--',lw=1,lab=root2)

   sbins, ns, mubins, nmu, DD,DR,RR=combine_count(root1,root2)
   srebins, murebins, xi2drebin=compute_xi02('All-NS',sbins, ns, mubins, nmu,
                               DD,DR,RR, 1.0,nscomb=nscomb,write=1)
   xi02=Xi_Legendre('All-NS',srebins, murebins, xi2drebin,outfold=outfoldcomb,write=1)
   plot_xi02(xi02,mark='k--',lab=root2)
   outplot='plots/All-NS.png'
   f_plot.append(outplot)
   pl.savefig(outplot)
   pl.close(3)

   pl.figure(1)
   pl.title('Individual region, red=N, blue=S',fontsize=32)
   outplot='plots/JN-All-xi02.png'
   f_plot.append(outplot)
   pl.savefig(outplot)

   pl.figure(2)
   plot_xi02(xi02,mark='k--',lab='All-NS')
   pl.title('COMB xi JN region, red=N, blue=S',fontsize=32)
   outplot='plots/JN-COMB-All-xi02.png'
   f_plot.append(outplot)
   pl.savefig(outplot)

def plot_CMASS():
   file='/global/u1/s/shadaba/Projects/Cosmomc/dr11_Comb_8mpc.txt'
   covfile='/global/u1/s/shadaba/Projects/Cosmomc/sdss_DR11_cov_xi02_8Mpch-0-24.mat'
   data=np.loadtxt(file)
   cov=np.loadtxt(covfile)
   diag_err=np.sqrt(np.diagonal(cov))

   rr=data[:-1,0]
   rr2=rr*rr

   pl.errorbar(rr,rr2*data[:-1,1],yerr=rr2*diag_err[:24],fmt='o',color='k',markersize=5,label='CMASS')
   pl.errorbar(rr,rr2*data[:-1,2],yerr=rr2*diag_err[24:],fmt='o',color='k',markersize=5)

   #pl.plot(rr,rr2*data[:,1],'ko-',linewidth=2,markersize=4,label='CMASS')
   #pl.plot(rr,rr2*data[:,2],'ko-',linewidth=2,markersize=4)
   pl.legend()

if __name__ == "__main__":
  #To compute xi02 and xi2d in rmu sampling
 if(1):
   print sys.argv
   if(len(sys.argv)!=4):
      xidir='example/' ; froot='out'
      root=xidir+'PairCount/'+froot
      xi2droot=xidir+'XI2D/'+froot
      xi02root=xidir+'XI02/'+froot
   else:
      root=sys.argv[1]
      xi2droot=sys.argv[2]
      xi02root=sys.argv[3]

   xi2d_wp_xi02_froot_jn(root,xi2droot,xi02root,nscomb=1,samp='rmu',NJN=0)

 print 'List of files written:',
 for ff in f_write:
     print ff

 print 'List of plots created:',
 for pp in f_plot:
     print pp

max_stars_repo_name: echaussidon/LSS

max_issues_repo_name: echaussidon/LSS

max_forks_repo_name: echaussidon/LSS

reason: import numpy,import scipy,from scipy


text:
"""
Tests specific to the patches module.
"""

from numpy.testing import assert_array_equal
from matplotlib.patches import Polygon

def test_Polygon_close():
    """
    Github issue #1018 identified a bug in the Polygon handling
    of the closed attribute; the path was not getting closed
    when set_xy was used to set the vertices.
    """
    # open set of vertices:
    xy = [[0,0], [0,1], [1,1]]
    # closed set:
    xyclosed = xy + [[0,0]]

    # start with open path and close it:
    p = Polygon(xy, closed=True)
    assert_array_equal(p.get_xy(), xyclosed)
    p.set_xy(xy)
    assert_array_equal(p.get_xy(), xyclosed)

    # start with closed path and open it:
    p = Polygon(xyclosed, closed=False)
    assert_array_equal(p.get_xy(), xy)
    p.set_xy(xyclosed)
    assert_array_equal(p.get_xy(), xy)

    # start with open path and leave it open:
    p = Polygon(xy, closed=False)
    assert_array_equal(p.get_xy(), xy)
    p.set_xy(xy)
    assert_array_equal(p.get_xy(), xy)

    # start with closed path and leave it closed:
    p = Polygon(xyclosed, closed=True)
    assert_array_equal(p.get_xy(), xyclosed)
    p.set_xy(xyclosed)
    assert_array_equal(p.get_xy(), xyclosed)

max_stars_repo_name: andreas-h/matplotlib

max_issues_repo_name: andreas-h/matplotlib

max_forks_repo_name: andreas-h/matplotlib

reason: from numpy


text:
import numpy as np
import os
from keras.models import load_model
from keras.applications import imagenet_utils
from tqdm import tqdm
import tifffile
import cv2
import glob
from densemapnet import densemapnet
import argparse

NO_DATA = -999.0
COMPLETENESS_THRESHOLD_PIXELS = 3.0

#GPU="0" # default GPU
GPU="-1" # no GPU

if __name__=="__main__":

    # parse arguments
    parser = argparse.ArgumentParser()
    parser.add_argument('test_folder', type=str)
    parser.add_argument('output_folder', type=str)
    parser.add_argument('model_file', type=str)
    args = parser.parse_args()

    # load the model
    height = 1024
    width = 1024
    bands = 3
    print(height, width, bands)
    os.environ["CUDA_VISIBLE_DEVICES"]=GPU
    settings = densemapnet.Settings()
    settings.xdim = width
    settings.ydim = height
    settings.channels = bands
    settings.nopadding = True
    settings.dropout_test = False
    settings.model_weights = args.model_file
    network = densemapnet.DenseMapNet(settings=settings)
    model = network.build_model()

    # predict disparities for all images in folder
    files = glob.glob(args.test_folder + '*LEFT_RGB.tif')
    nfiles = len(files)
    print('Number of files = ', nfiles)
    for i in tqdm(range(nfiles)):

        name = files[i]
        pos = name.find('LEFT_RGB')
        left_name = name
        right_name = name[0:pos] + 'RIGHT_RGB.tif'
        name = os.path.basename(name)
        pos = name.find('LEFT_RGB')
        dsp_name = args.output_folder + name[0:pos] + 'LEFT_DSP.tif'
        viz_name = args.output_folder + name[0:pos] + 'STEREO_GRAY.tif'
        left = tifffile.imread(left_name)
        right = tifffile.imread(right_name)
        left = np.expand_dims(left,axis=0)
        right = np.expand_dims(right,axis=0)

        # scale image values to [0,1]
        left = (left - 127.5)/255.0
        right = (right - 127.5)/255.0

        disp = model.predict([left, right])[0,:,:,:]
        tifffile.imsave(dsp_name, disp, compress=6)

        # save grayscale version of image for visual inspection
        disp = disp - disp.min()
        disp = ((disp / disp.max()) * 255.0).astype(np.uint8)
        disp = cv2.cvtColor(disp,cv2.COLOR_GRAY2RGB)
        tifffile.imsave(viz_name, disp, compress=6)