text:
import numpy as np
import pandas as pd
import time, os
from sklearn.model_selection import GridSearchCV
import matplotlib.pyplot as plt

class test_model:
    def __init__(self):
        # set test and train score
        self.test_acc = []
        self.train_acc = []

    def score_data(self, model, datadict):
        # get input
        self.model = model
        self.data = datadict
        # fit model
        self.model.fit(self.data['X_train'], self.data['y_train'])
        # train score
        train_score = self.model.score(self.data['X_train'], self.data['y_train'])
        print(f'\tTrain Score: \t\t{train_score}')
        self.train_acc.append(train_score)
        # test score
        test_score = self.model.score(self.data['X_test'], self.data['y_test'])
        print(f'\tTest Score: \t\t{test_score}')
        self.test_acc.append(test_score)
        # print comparison
        print(f'\tOverfit (train - test): \t{train_score - test_score}')

def test_train_score_capture(model, data, train_acc, test_acc):
    # fit model
    model.fit(data['X_train'], data['y_train'])
    # train score
    train_score = model.score(data['X_train'], data['y_train'])
    print(f"\tTrain Score: {train_score}")
    train_acc.append(train_score)
    # test score
    test_score = model.score(data['X_test'], data['y_test'])
    print(f"\tTest Score: {test_score}")
    test_acc.append(test_score)
    return train_acc, test_acc

class ml_timer:
    def __init__(self):
        # initate session start time
        self.start = time.time()
        # initiate holder for times
        self.session_times = []
        print('timer starts')

    def param_start(self):
        # update current session start time
        self.current_session_start = time.time()

    def param_end(self):
        end_time = time.time()
        elapsed_time = (end_time - self.current_session_start)/60
        print(f'\telapsed time {elapsed_time:.3f} min')
        self.session_times.append(elapsed_time)

    def session_end(self):
        self.end = time.time()

    def get_time(self):
        self.runtime = (self.end - self.start)/60
        print(f'total time: {self.runtime:.3f} min')
        return self.runtime

class GridSearchCVSOP:
    def __init__(self, model, grid):
        self.grid = grid
        self.model = model

    def run(self, X, y, **kwargs):
        # process kwargs
        n_jobs = kwargs.pop('n_jobs', -1)
        cv = kwargs.pop('cv', 3)
        scoring = kwargs.pop('scoring', 'accuracy')
        error_score = kwargs.pop('error_score', 0)
        verbose = kwargs.pop('verbose', 3)
        # instantiate grid search
        self.grid_search = GridSearchCV(estimator=self.model,
                                param_grid=self.grid,
                                n_jobs=n_jobs, verbose=verbose,
                                cv=cv, scoring=scoring,
                                error_score=error_score)
        self.grid_result = self.grid_search.fit(X, y)
        return self.grid_result, self.grid_search

    def print_result(self):
        print("Best: %f using %s" % (self.grid_result.best_score_,
                                    self.grid_result.best_params_))
        means = self.grid_result.cv_results_['mean_test_score']
        stds = self.grid_result.cv_results_['std_test_score']
        params = self.grid_result.cv_results_['params']
        for mean, stdev, param in zip(means, stds, params):
            print("%f (%f) with: %r" % (mean, stdev, param))
        return means, stds, params

def load_nutcracker_csv(dir_datafolder):
    datatype = ['X_train','X_test','y_train','y_test']
    print(f'loading {len(datatype)} files')
    datadict = dict()
    for i, dname in enumerate(datatype):
        print(f'loading file: {i}', end='\r')
        filename = 'nutcracker' + '_' + dname + '.csv'
        filepath = os.path.join(dir_datafolder, filename)
        data = pd.read_csv(filepath, header=None, index_col=False)
        datadict[dname] = data.to_numpy()
    print('\nloading completed')
    return datadict

class ModelEvaluation:
    def __init__(self, model, data_dir):
        self.model = model
        self.data_dir = data_dir

    def load_data(self):
        if not hasattr(self, 'data'):
            self.data = load_nutcracker_csv(self.data_dir)

    def cross_val_score(self, cv=5):
        timer = ml_timer()
        if not hasattr(self, 'data'):
            self.load_data()
        from sklearn.model_selection import cross_val_score
        scores = cross_val_score(self.model,
                    self.data['X_train'],
                    self.data['y_train'],
                    cv=cv)
        timer.session_end()
        self.runtime_crossval = timer.get_time()
        print(f'cross validation scores: {scores}')
        print(f'validation score (mean):{np.mean(scores)}')
        print(f'validation score (std):{np.std(scores)}')
        self.cross_val_score_ = scores
        return self.cross_val_score_

    def fitmodel(self):
        if not hasattr(self, 'data'):
            self.load_data()
        self.model.fit(self.data['X_train'], self.data['y_train'])
        return self.model

    def predict(self):
        if not hasattr(self, 'data'):
            self.load_data()
        timer = ml_timer()
        self.y_pred_test = self.model.predict(self.data['X_test'])
        timer.session_end()
        self.runtime_predict = timer.get_time()
        self.y_pred_train = self.model.predict(self.data['X_train'])
        return self.y_pred_test

    def accuracy_score(self):
        if not hasattr(self, 'data'):
            self.load_data()
        self.score_train = self.model.score(self.data["X_train"], self.data['y_train'])
        print(f'accuracy score on train: {self.score_train}')
        self.score_test = self.model.score(self.data['X_test'], self.data['y_test'])
        print(f'accuracy score on test: {self.score_test}')
        return self.score_train, self.score_test

    def confusion_matrix(self):
        if not hasattr(self, 'y_pred_test'):
            self.predict()
        # fitmodel and predict must proceed this.
        # define dataframe labels
        columns = ['Predicted normal', 'Predicted alcohol']
        indexname = ['True normal', 'True alcohol']
        from sklearn.metrics import confusion_matrix
        # run confusion matrix - test
        self.conf_matrix_test = confusion_matrix(self.data['y_test'], self.y_pred_test,
                                            normalize='true')
        conf_matrix_test_df = pd.DataFrame(self.conf_matrix_test, columns=columns)
        conf_matrix_test_df.index = indexname
        print('\nconfusion matrix: test data')
        print(conf_matrix_test_df)
        # run confusion matrix - train
        self.conf_matrix_train = confusion_matrix(self.data['y_train'],
                                            self.y_pred_train,
                                            normalize='true')
        conf_matrix_train_df = pd.DataFrame(self.conf_matrix_train, columns=columns)
        conf_matrix_train_df.index = indexname
        print('\nconfusion matrix: train data')
        print(conf_matrix_train_df)
        return self.conf_matrix_test, self.conf_matrix_train

    def display_confusion_matrix(self):
        if not hasattr(self, 'conf_matrix_test'):
            self.confusion_matrix()
        # confusion_matrix must proceed this
        display_labels = ['normal', 'alcohol']
        from sklearn.metrics import ConfusionMatrixDisplay
        print('\nconfusion matrix for test')
        plt.figure()
        ConfusionMatrixDisplay(self.conf_matrix_test, display_labels=display_labels).plot()
        plt.show()
        print('\nconfusion matrix for train')
        plt.figure()
        ConfusionMatrixDisplay(self.conf_matrix_train, display_labels=display_labels).plot()
        plt.show()

    def classification_report(self):
        if not hasattr(self, 'y_pred_test'):
            self.predict()
        from sklearn.metrics import classification_report
        self.eval_score_report = classification_report(self.data['y_test'],
                                                        self.y_pred_test)
        print(self.eval_score_report)

    def print_evaluation_scores(self):
        if not hasattr(self, 'y_pred_test'):
            self.predict()
        from sklearn.metrics import precision_score
        from sklearn.metrics import recall_score
        from sklearn.metrics import f1_score
        self.precision_score = precision_score(self.data["y_test"], self.y_pred_test)
        self.recall_score = recall_score(self.data["y_test"], self.y_pred_test)
        self.f1_score = f1_score(self.data["y_test"], self.y_pred_test)
        print(f'precision_score = {self.precision_score}')
        print(f'recall_score = {self.recall_score}')
        print(f'f1_score = {self.f1_score}')

    def test_data_class_proba(self):
        if not hasattr(self, 'data'):
            self.load_data()
        false_proba = np.count_nonzero(self.data['y_test']) / self.data['y_test'].shape[0]
        true_proba = 1.0 - false_proba
        print(f'test set normal case probability: {false_proba}')
        print(f'test set alcohol case probability: {true_proba}')
        self.real_proba = dict()
        self.real_proba['false_proba'] = false_proba
        self.real_proba['true_proba'] = true_proba

    def predict_proba(self):
        if not hasattr(self, 'data'):
            self.load_data()
        self.y_proba_test = self.model.predict_proba(self.data['X_test'])[:,1]
        self.y_proba_train = self.model.predict_proba(self.data['X_train'])[:,1]
        return self.y_proba_test, self.y_proba_train

    def proba_thresholds(self):
        if not hasattr(self, 'data'):
            self.load_data()
        from sklearn.metrics import accuracy_score
        from sklearn.metrics import precision_score
        from sklearn.metrics import recall_score
        # Vary thresholds by 0.05 from 0.05 to 0.95
        thresholds = np.arange(0.05, 1, 0.05)
        precisions = list()
        recalls = list()
        neg_recalls = list()
        for threshold in thresholds:
            # Apply threshold
            y_threshold = np.where(self.y_proba_test > threshold, 1, 0)
            # Calculate precision and recall
            precision = precision_score(self.data['y_test'], y_threshold)
            recall = recall_score(self.data['y_test'], y_threshold)
            neg_recall = recall_score(1-self.data['y_test'], 1-y_threshold)
            # Append to list
            precisions.append(precision)
            recalls.append(recall)
            neg_recalls.append(neg_recall)
        # Visualize the result
        plt.figure()
        plt.plot(thresholds, precisions, label='precision', marker='o')
        plt.plot(thresholds, recalls, label='recall', marker='o')
        plt.xlim(0, 1)
        plt.xlabel('threshold')
        plt.ylabel('score')
        plt.legend()
        plt.show()
        return precisions, recalls, neg_recalls

    def roc_auc(self):
        if hasattr(self, 'y_proba_train'):
            self.predict_proba()
        from sklearn.metrics import roc_curve, roc_auc_score
        # get roc auc train
        fprs_train, tprs_train, thresholds_train = roc_curve(self.data['y_train'], self.y_proba_train)
        roc_auc_train = roc_auc_score(self.data['y_train'], self.y_proba_train)
        # get roc auc test
        fprs_test, tprs_test, thresholds_test = roc_curve(self.data['y_test'], self.y_proba_test)
        roc_auc_test = roc_auc_score(self.data['y_test'], self.y_proba_test)
        # Plot the ROC curve.
        plt.figure()
        plt.plot(fprs_train, tprs_train, color='gray', lw=5, label='train', linestyle=' ', marker='.')
        plt.plot(fprs_test, tprs_test, lw=1, color='red', label='test')
        plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--', label='expected')
        plt.xlabel('False Positive Rate')
        plt.ylabel('True Positive Rate')
        plt.title('ROC and AUC')
        plt.legend(loc="best")
        plt.show()
        print(f"Test AUC score: {roc_auc_test}")
        print(f"Train AUC score: {roc_auc_train}")
        self.roc_auc_test = roc_auc_test
        self.roc_auc_train = roc_auc_train

    def update_data_dir(self, data_dir):
        if not hasattr(self, data_dir):
            self.data_dir = data_dir
        else:
            if self.data_dir != data_dir:
                self.data_dir_original = self.data_dir
                self.data_dir = data_dir

    def save(self, savedir):
        # remove data from object to save space
        if hasattr(self, 'data'):
            delattr(self, 'data')
        if hasattr(self, 'y_pred_test'):
            delattr(self, 'y_pred_test')
        if hasattr(self, 'y_pred_train'):
            delattr(self, 'y_pred_train')
        if hasattr(self, 'y_proba_test'):
            delattr(self, 'y_proba_test')
        if hasattr(self, 'y_proba_train'):
            delattr(self, 'y_proba_train')
        # get model name
        model_type = type(self.model)
        model_type_str  = str(model_type)
        model_name_components = model_type_str.split('.')
        model_name = model_name_components[-1].replace("'>","")
        # save
        import pickle, os
        savepath = os.path.join(savedir, model_name+'_eval.pickle')
        pickle.dump(self, open(savepath, 'wb'))

    def excel_input_array(self):
        report = [np.mean(self.cross_val_score_),
                    np.std(self.cross_val_score_),
                    self.score_train,
                    self.score_test,
                    self.precision_score,
                    self.recall_score,
                    self.f1_score]
        if hasattr(self, 'roc_auc_train'):
            report.append(self.roc_auc_train)
            report.append(self.roc_auc_test)
        report.append(self.runtime_crossval)
        report.append(self.runtime_predict)
        print(report)
        print(self.model)

    def standard(self, save_dir):
        if not hasattr(self.model, 'predict_proba'):
            print('confirm: model has no predict_proba')
        print('\nloading data from directory')
        self.load_data()
        print('\nfit model...')
        self.fitmodel()
        print('predict model...')
        self.predict()
        print('\naccuracy scores:')
        self.accuracy_score()
        print('\nconfusion matrix:')
        self.confusion_matrix()
        self.display_confusion_matrix()
        print('\nclassification report:')
        self.classification_report()
        self.print_evaluation_scores()
        if hasattr(self.model, 'predict_proba'):
            print('\nreal data class proba:')
            self.test_data_class_proba()
            print('\n prediction proba:')
            self.predict_proba()
            print('\nproba threshold analysis:')
            self.proba_thresholds()
            print('\nROC AUC analysis:')
            self.roc_auc()
            print('\nSaving model...')
        else:
            print('\nthis model does not have predict_proba attr')
        print('\nruning cross validation scores (this takes a while):')
        self.cross_val_score(5)
        print('\nexcel record:')
        self.excel_input_array()
        print('\nsave model')
        self.save(save_dir)


text:
import numpy as np
import matplotlib.pyplot as plt
import math

def squish(sample, squishFactor):
    #check for input errors
    if squishFactor < 0:
        raise EOFError

    #fix a clipping signal
    if sample > 1:
        return 1
    elif sample < -1:
        return -1

    #compress signal between -1 and 1
    exp = 1/squishFactor
    modified = 0
    if sample > 0:
        modified = -(1 - sample)**exp + 1
    else:
        modified = -(-(1 + sample)**exp + 1)
    return modified

#constants

xMax = math.pi
yMax = 2.0
dx = 0.01

#X values
xArray = np.arange(-xMax, xMax, dx)
signal = []
for x in xArray:
    signal.append(math.sin(x))

plt.plot(xArray, signal)

# squishFactors = [0.01, 0.1, .5, 2, 10, 100]
squishFactors = [1, 10, 100]

for sf in squishFactors:
    modifiedSignal = []
    for sample in signal:
        modifiedSignal.append(squish(sample, sf))
    plt.plot(xArray, modifiedSignal)

# for sf in squishFactors:
#     modifiedSignal = []
#     for sample in signal:
#         if sample < .5 :
#             sample = sample + sample*sf
#         else :
#             sample = sample + (1-sample)*sf

#         modifiedSignal.append(sample)
#     plt.plot(xArray, modifiedSignal)

# for sf in squishFactors:
#     modifiedSignal = []
#     for sample in signal:
#         compression = sample**1.25
#         newSample = sample*sf - compression*(sf-1)

#         modifiedSignal.append(newSample)
#     plt.plot(xArray, modifiedSignal)

plt.xlabel("x")
plt.ylabel("y")
plt.title("Spenser's Title")
plt.legend()
plt.axis([-xMax, xMax, -yMax, yMax])
plt.grid(color='g', linestyle='--', linewidth=1, axis='both', markevery=1)
plt.show()

ext: py

max_stars_repo_name: spensbot/crispy

max_issues_repo_name: spensbot/crispy

max_forks_repo_name: spensbot/crispy

reason: import numpy


text:
# Note: The codes were originally created by Prof. Jack Baker in the MATLAB

import numpy as np
from scipy.stats import norm
from scipy.interpolate import interp1d

# CY_14 attenuation equation - Update of the Chiou and Youngs NGA Model for
# the Average Horizontal Component of Peak Ground Motion and Response
# Spectra, Chiou and Youngs (2014)

##############
### Inputs ###
##############

# M     = Moment Magnitude
# T     = Period (sec); Use Period = -1 for PGV computation
#                 Use 1000 for output the array of Sa with original period
#                 (no interpolation)
# Rup   = Closest distance (km) to the ruptured plane
# Rjb   = Joyner-Boore distance (km); closest distance (km) to surface
#       projection of rupture plane
# Rx    =    Horizontal distance from top of rupture measured perpendicular
#       to fault strike (km). See Figures a, b and c for illustation% Ztor  = Depth(km) to the top of ruptured plane
# delta = Fault dip angle (in degree)
# lamda = Rake angle      (in degree)
# region        = 0 for global (incl. Taiwan)
#               = 1 for California
#               = 2 for Japan
#               = 3 for China
#               = 4 for Italy
#               = 5 for Turkey
# Z10            = Basin depth (km); depth from the groundsurface to the
#                   1km/s shear-wave horizon.
#               = 999 if unknown
# Vs30          = shear wave velocity averaged over top 30 m in m/s
#               = ref: 1130
# FVS30         = 1 for Vs30 is inferred from geology
#               = 0 for measured  Vs30
# d_DPP         = DPP centered on the site and earthquake specific average
#               DPP, = 0 for median calc
#               %d_DPP functionality added BAB 26 July 2020

###############
### Outputs ###
###############

# Sa: Median spectral acceleration prediction
# sigma: logarithmic standard deviation of spectral acceleration
#          prediction

def CY_2014_sub(M, ip, R_RUP, R_JB, Rx, Ztor, delta, F_RV, F_NM, HW, Z10, Vs30, FVS30, region, d_DPP):

    dict = {}
    dict["c2"] = [1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06, 1.06]
    dict["c4"] = [-2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1, -2.1]
    dict["c4_a"] = [-0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5]
    dict["c_RB"] = [50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50]
    dict["c8"] = [0.2154, 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.0991, 0.1982, 0.2154, 0.2154, 0.2154, 0.2154, 0.2154, 0.2154, 0.2154, 0.2154]
    dict["c8_a"] = [0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695, 0.2695]
    dict["c1"] = [ 2.3549, -1.5065, -1.5065, -1.4798, -1.2972, -1.1007, -0.9292, -0.658 , -0.5613, -0.5342, -0.5462, -0.5858, -0.6798, -0.8663, -1.0514, -1.3794, -1.6508, -2.1511, -2.5365, -3.0686, -3.4148, -3.9013, -4.2466, -4.5143, -5.0009, -5.3461]
    dict["c1_a"] = [0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165, 0.165, 0.165, 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.165 , 0.1645, 0.1168, 0.0732, 0.0484, 0.022 , 0.0124]
    dict["c1_b"] = [-0.0626, -0.255 , -0.255 , -0.255 , -0.255 , -0.255 , -0.255 , -0.254 , -0.253 , -0.252 , -0.25  , -0.248 , -0.2449, -0.2382, -0.2313, -0.2146, -0.1972, -0.162 , -0.14  , -0.1184, -0.11  , -0.104 , -0.102 , -0.101 , -0.101 , -0.1   ]
    dict["c1_c"] = [-0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.165 , -0.1645, -0.1168, -0.0732, -0.0484, -0.022 , -0.0124]
    dict["c1_d"] = [0.0626, 0.255 , 0.255 , 0.255 , 0.255 , 0.255 , 0.255 , 0.254 , 0.253 , 0.252 , 0.25  , 0.248 , 0.2449, 0.2382, 0.2313, 0.2146, 0.1972, 0.162 , 0.14  , 0.1184, 0.11  , 0.104 , 0.102 , 0.101 , 0.101 , 0.1   ]
    dict["c_n"] = [ 3.3024, 16.0875, 16.0875, 15.7118, 15.8819, 16.4556, 17.6453, 20.1772, 19.9992, 18.7106, 16.6246, 15.3709, 13.7012, 11.2667, 9.1908,  6.5459,  5.2305,  3.7896,  3.3024,  2.8498,  2.5417, 2.1488,  1.8957,  1.7228,  1.5737,  1.5265]
    dict["c_m"] = [5.423 , 4.9993, 4.9993, 4.9993, 4.9993, 4.9993, 4.9993, 5.0031, 5.0172, 5.0315, 5.0547, 5.0704, 5.0939, 5.1315, 5.167 , 5.2317, 5.2893, 5.4109, 5.5106, 5.6705, 5.7981, 5.9983, 6.1552, 6.2856, 6.5428, 6.7415]
    dict["c3"] = [2.3152, 1.9636, 1.9636, 1.9636, 1.9636, 1.9636, 1.9636, 1.9636, 1.9636, 1.9795, 2.0362, 2.0823, 2.1521, 2.2574, 2.344 , 2.4709, 2.5567, 2.6812, 2.7474, 2.8161, 2.8514, 2.8875, 2.9058, 2.9169, 2.932 , 2.9396]
    dict["c5"] = [5.8096, 6.4551, 6.4551, 6.4551, 6.4551, 6.4551, 6.4551, 6.4551, 6.8305, 7.1333, 7.3621, 7.4365, 7.4972, 7.5416, 7.56  , 7.5735, 7.5778, 7.5808, 7.5814, 7.5817, 7.5818, 7.5818, 7.5818, 7.5818, 7.5818, 7.5818]
    dict["c_HM"] = [3.0514, 3.0956, 3.0956, 3.0963, 3.0974, 3.0988, 3.1011, 3.1094, 3.2381, 3.3407, 3.43  , 3.4688, 3.5146, 3.5746, 3.6232, 3.6945, 3.7401, 3.7941, 3.8144, 3.8284, 3.833 , 3.8361, 3.8369, 3.8376, 3.838 , 3.838 ]
    dict["c6"] = [0.4407, 0.4908, 0.4908, 0.4925, 0.4992, 0.5037, 0.5048, 0.5048, 0.5048, 0.5048, 0.5045, 0.5036, 0.5016, 0.4971, 0.4919, 0.4807, 0.4707, 0.4575, 0.4522, 0.4501, 0.45  , 0.45  , 0.45  , 0.45  , 0.45  , 0.45  ]
    dict["c7"] = [0.0324, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.0352, 0.016 , 0.0062, 0.0029, 0.0007, 0.0003]
    dict["c7_b"] = [ 0.0097,  0.0462,  0.0462,  0.0472,  0.0533,  0.0596,  0.0639, 0.063 ,  0.0532,  0.0452,  0.0345,  0.0283, 0.0202,  0.009 , -0.0004, -0.0155, -0.0278, -0.0477, -0.0559, -0.063 , -0.0665, -0.0516, -0.0448, -0.0424, -0.0348, -0.0253]
    dict["c8_b"] = [5.    , 0.4833, 0.4833, 1.2144, 1.6421, 1.9456, 2.181 , 2.6087, 2.9122, 3.1045, 3.3399, 3.4719, 3.6434, 3.8787, 4.0711, 4.3745, 4.6099, 5.0376, 5.3411, 5.7688, 6.0723, 6.5   , 6.8035, 7.0389, 7.4666, 7.77  ]
    dict["c9"] = [0.3079, 0.9228, 0.9228, 0.9296, 0.9396, 0.9661, 0.9794, 1.026 , 1.0177, 1.0008, 0.9801, 0.9652, 0.9459, 0.9196, 0.8829, 0.8302, 0.7884, 0.6754, 0.6196, 0.5101, 0.3917, 0.1244, 0.0086, 0.    , 0.    , 0.    ]
    dict["c9_a"] = [0.1   , 0.1202, 0.1202, 0.1217, 0.1194, 0.1166, 0.1176, 0.1171, 0.1146, 0.1128, 0.1106, 0.115 , 0.1208, 0.1208, 0.1175, 0.106 , 0.1061, 0.1   , 0.1   , 0.1   , 0.1   , 0.1   , 0.1   , 0.1   , 0.1   , 0.1   ]
    dict["c9_b"] = [6.5   , 6.8607, 6.8607, 6.8697, 6.9113, 7.0271, 7.0959, 7.3298, 7.2588, 7.2372, 7.2109, 7.2491, 7.2988, 7.3691, 6.8789, 6.5334, 6.526 , 6.5   , 6.5   , 6.5   , 6.5   , 6.5   , 6.5   , 6.5   , 6.5   , 6.5   ]
    dict["c11"] = [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]
    dict["c11_b"] = [-0.3834, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.4536, -0.444 , -0.3539, -0.2688, -0.1793, -0.1428, -0.1138, -0.1062, -0.102 , -0.1009, -0.1003, -0.1001, -0.1001, -0.1, -0.1]
    dict["c_g1"] = [-0.001852, -0.007146, -0.007146, -0.007249, -0.007869, -0.008316, -0.008743, -0.009537, -0.00983 , -0.009913, -0.009896, -0.009787, -0.009505, -0.008918, -0.008251, -0.007267, -0.006492, -0.005147, -0.004277, -0.002979, -0.002301, -0.001344, -0.001084, -0.00101 , -0.000964, -0.00095 ]
    dict["c_g2"] = [-0.007403, -0.006758, -0.006758, -0.006758, -0.006758, -0.006758, -0.006758, -0.00619 , -0.005332, -0.004732, -0.003806, -0.00328 , -0.00269 , -0.002128, -0.001812, -0.001274, -0.001074, -0.001115, -0.001197, -0.001675, -0.002349, -0.003306, -0.003566, -0.00364 , -0.003686, -0.0037  ]
    dict["c_g3"] = [4.3439, 4.2542, 4.2542, 4.2386, 4.2519, 4.296 , 4.3578, 4.5455, 4.7603, 4.8963, 5.0644, 5.1371, 5.188 , 5.2164, 5.1954, 5.0899, 4.7854, 4.3304, 4.1667, 4.0029, 3.8949, 3.7928, 3.7443, 3.709 , 3.6632, 3.623 ]
    dict["phi1"] = [-0.7936, -0.521 , -0.521 , -0.5055, -0.4368, -0.3752, -0.3469, -0.3747, -0.444 , -0.4895, -0.5477, -0.5922, -0.6693, -0.7766, -0.8501, -0.9431, -1.0044, -1.0602, -1.0941, -1.1142, -1.1154, -1.1081, -1.0603, -0.9872, -0.8274, -0.7053]
    dict["phi2"] = [-0.0699, -0.1417, -0.1417, -0.1364, -0.1403, -0.1591, -0.1862, -0.2538, -0.2943, -0.3077, -0.3113, -0.3062, -0.2927, -0.2662, -0.2405, -0.1975, -0.1633, -0.1028, -0.0699, -0.0425, -0.0302, -0.0129, -0.0016,  0.    ,  0.    ,  0.    ]
    dict["phi3"] = [-0.008444, -0.00701 , -0.00701 , -0.007279, -0.007354, -0.006977, -0.006467, -0.005734, -0.005604, -0.005696, -0.005845, -0.005959, -0.006141, -0.006439, -0.006704, -0.007125, -0.007435, -0.00812 , -0.008444, -0.007707, -0.004792, -0.001828, -0.001523, -0.00144 , -0.001369, -0.001361]
    dict["phi4"] = [5.41000e+00, 1.02151e-01, 1.02151e-01, 1.08360e-01, 1.19888e-01, 1.33641e-01, 1.48927e-01, 1.90596e-01, 2.30662e-01, 2.53169e-01, 2.66468e-01, 2.65060e-01, 2.55253e-01, 2.31541e-01, 2.07277e-01, 1.65464e-01, 1.33828e-01, 8.51530e-02, 5.85950e-02, 3.17870e-02, 1.97160e-02, 9.64300e-03, 5.37900e-03, 3.22300e-03, 1.13400e-03, 5.15000e-04]
    dict["phi5"] = [0.0202, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.001 , 0.004 , 0.01  , 0.034 , 0.067 , 0.143 , 0.203, 0.277, 0.309 , 0.321 , 0.329 , 0.33  ]
    dict["phi6"] = [300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300]
    dict["tau1"] = [0.3894, 0.4   , 0.4   , 0.4026, 0.4063, 0.4095, 0.4124, 0.4179, 0.4219, 0.4244, 0.4275, 0.4292, 0.4313, 0.4341, 0.4363, 0.4396, 0.4419, 0.4459, 0.4484, 0.4515, 0.4534, 0.4558, 0.4574, 0.4584, 0.4601, 0.4612]
    dict["tau2"] = [0.2578, 0.26  , 0.26  , 0.2637, 0.2689, 0.2736, 0.2777, 0.2855, 0.2913, 0.2949, 0.2993, 0.3017, 0.3047, 0.3087, 0.3119, 0.3165, 0.3199, 0.3255, 0.3291, 0.3335, 0.3363, 0.3398, 0.3419, 0.3435, 0.3459, 0.3474]
    dict["sigma1"] = [0.4785, 0.4912, 0.4912, 0.4904, 0.4988, 0.5049, 0.5096, 0.5179, 0.5236, 0.527 , 0.5308, 0.5328, 0.5351, 0.5377, 0.5395, 0.5422, 0.5433, 0.5294, 0.5105, 0.4783, 0.4681, 0.4617, 0.4571, 0.4535, 0.4471, 0.4426]
    dict["sigma2"] = [0.3629, 0.3762, 0.3762, 0.3762, 0.3849, 0.391 , 0.3957, 0.4043, 0.4104, 0.4143, 0.4191, 0.4217, 0.4252, 0.4299, 0.4338, 0.4399, 0.4446, 0.4533, 0.4594, 0.468 , 0.4681, 0.4617, 0.4571, 0.4535, 0.4471, 0.4426]
    dict["sigma3"] = [0.7504, 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.8   , 0.7999, 0.7997, 0.7988, 0.7966, 0.7792, 0.7504, 0.7136, 0.7035, 0.7006, 0.7001, 0.7   , 0.7   , 0.7   ]
    dict["sigma2_JP"] = [0.3918, 0.4528, 0.4528, 0.4551, 0.4571, 0.4642, 0.4716, 0.5022, 0.523 , 0.5278, 0.5304, 0.531 , 0.5312, 0.5309, 0.5307, 0.531 , 0.5313, 0.5309, 0.5302, 0.5276, 0.5167, 0.4917, 0.4682, 0.4517, 0.4167, 0.3755]
    dict["gamma_JP_IT"] = [2.2306, 1.5817, 1.5817, 1.574 , 1.5544, 1.5502, 1.5391, 1.4804, 1.4094, 1.3682, 1.3241, 1.3071, 1.2931, 1.315 , 1.3514, 1.4051, 1.4402, 1.528 , 1.6523, 1.8872, 2.1348, 3.5752, 3.8646, 3.7292, 2.3763, 1.7679]
    dict["gamma_Wn"] = [0.335 , 0.7594, 0.7594, 0.7606, 0.7642, 0.7676, 0.7739, 0.7956, 0.7932, 0.7768, 0.7437, 0.7219, 0.6922, 0.6579, 0.6362, 0.6049, 0.5507, 0.3582, 0.2003, 0.0356, 0.    , 0.    , 0.    , 0.    , 0.    , 0.    ]
    dict["phi1_JP"] = [-0.7966, -0.6846, -0.6846, -0.6681, -0.6314, -0.5855, -0.5457, -0.4685, -0.4985, -0.5603, -0.6451, -0.6981, -0.7653, -0.8469, -0.8999, -0.9618, -0.9945, -1.0225, -1.0002, -0.9245, -0.8626, -0.7882, -0.7195, -0.656 , -0.5202, -0.4068]
    dict["phi5_JP"] = [0.9488, 0.459 , 0.459 , 0.458 , 0.462 , 0.453 , 0.436 , 0.383 , 0.375 , 0.377 , 0.379 , 0.38  , 0.384 , 0.393 , 0.408 , 0.462 , 0.524 , 0.658 , 0.78  , 0.96  , 1.11  , 1.291 , 1.387 , 1.433 , 1.46  , 1.464 ]
    dict["phi6_JP"] = [800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800, 800]

    if region == 2:
        dict["sigma2"] = dict["sigma2_JP"]
        dict["phi1"] = dict["phi1_JP"]
        dict["phi5"] = dict["phi5_JP"]
        phi6 = phi6_JP

#     fmag
    term6 = dict["c2"][ip]*(M-6)
    term7 = (dict["c2"][ip]-dict["c3"][ip])/dict["c_n"][ip]*np.log(1+np.exp(dict["c_n"][ip]*(dict["c_m"][ip]-M)))

#     Distance Scaling and attenuation term
    term8 = dict["c4"][ip] * np.log(R_RUP + dict["c5"][ip] * np.cosh(dict["c6"][ip] * max(M-dict["c_HM"][ip],0)))
    term9 = (dict["c4_a"][ip] - dict["c4"][ip]) * np.log(np.sqrt(R_RUP**2 + dict["c_RB"][ip]**2))
    term10 = (dict["c_g1"][ip] + dict["c_g2"][ip] / (np.cosh(max(M-dict["c_g3"][ip],0)))) * R_RUP

    if (region == 2 or region == 4):
        if (M>6 and M<6.9):
            term10 = dict["gamma_JP_IT"][ip] * term10

    if region == 3:
        term10 = dict["gamma_Wn"][ip] * term10

#     Style of faulting term
    term2=(dict["c1_a"][ip] + dict["c1_c"][ip] / np.cosh(2 * max(M-4.5,0))) * F_RV
    term3=(dict["c1_b"][ip] + dict["c1_d"][ip] / np.cosh(2 * max(M-4.5,0))) * F_NM

#   Ztor term
    if F_RV == 1:
        E_Ztor = (max(2.704 - 1.226*max(M-5.849,0), 0)) ** 2
    else:
        E_Ztor = (max(2.673 - 1.136*max(M-4.970,0), 0)) ** 2

    if Ztor == 999:
        Ztor = E_Ztor

    delta_ZTOR = Ztor - E_Ztor

    term4 = (dict["c7"][ip] + dict["c7_b"][ip] / np.cosh(2 * max(M-4.5,0))) * delta_ZTOR

#     Hanging wall term
    term12 = dict["c9"][ip] * HW * np.cos(delta) * (dict["c9_a"][ip] + (1 - dict["c9_a"][ip]) * np.tanh(Rx/dict["c9_b"][ip])) * (1 - np.sqrt(R_JB**2 + Ztor**2) / (R_RUP + 1))

#     Basin Depth term
#   Z1.0 (m) ~ Vs30 (m/s) relationship

    if region != 2: # in California and non-Japan region
        z_1 = np.exp(-7.15 / 4 * np.log((Vs30**4 + 570.94**4) / (1360**4 + 570.94**4)))
    else:
        z_1 = np.exp(-5.23 / 2 * np.log((Vs30**2 + 412.39**2) / (1360**2 + 412.39**2)))

    if Z10 ==999:
        d_Z1 = 0
    else:
        d_Z1 = Z10 * 1000 - z_1

#     Dip term
    term5 = (dict["c11"][ip] + dict["c11_b"][ip] / np.cosh(2 * max(M-4.5,0))) * (np.cos(delta)**2)

#     Directivity
    term11 = dict["c8"][ip] * max(1 - max(R_RUP-40,0)/30,0) * min(max(M-5.5,0)/0.8,1) * np.exp(-dict["c8_a"][ip] * (M-dict["c8_b"][ip])**2) * d_DPP
    term1 = dict["c1"][ip]
    ln_yrefij = term1 + term2 + term3 + term4 + term5 + term6 + term7 + term8 + term9 + term10 + term11 + term12
    yrefij = np.exp(ln_yrefij)

#     Site response
    term14 = dict["phi1"][ip] * min(np.log(Vs30/1130),0)
    term15 = dict["phi2"][ip] * (np.exp(dict["phi3"][ip]*(min(Vs30,1130)-360)) - np.exp(dict["phi3"][ip] * (1130-360))) * np.log((yrefij + dict["phi4"][ip]) / dict["phi4"][ip])
    term16 = dict["phi5"][ip] * (1 - np.exp(-d_Z1 / dict["phi6"][ip]))

    Sa = yrefij * np.exp(term14 + term15 + term16)

#     Compute standard deviation
    Finferred = (FVS30==0) # equal 1 if Vs30 is inferred from geology.
    Fmeasured = (FVS30==1) # equal 1 if Vs30 is measured.

    NL0 = dict["phi2"][ip] * (np.exp(dict["phi3"][ip] * (min(Vs30,1130)-360)) - np.exp(dict["phi3"][ip] * (1130-360))) * (yrefij / (yrefij + dict["phi4"][ip]))

    sigmaNL0 = (dict["sigma1"][ip] + (dict["sigma2"][ip] - dict["sigma1"][ip]) / 1.5 * (min(max(M,5),6.5) - 5)) * np.sqrt((dict["sigma3"][ip] * Finferred + 0.7 * Fmeasured) + (1 + NL0)**2)

    tau = dict["tau1"][ip] + (dict["tau2"][ip] - dict["tau1"][ip]) / 1.5 * (min(max(M,5),6.5) - 5)
    sigma = np.sqrt((1 + NL0)**2 * tau**2 + sigmaNL0**2)

    return Sa, sigma

def gmpe_CY_2014(M, T, Rup, Rjb, Rx, Ztor, delta, lambda0, Z10, Vs30, Fhw, FVS30, region, *d_DPP):
    period = [-1, 0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.075, 0.1, 0.12, 0.15, 0.17, 0.2, 0.25, 0.3, 0.4, 0.5, 0.75, 1, 1.5, 2, 3, 4, 5, 7.5, 10]

    delta = delta * np.pi / 180.0
    frv = (lambda0 >= 30 and lambda0 <= 150) # frv: 1 for lambda between 30 and 150, 0 otherwise
    fnm = (lambda0 >= -120 and lambda0 <= -60)  # fnm: 1 for lambda between -120 and -60, 0 otherwise

    if Fhw == 1:
        HW = 1
    elif Fhw == 0:
        HW = 0
    else:
        HW = Rx>=0

    arg = list(d_DPP)
    if len(arg) < 1: # for backward compatibility with prior versions that omit this term
        d_DPP = 0 # for median calculatio, d_DPP=0.

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

    if (len(T) == 1 and T == 1000): # Compute Sa and sigma with pre-defined period
        Sa = np.zeros((1,len(period)-2))
        sigma = np.zeros((1,len(period)-2))
        period1 = period[2:]
        for ip in range(2, len(period)): # 3:length(period)
            Sa[ip-2],sigma[ip-2] = CY_2014_sub(M, ip, Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30,region, d_DPP)
            [PGA,sigma_PGA] = CY_2014_sub(M,1,Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30,region, d_DPP)
            if (Sa[ip-2]<PGA and period1[ip-2]<=0.3):
                Sa[ip-2]= PGA
    else:                            # Compute Sa and sigma with user-defined period
        Sa = np.zeros((1, len(T)))
        sigma = np.zeros((1, len(T)))
        period1 = T
        for i in range(len(T)):
            Ti = T[i]
            if len(np.where(abs(np.array(period)-Ti) < 0.0001)[0]) == 0: # The user defined period requires interpolation
                T_low = max(np.array(period)[np.array(period) < Ti])
                T_high = min(np.array(period)[np.array(period) > Ti])
                ip_low  = period.index(T_low)
                ip_high = period.index(T_high)

                Sa_low, sigma_low = CY_2014_sub(M, ip_low, Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30, region, d_DPP)
                Sa_high, sigma_high = CY_2014_sub(M, ip_high, Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30, region, d_DPP)
                PGA, sigma_PGA = CY_2014_sub(M,1,Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30, region, d_DPP)

                x = [np.log(T_low), np.log(T_high)]
                Y_sa = [np.log(Sa_low), np.log(Sa_high)]
                Y_sigma = [sigma_low, sigma_high]
                interpF1 = interp1d(x, Y_sa)
                Sa[i] = np.exp(interpF1(np.log(Ti)))
                interpF2 = interp1d(x, Y_sigma)
                sigma[i] = interpF2(np.log(Ti))
                if (Sa[i] < PGA and period1[i] <= 0.3):
                    Sa[i]= PGA
            else:
                ip_T = np.where(abs(np.array(period) - Ti) < 0.0001)[0][0]
                Sa[i], sigma[i] = CY_2014_sub(M, ip_T, Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30, region, d_DPP)
                PGA, sigma_PGA = CY_2014_sub(M,1,Rup, Rjb, Rx, Ztor, delta, frv, fnm, HW, Z10, Vs30, FVS30, region, d_DPP)
                if (Sa[i]<PGA and period1[i] <= 0.3):
                    Sa[i] = PGA

    return Sa[0][0], sigma[0][0], period1

max_stars_repo_path: Chapter 6-7/gmpe_CY_2014.py

max_issues_repo_path: Chapter 6-7/gmpe_CY_2014.py

max_forks_repo_path: Chapter 6-7/gmpe_CY_2014.py

reason: import numpy,from scipy


text:
import sympy

ans = sympy.prime(10001)

print(ans)

# Copyright Junipyr. All rights reserved.
# https://github.com/Junipyr

max_stars_repo_path: Problem 007.py

max_stars_repo_name: Junipyr/Project-Euler

max_issues_repo_path: Problem 007.py

max_issues_repo_name: Junipyr/Project-Euler

max_forks_repo_path: Problem 007.py

max_forks_repo_name: Junipyr/Project-Euler

reason: import sympy


text:
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
from constants import _OPTIMIZER, _XTOL, _FTOL, _GTOL, _STEP
from affine import Rigid

from nipy.neurospin.image import apply_affine
from nipy.neurospin.image.image_module import cspline_transform, cspline_sample4d
from nipy.neurospin.utils.optimize import fmin_steepest

import numpy as np
from scipy.optimize import fmin as fmin_simplex, fmin_powell, fmin_cg, fmin_bfgs

_SLICE_ORDER = 'ascending'
_INTERLEAVED = False
_SLICE_AXIS = 2
_SPEEDUP = 4
_WITHIN_LOOPS = 2
_BETWEEN_LOOPS = 5

def interp_slice_order(Z, slice_order):
    Z = np.asarray(Z)
    nslices = len(slice_order)
    aux = np.asarray(list(slice_order)+[slice_order[0]+nslices])
    Zf = np.floor(Z).astype('int')
    w = Z - Zf
    Zal = Zf % nslices
    Za = Zal + w
    """
    Za = Z % nslices
    Zal = Za.astype('int')
    w = Za - Zal
    """
    ret = (1-w)*aux[Zal] + w*aux[Zal+1]
    ret += (Z-Za)
    return ret

def grid_coords(xyz, affine, from_world, to_world):
    Tv = np.dot(from_world, np.dot(affine, to_world))
    XYZ = apply_affine(Tv, xyz)
    return XYZ[:,0], XYZ[:,1], XYZ[:,2]

class Image4d(object):
    """
    Class to represent a sequence of 3d scans acquired on a
    slice-by-slice basis.
    """
    def __init__(self, array, to_world, tr, tr_slices=None, start=0.0,
                 slice_order=_SLICE_ORDER, interleaved=_INTERLEAVED,
                 slice_axis=_SLICE_AXIS):
        """
        Configure fMRI acquisition time parameters.

        tr  : inter-scan repetition time, i.e. the time elapsed
              between two consecutive scans
        tr_slices : inter-slice repetition time, same as tr for slices
        start   : starting acquisition time respective to the implicit
                  time origin
        slice_order : string or array
        """
        self.array = array
        self.to_world = to_world
        nslices = array.shape[slice_axis]

        # Default slice repetition time (no silence)
        if tr_slices == None:
            tr_slices = tr/float(nslices)

        # Set slice order
        if isinstance(slice_order, str):
            if not interleaved:
                aux = range(nslices)
            else:
                p = nslices/2
                aux = []
                for i in range(p):
                    aux.extend([i,p+i])
                if nslices%2:
                    aux.append(nslices-1)
            if slice_order == 'descending':
                aux.reverse()
            slice_order = aux

        # Set timing values
        self.nslices = nslices
        self.tr = float(tr)
        self.tr_slices = float(tr_slices)
        self.start = float(start)
        self.slice_order = np.asarray(slice_order)
        self.interleaved = bool(interleaved)
        ## assume that the world referential is 'scanner' as defined
        ## by the nifti norm
        self.reversed_slices = to_world[slice_axis][slice_axis]<0

    def z_to_slice(self, z):
        """
        Account for the fact that slices may be stored in reverse
        order wrt the scanner coordinate system convention (slice 0 ==
        bottom of the head)
        """
        if self.reversed_slices:
            return self.nslices - 1 - z
        else:
            return z

    def from_time(self, zv, t):
        """
        tv = from_time(zv, t)
        zv, tv are grid coordinates; t is an actual time value.
        """
        corr = self.tr_slices*interp_slice_order(self.z_to_slice(zv), self.slice_order)
        return (t-self.start-corr)/self.tr

class Realign4d(object):

    def __init__(self,
                 im4d,
                 speedup=_SPEEDUP,
                 optimizer=_OPTIMIZER,
                 transforms=None):
        self.optimizer = optimizer
        dims = im4d.array.shape
        self.dims = dims
        self.nscans = dims[3]
        # Define mask
        speedup = max(1, int(speedup))
        xyz = np.mgrid[0:dims[0]:speedup, 0:dims[1]:speedup, 0:dims[2]:speedup]
        xyz = np.rollaxis(xyz, 0, 4)
        self.xyz = np.reshape(xyz, [np.prod(xyz.shape[0:-1]), 3])
        masksize = self.xyz.shape[0]
        self.data = np.zeros([masksize, self.nscans], dtype='double')
        # Initialize space/time transformation parameters
        self.to_world = im4d.to_world
        self.from_world = np.linalg.inv(self.to_world)
        if transforms == None:
            self.transforms = [Rigid() for scan in range(self.nscans)]
        else:
            self.transforms = transforms
        self.from_time = im4d.from_time
        self.timestamps = im4d.tr*np.arange(self.nscans)
        # Compute the 4d cubic spline transform
        self.cbspline = cspline_transform(im4d.array)

    def resample_inmask(self, t):
        X, Y, Z = grid_coords(self.xyz, self.transforms[t],
                              self.from_world, self.to_world)
        T = self.from_time(Z, self.timestamps[t])
        cspline_sample4d(self.data[:,t], self.cbspline, X, Y, Z, T)

    def resample_all_inmask(self):
        for t in range(self.nscans):
            print('Resampling scan %d/%d' % (t+1, self.nscans))
            self.resample_inmask(t)

    def init_motion_detection(self, t):
        """
        The idea is to compute the global variance using the following
        decomposition:

        V = (n-1)/n V1 + (n-1)/n^2 (x1-m1)^2
          = alpha + beta d2,

        with alpha=(n-1)/n V1, beta = (n-1)/n^2, d2 = (x1-m1)^2.

        Only the second term is variable when one image moves while
        all other images are fixed.
        """
        self.resample_inmask(t)
        fixed = range(self.nscans)
        fixed.remove(t)
        aux = self.data[:, fixed]
        self.m1 = aux.mean(1)
        self.d2 = np.zeros(np.shape(self.m1))
        self.alpha = ((self.nscans-1.0)/self.nscans)*aux.var(1).mean()
        self.beta = (self.nscans-1.0)/self.nscans**2

    def msid(self, t):
        """
        Mean square intensity difference
        """
        self.resample_inmask(t)
        self.d2[:] = self.data[:,t]
        self.d2 -= self.m1
        self.d2 **= 2
        return self.d2.mean()

    def variance(self, t):
        return self.alpha + self.beta*self.msid(t)

    def safe_variance(self, t):
        """
        No need to invoke self.init_motion_detection.
        """
        self.resample_inmask(t)
        self.m = self.data.mean(1)
        self.m2 = (self.data**2).mean(1)
        self.m **= 2
        self.m2 -= self.m
        return self.m2.mean()

    def correct_motion(self):
        optimizer = self.optimizer

        def callback(pc):
            self.transforms[t].param = pc
            print(self.transforms[t])

        if optimizer=='powell':
            tols = {'xtol': _XTOL, 'ftol': _FTOL}
            fmin = fmin_powell
        elif optimizer=='steepest':
            tols = {'xtol': _XTOL, 'ftol': _FTOL, 'step':_STEP}
            fmin = fmin_steepest
        elif optimizer=='cg':
            tols = {'gtol': _GTOL}
            fmin = fmin_cg
        elif optimizer=='bfgs':
            tols = {'gtol': _GTOL}
            fmin = fmin_bfgs
        else: # simplex method
            tols = {'xtol': _XTOL, 'ftol': _FTOL}
            fmin = fmin_simplex

        # Resample data according to the current space/time transformation
        self.resample_all_inmask()

        # Optimize motion parameters
        for t in range(self.nscans):
            print('Correcting motion of scan %d/%d...' % (t+1, self.nscans))
            def loss(pc):
                self.transforms[t].param = pc
                return self.msid(t)
            self.init_motion_detection(t)
            self.transforms[t].param = fmin(loss, self.transforms[t].param,
                                            callback=callback, **tols)

        # At this stage, transforms map an implicit 'ideal' grid to
        # the 'acquisition' grid. We redefine the ideal grid as being
        # conventionally aligned with the first scan.
        T0inv = self.transforms[0].inv()
        for t in range(self.nscans):
            self.transforms[t] = self.transforms[t]*T0inv

    def resample(self):
        print('Gridding...')
        dims = self.dims
        XYZ = np.mgrid[0:dims[0], 0:dims[1], 0:dims[2]]
        XYZ = np.rollaxis(XYZ, 0, 4)
        XYZ = np.reshape(XYZ, [np.prod(XYZ.shape[0:-1]), 3])
        res = np.zeros(dims)
        for t in range(self.nscans):
            print('Fully resampling scan %d/%d' % (t+1, self.nscans))
            X, Y, Z = grid_coords(XYZ, self.transforms[t],
                                  self.from_world, self.to_world)
            T = self.from_time(Z, self.timestamps[t])
            cspline_sample4d(res[:,:,:,t], self.cbspline, X, Y, Z, T)
        return res

def resample4d(im4d, transforms=None):
    """
    corr_im4d_array = resample4d(im4d, transforms=None)
    """
    r = Realign4d(im4d, transforms=transforms)
    return r.resample()

def _realign4d(im4d,
               loops=_WITHIN_LOOPS,
               speedup=_SPEEDUP,
               optimizer=_OPTIMIZER):
    """
    transforms = _realign4d(im4d, loops=2, speedup=4, optimizer='powell')

    Parameters
    ----------
    im4d : Image4d instance

    """
    r = Realign4d(im4d, speedup=speedup, optimizer=optimizer)
    for loop in range(loops):
        r.correct_motion()
    return r.transforms

def realign4d(runs,
              within_loops=_WITHIN_LOOPS,
              between_loops=_BETWEEN_LOOPS,
              speedup=_SPEEDUP,
              optimizer=_OPTIMIZER,
              align_runs=True):
    """

    Parameters
    ----------

    runs : list of Image4d objects

    Returns
    -------
    transforms : list
                 nested list of rigid transformations

    transforms map an 'ideal' 4d grid (conventionally aligned with the
    first scan of the first run) to the 'acquisition' 4d grid for each
    run
    """

    # Single-session case
    if not hasattr(runs, '__iter__'):
        runs = [runs]
    nruns = len(runs)
    if nruns == 1:
        align_runs = False

    # Correct motion and slice timing in each sequence separately
    transforms = [_realign4d(run, loops=within_loops,
                             speedup=speedup, optimizer=optimizer) for run in runs]
    if not align_runs:
        return transforms, transforms, None

    # Correct between-session motion using the mean image of each corrected run
    corr_runs = [resample4d(runs[i], transforms=transforms[i]) for i in range(nruns)]
    aux = np.rollaxis(np.asarray([c.mean(3) for c in corr_runs]), 0, 4)
    ## Fake time series with zero inter-slice time
    ## FIXME: check that all runs have the same to-world transform
    mean_img = Image4d(aux, to_world=runs[0].to_world, tr=1.0, tr_slices=0.0)
    transfo_mean = _realign4d(mean_img, loops=between_loops, speedup=speedup,
                              optimizer=optimizer)

    # Compose transformations for each run
    ctransforms = [None for i in range(nruns)]
    for i in range(nruns):
        ctransforms[i] = [t*transfo_mean[i] for t in transforms[i]]
    return ctransforms, transforms, transfo_mean

max_stars_repo_name: fperez/nipy

max_issues_repo_name: fperez/nipy

max_forks_repo_name: fperez/nipy


text:
#! /usr/bin/env python

# SPDX-FileCopyrightText: Copyright 2021, Siavash Ameli <sameli@berkeley.edu>
# SPDX-License-Identifier: BSD-3-Clause
# SPDX-FileType: SOURCE
#
# This program is free software: you can redistribute it and/or modify it
# under the terms of the license found in the LICENSE.txt file in the root
# directory of this source tree.

# =======
# Imports
# =======

import sys
import time
import numpy
import scipy.sparse
from imate.sample_matrices import band_matrix, band_matrix_logdet
from imate import logdet

# ==============
# relative error
# ==============

def relative_error(estimate, exact):
    """
    Compute the relative error of an estimate, in percent.
    """

    tol = 1e-15
    if numpy.abs(exact) < tol:
        if numpy.abs(estimate - exact) < tol:
            relative_error = 0.0
        else:
            relative_error = numpy.inf
    else:
        relative_error = numpy.abs((estimate - exact) / exact) * 100.0

    return relative_error

# ===================
# test logdet methods
# ===================

def _test_logdet_methods(K, matrix, gram, exponent, assume_matrix):
    """
    Computes the log-determinant of matrix ``K`` with multiple method.

    :param K: Invertible matrix.
    :type K: numpy.ndarray
    """

    # slq method settings
    min_num_samples = 100
    max_num_samples = 200
    lanczos_degree = 30
    error_rtol = 1e-2

    # Use direct method
    time00 = time.time()
    if scipy.sparse.issparse(K):
        logdet0 = numpy.real(numpy.log(
            numpy.linalg.det(K.toarray()).astype(numpy.complex128)))
    else:
        logdet0 = numpy.real(numpy.log(
            numpy.linalg.det(K).astype(numpy.complex128)))
    logdet0 *= exponent
    if gram:
        logdet0 = 2.0 * logdet0
    time01 = time.time()

    # Use eigenvalue method
    time10 = time.time()
    logdet1, _ = logdet(K, method='eigenvalue', gram=gram,
                        assume_matrix=assume_matrix, exponent=exponent,
                        non_zero_eig_fraction=0.95)
    time11 = time.time()

    # Use Cholesky method
    time20 = time.time()
    logdet2, _ = logdet(K, method='cholesky', gram=gram, exponent=exponent,
                        cholmod=None)
    time21 = time.time()

    # Use Stochastic Lanczos Quadrature method
    time30 = time.time()
    logdet3, _ = logdet(K, method='slq', min_num_samples=min_num_samples,
                        max_num_samples=max_num_samples, orthogonalize=-1,
                        lanczos_degree=lanczos_degree, error_rtol=error_rtol,
                        gram=gram, exponent=exponent, verbose=False)
    time31 = time.time()

    # Elapsed times
    elapsed_time0 = time01 - time00
    elapsed_time1 = time11 - time10
    elapsed_time2 = time21 - time20
    elapsed_time3 = time31 - time30

    # Exact solution of logdet for band matrix
    if exponent == 1:
        logdet_exact = band_matrix_logdet(matrix['a'], matrix['b'],
                                          matrix['size'], gram)
        if not gram:
            logdet_exact = 2.0 * logdet_exact
    else:
        logdet_exact = logdet0

    # error
    error0 = relative_error(logdet0, logdet_exact)
    error1 = relative_error(logdet1, logdet_exact)
    error2 = relative_error(logdet2, logdet_exact)
    error3 = relative_error(logdet3, logdet_exact)

    # Print results
    print('')
    print('-------------------------------------------------------------')
    print('Method      Options                     logdet   error   time')
    print('----------  ------------------------  --------  ------  -----')
    print('direct      N/A                       %+8.3f  %5.2f%%  %5.2f'
          % (logdet0, error0, elapsed_time0))
    print('eigenvalue  N/A                       %+8.3f  %5.2f%%  %5.2f'
          % (logdet1, error1, elapsed_time1))
    print('cholesky    N/A                       %+8.3f  %5.2f%%  %5.2f'
          % (logdet2, error2, elapsed_time2))
    print('slq         N/A                       %+8.3f  %5.2f%%  %5.2f'
          % (logdet3, error3, elapsed_time3))
    print('-------------------------------------------------------------')
    print('')

# ===========
# test logdet
# ===========

def test_logdet():
    """
    A test for :mod:`logdetInv.logdet` sub-package.
    """

    matrix = {
        'a': 2.0,
        'b': 1.0,
        'size': 50
    }

    exponents = [0, 1, 2]
    dtypes = [r'float32', r'float64']
    grams = [True, False]
    sparses = [True, False]

    for dtype in dtypes:
        for gram in grams:

            if gram:
                assume_matrix = 'gen'
            else:
                assume_matrix = 'sym'

            # When gram is True:
            #     1. We generate a 2-band nonsymmetric matrix K (hence we set
            #        gram=False in band_matrix).
            #     2. We compute logdet of K.T @ K using only K (hence we set
            #        gram=True in logdet method).
            #
            # When gram is False:
            #     1. We generate a 3-band symmetric matrix K (hence we set
            #        gram=True in band_matrix).
            #     2. We compute logdet of K using K (hence we set
            #        gram=False in logdet method).
            K = band_matrix(matrix['a'], matrix['b'], matrix['size'],
                            gram=(not gram), dtype=dtype)

            for sparse in sparses:
                if not sparse:
                    K = K.toarray()

                for exponent in exponents:
                    print('dtype: %s, ' % (dtype) +
                          'sparse: %5s, ' % (sparse) +
                          'gram: %5s, ' % (gram) +
                          'exponent: %0.4f,\n' % (exponent) +
                          'assume_matrix: %s.' % (assume_matrix))

                    _test_logdet_methods(K, matrix, gram, exponent,
                                         assume_matrix)

# ===========
# System Main
# ===========

if __name__ == "__main__":
    sys.exit(test_logdet())

max_stars_repo_name: ameli/TraceInv

max_issues_repo_name: ameli/TraceInv

max_forks_repo_name: ameli/TraceInv

reason: import numpy,import scipy


text:
import numpy as np
import torch
import torch.nn as nn
from torch.autograd import Variable
import torch.nn.functional as F
from time import time
import numpy as np
#from .pointnet2.pointnet2_modules import PointnetFPModule, PointnetSAModuleMSG
from multi_model.utils.pn2_utils.modules import PointNetSAModule, PointnetFPModule
from multi_model.utils.pn2_utils.nn import SharedMLP

class PointNet2Seg(nn.Module):
    """PointNet++ part segmentation with single-scale grouping

    PointNetSA: PointNet Set Abstraction Layer
    PointNetFP: PointNet Feature Propagation Layer

    Args:
        score_classes (int): the number of grasp score classes
        num_centroids (tuple of int): the numbers of centroids to sample in each set abstraction module
        radius (tuple of float): a tuple of radius to query neighbours in each set abstraction module
        num_neighbours (tuple of int): the numbers of neighbours to query for each centroid
        sa_channels (tuple of tuple of int): the numbers of channels within each set abstraction module
        fp_channels (tuple of tuple of int): the numbers of channels for feature propagation (FP) module
        num_fp_neighbours (tuple of int): the numbers of nearest neighbor used in FP
        seg_channels (tuple of int): the numbers of channels in segmentation mlp
        dropout_prob (float): the probability to dropout input features

    References:
        https://github.com/charlesq34/pointnet2/blob/master/models/pointnet2_part_seg.py

    """
    _SA_MODULE = PointNetSAModule
    _FP_MODULE = PointnetFPModule

    def __init__(self, input_chann = 3, k_score = 1, k_obj = 2, add_channel_flag=False, dropout_prob=0.5):
        super(PointNet2Seg, self).__init__()
        self.k_score = k_score

        num_centroids     = (5120, 1024, 256)
        radius            = (0.02, 0.08, 0.32)
        num_neighbours    = (64, 64, 64)
        sa_channels       = ((128, 128, 256), (256, 256, 512), (512, 512, 1024))
        fp_channels       = ((1024, 1024), (512, 512), (256, 256, 256))
        num_fp_neighbours = (3, 3, 3)
        seg_channels      = (512, 256, 256, 128)

        num_sa_layers = len(num_centroids)
        num_fp_layers = len(fp_channels)

        # Set Abstraction Layers
        feature_channels = input_chann - 3 # 0
        self.sa_modules = nn.ModuleList()
        for ind in range(num_sa_layers):
            sa_module = self._SA_MODULE(in_channels=feature_channels,
                                        mlp_channels=sa_channels[ind],
                                        num_centroids=num_centroids[ind],
                                        radius=radius[ind],
                                        num_neighbours=num_neighbours[ind],
                                        use_xyz=True)
            self.sa_modules.append(sa_module)
            feature_channels = sa_channels[ind][-1]

        inter_channels = [input_chann - 3]# [0]
        inter_channels.extend([x[-1] for x in sa_channels])
        # Feature Propagation Layers
        self.fp_modules = nn.ModuleList()
        feature_channels = inter_channels[-1]
        for ind in range(num_fp_layers):
            fp_module = self._FP_MODULE(in_channels=feature_channels + inter_channels[-2 - ind],
                                        mlp_channels=fp_channels[ind],
                                        num_neighbors=num_fp_neighbours[ind])
            self.fp_modules.append(fp_module)
            feature_channels = fp_channels[ind][-1]

        # MLP
        if not add_channel_flag:
            self.mlp = SharedMLP(feature_channels, seg_channels, ndim=1, dropout_prob=dropout_prob)
        else:
            self.mlp = SharedMLP(feature_channels*3, seg_channels, ndim=1, dropout_prob=dropout_prob)

        self.conv_score = nn.Conv1d(seg_channels[-1], self.k_score, 1)
        self.bn_score = nn.BatchNorm1d(self.k_score)
        self.sigmoid = nn.Sigmoid()

    def forward(self, points, add_channel1=None, add_channel2=None):
        B,C,N = points.size()

        xyz = points[:,:3,:]
        feature = points[:,3:6,:]

        # save intermediate results
        inter_xyz = [xyz]
        inter_feature = [feature]
        # Set Abstraction Layers
        for sa_module in self.sa_modules:
            xyz, feature = sa_module(xyz, feature)
            inter_xyz.append(xyz)
            inter_feature.append(feature)

        # Feature Propagation Layers
        sparse_xyz = xyz
        sparse_feature = feature
        for fp_ind, fp_module in enumerate(self.fp_modules):
            dense_xyz = inter_xyz[-2 - fp_ind]
            dense_feature = inter_feature[-2 - fp_ind]
            fp_feature = fp_module(dense_xyz, sparse_xyz, dense_feature, sparse_feature)
            sparse_xyz = dense_xyz
            sparse_feature = fp_feature

        if add_channel1 is not None and add_channel2 is not None:
            add_channel1 = add_channel1.view(B,1,N).repeat(1,sparse_feature.shape[1],1)
            add_channel2 = add_channel2.view(B,1,N).repeat(1,sparse_feature.shape[1],1)
            sparse_feature = torch.cat((sparse_feature, add_channel1.float(), add_channel2.float()), dim=1)
        # MLP
        x = self.mlp(sparse_feature)
        x_score = self.bn_score(self.conv_score(x))
        x_score = x_score.transpose(2,1).contiguous()
        x_score = self.sigmoid(x_score).view(B, N)

        return sparse_feature, x_score

class PointNet2TwoStage(nn.Module):
    def __init__(self, num_points, input_chann, k_cls, k_reg, k_reg_theta, add_channel_flag=False):
        super(PointNet2TwoStage, self).__init__()
        self.num_points = num_points
        self.k_reg = k_reg
        self.k_cls = k_cls
        self.k_reg_no_anchor = self.k_reg // self.k_cls
        self.k_reg_theta = k_reg_theta

        if not add_channel_flag:
            self.conv = nn.Conv1d(256, 1024, 1)
        else:
            self.conv = nn.Conv1d(256*3, 1024, 1)

        self.bn = nn.BatchNorm1d(1024)

        #self.conv_cls1 = nn.Conv1d(128, 1024, 1)#128128+1024
        self.conv_cls2 = nn.Conv1d(1024, 256, 1)#128+1024
        self.conv_cls3 = nn.Conv1d(256, 128, 1)#128+1024
        self.linear_cls = torch.nn.Linear(128, self.k_cls)
        self.conv_cls4 = nn.Conv1d(128, self.k_cls, 1)
        #self.bn_cls1 = nn.BatchNorm1d(1024)
        self.bn_cls2 = nn.BatchNorm1d(256)
        self.bn_cls3 = nn.BatchNorm1d(128)
        self.bn_cls4 = nn.BatchNorm1d(self.k_cls)

        #self.conv_reg1 = nn.Conv1d(128, 1024, 1)#128+1024
        self.conv_reg2 = nn.Conv1d(1024, 256, 1)#+1024
        self.conv_reg3 = nn.Conv1d(256, 128, 1)#+1024
        self.conv_reg4 = nn.Conv1d(128, self.k_reg, 1)
        #self.bn_reg1 = nn.BatchNorm1d(1024)
        self.bn_reg2 = nn.BatchNorm1d(256)
        self.bn_reg3 = nn.BatchNorm1d(128)
        self.bn_reg4 = nn.BatchNorm1d(self.k_reg)

        #self.conv_reg_theta = nn.Conv1d(int(self.k_reg/7*3), self.k_reg_theta, 1)
        #self.bn_reg_theta = nn.BatchNorm1d(self.k_reg_theta)

        self.mp1 = nn.MaxPool1d(num_points)
        self.ap = torch.nn.AdaptiveAvgPool1d(1)
        self.sigmod = nn.Sigmoid()

    def forward(self, xyz, feature):
        #x = F.avg_pool1d(xyz.float(),self.num_points)
        mp_x = self.mp1(xyz) #[len(true_mask), 128, 1]
        ####mp_x = xyz.view(-1, 128, 1)
        if feature is not None:
            #x[:,:128,:] = x[:,:128,:] + feature.view(feature.shape[0], feature.shape[1],1)
            mp_x = torch.cat((mp_x, feature.view(feature.shape[0], feature.shape[1],1)), dim=1)
            #mp_x = feature.view(feature.shape[0], feature.shape[1],1)
        x = mp_x

        x = F.relu(self.bn(self.conv(x)))

        x_cls = F.relu(self.bn_cls2(self.conv_cls2(x)))
        x_cls = F.relu(self.bn_cls3(self.conv_cls3(x_cls)))
        x_cls = self.bn_cls4(self.conv_cls4(x_cls))

        B,C,_ = x_cls.size()
        x_cls = x_cls.view(B,C)

        x_reg = F.relu(self.bn_reg2(self.conv_reg2(x)))
        x_reg = F.relu(self.bn_reg3(self.conv_reg3(x_reg)))
        x_reg = self.bn_reg4(self.conv_reg4(x_reg))

        x_reg = x_reg.view(B,-1,self.k_reg_no_anchor)
        x_reg[:,:,7:] = self.sigmod(x_reg[:,:,7:])
        '''
        x_reg = x_reg.view(B,-1,7)
        x_reg_theta = self.bn_reg_theta(self.conv_reg_theta(x_reg[:,:,3:6].contiguous().view(B,-1,1)))
        x_reg_new = torch.cat([x_reg, x_reg_theta], dim=2)
        x_reg_new = x_reg_new[:,:,[0,1,2,3,4,5,7,6]]
        return x_cls, x_reg_new, mp_x
        '''
        return x_cls, x_reg, mp_x

class PointNet2Refine(nn.Module):
    def __init__(self, num_points = 2500, input_chann = 3, k_cls = 2, k_reg = 8):
        super(PointNet2Refine, self).__init__()
        self.num_points = num_points
        self.k_reg = k_reg
        self.k_cls = k_cls

        self.conv_formal = nn.Conv1d(384, 1024, 1)
        self.bn_formal = nn.BatchNorm1d(1024)

        # self.conv_formal_cls2 = nn.Conv1d(1024, 128, 1)
        # self.conv_formal_cls3 = nn.Conv1d(128, self.k_cls, 1)
        # self.bn_formal_cls2 = nn.BatchNorm1d(128)
        # self.bn_formal_cls3 = nn.BatchNorm1d(self.k_cls)
        self.conv_formal_cls2 = nn.Conv1d(1024, 128, 1)
        self.conv_formal_cls3 = nn.Conv1d(128, self.k_cls, 1)
        self.bn_formal_cls2 = nn.BatchNorm1d(128)
        self.bn_formal_cls3 = nn.BatchNorm1d(self.k_cls)

        self.conv_formal_reg2 = nn.Conv1d(1024, 128, 1)
        self.conv_formal_reg3 = nn.Conv1d(128, self.k_reg, 1)
        self.bn_formal_reg2 = nn.BatchNorm1d(128)
        self.bn_formal_reg3 = nn.BatchNorm1d(self.k_reg)

        self.mp1 = nn.MaxPool1d(num_points)
        self.ap = torch.nn.AdaptiveAvgPool1d(1)
        self.sigmoid = nn.Sigmoid()

    def forward(self, gripper_feature, group_feature):
        '''
        gripper_feature : [B, region_num, 128]
        group_feature: [B, 128]
        '''
        gripper_feature = self.mp1(gripper_feature) #[B, 128, 1]
        # x: [B, 128, 1]
        x = gripper_feature
        if group_feature is not None:
            # x: [B, 256, 1]
            x = torch.cat((x, group_feature.view(group_feature.shape[0], \
                                group_feature.shape[1],1)), dim=1)
        x = F.relu(self.bn_formal(self.conv_formal(x)))

        # x_cls = F.relu(self.bn_formal_cls2(self.conv_formal_cls2(x)))
        # x_cls = self.bn_formal_cls3(self.conv_formal_cls3(x_cls))
        # x_cls = x_cls.view(x_cls.shape[0], x_cls.shape[1])

        x_cls = F.relu(self.bn_formal_cls2(self.conv_formal_cls2(x)))
        x_cls = self.bn_formal_cls3(self.conv_formal_cls3(x_cls))
        x_cls = x_cls.view(x_cls.shape[0], x_cls.shape[1])

        x_reg = F.relu(self.bn_formal_reg2(self.conv_formal_reg2(x)))
        x_reg = self.bn_formal_reg3(self.conv_formal_reg3(x_reg))
        x_reg = x_reg.view(x_reg.shape[0], x_reg.shape[1])
        #x_reg[:,-1] = self.sigmoid(x_reg[:,-1])

        return x_cls, x_reg
        # return x_reg

if __name__ == '__main__':
    pass


text:
# coding: utf-8
# CS FutureMobility Tool
# See full license in LICENSE.txt.

import numpy as np
import pandas as pd
import openmatrix as omx
import mode_choice.model_defs as md
import mode_choice.matrix_utils as mtx
import mode_choice.mode_choice as mode_choice
import copy
import re
import os
from mode_choice.mc_table_container import table_container

def clean_vehicle_costs(mc_obj, clean_savings):
    ''' Decrease the cost per mile for clean vehicles

    :param mc_obj: the mode choice object containing parameters and inputs
    :param clean_savings: adjustment value (cents per mile) added to op cost
    '''
    mc_obj.cost_per_mile += clean_savings

def land_use_growth_shift(mc_obj, growth_shift_factor, zone_factor_file):
    ''' shifts growth in trips to zones identified in zone factor file

    Growth is shifted to zones identified in the input file from Boston
    zones not included in the file

    :param mc_obj: the mode choice object containing parameters and inputs
    :param growth_shift_factor: amount of growth to shift
    :param zone_factor_file: file containing zones to shift growth to
    '''
    modified_trips = mc_obj.config.scen_path + r"growth_shift_trip_tables.omx"
    factor = growth_shift_factor

    # if file is already in place, just use it, else derive it from baseline
    if os.path.isfile(modified_trips):
        mc_obj.pre_MC_trip_table = mtx.store_omx_as_dict(modified_trips)
    else:
        pre_MC_trip_file_2016 = mc_obj.config.data_path + r"..\2016\pre_MC_trip_6_purposes.omx"
        pre_MC_trip_file_future = mc_obj.config.data_path + r".\pre_MC_trip_6_purposes.omx"

        dense_taz = pd.read_csv(mc_obj.config.param_path + zone_factor_file).sort_values('ID_FOR_CS')[['ID_FOR_CS']]
        dense_taz['dense'] = 1
        all_taz = pd.read_csv(mc_obj.config.taz_file)
        dense_taz_list = all_taz[['ID_FOR_CS','TOWN']].merge(dense_taz,on = 'ID_FOR_CS',how = 'left').fillna(0)[:md.max_zone]
        with omx.open_file(pre_MC_trip_file_2016) as f1 , omx.open_file(pre_MC_trip_file_future) as f2, omx.open_file(modified_trips,'w') as fout:
            for name in f1.list_matrices():
                tt_2016 = np.array(f1[name])[:md.max_zone,:md.max_zone]
                tt_future = np.array(f2[name])[:md.max_zone,:md.max_zone]
                diff = np.maximum(tt_future - tt_2016,0)
                prod_sum_future = tt_future.sum(axis = 1)
                prod_sum_diff = diff.sum(axis = 1)
                prod_sum_diff[(dense_taz_list['TOWN']!='BOSTON,MA')]=0
                sum_to_transfer = prod_sum_diff[(dense_taz_list['TOWN']=='BOSTON,MA')&(dense_taz_list['dense']==0)].sum() * factor
                prod_sum_diff[(dense_taz_list['TOWN']=='BOSTON,MA')&(dense_taz_list['dense']==0)] = (
                    prod_sum_diff[(dense_taz_list['TOWN']=='BOSTON,MA')&(dense_taz_list['dense']==0)] * (factor - 1))
                prod_sum_diff[(dense_taz_list['TOWN']=='BOSTON,MA')&(dense_taz_list['dense']==1)] = (
                        sum_to_transfer * (prod_sum_future[(dense_taz_list['TOWN']=='BOSTON,MA')&(dense_taz_list['dense']==1)]
                        / prod_sum_future[(dense_taz_list['TOWN']=='BOSTON,MA')&(dense_taz_list['dense']==1)].sum()))
                # distribute prod_sum to each attraction zone
                tt_new = pd.DataFrame(tt_future).divide(pd.Series(prod_sum_future), axis = 0).fillna(0).multiply(prod_sum_future + prod_sum_diff,axis = 0).values
                #print('writing {0} trips to table {1}' .format(tt_new.sum(), name))
                fout[name] = tt_new
        mc_obj.pre_MC_trip_table = mtx.store_omx_as_dict(modified_trips)

def compact_land_use(mc_obj, zone_factor_file, ):
    ''' shifts growth in trips according to zones and factors defined in file

    :param mc_obj: the mode choice object containing parameters and inputs
    :param zone_factor_file: file containing zones and growth factor
    '''
    modified_trips = mc_obj.config.scen_path + r"growth_shift_trip_tables.omx"

    pre_MC_trip_file_2016 = mc_obj.config.data_path + r"..\2016\pre_MC_trip_6_purposes.omx"
    pre_MC_trip_file_future = mc_obj.config.data_path + r".\pre_MC_trip_6_purposes.omx"

    zone_factors = pd.read_csv(mc_obj.config.param_path + zone_factor_file).sort_values('ID')

    all_taz = pd.read_csv(mc_obj.data_paths['taz_file'])
    taz_factors = all_taz[['ID','TOWN']].merge(zone_factors,on = 'ID',how = 'left').fillna(1)[:md.max_zone]

# ==DEBUG===========================================================================
#         margs_2016 = {}
#         margs_2040_base = {}
#         margs_2040_new = {}
#
#         margs_2016['ID'] = taz_factors['ID']
#         margs_2016['BOSTON_NB'] = taz_factors['BOSTON_NB']
#         margs_2040_base['ID'] = taz_factors['ID']
#         margs_2040_base['BOSTON_NB'] = taz_factors['BOSTON_NB']
#         margs_2040_new['ID'] = taz_factors['ID']
#         margs_2040_new['BOSTON_NB'] = taz_factors['BOSTON_NB']
# =============================================================================

    with omx.open_file(pre_MC_trip_file_2016) as f1 , omx.open_file(pre_MC_trip_file_future) as f2, omx.open_file(modified_trips,'w') as fout:
        for name in f1.list_matrices():

            tt_2016 = np.array(f1[name])[:md.max_zone,:md.max_zone]
            tt_future = np.array(f2[name])[:md.max_zone,:md.max_zone]

            # only apply to home-based productions
            if (name[:3]!='NHB'):
                diff = np.maximum(tt_future - tt_2016,0)
                adj_diff = np.multiply(diff.T, (np.array(taz_factors['POP_FACTOR']) - 1)).T
                tt_newfuture = tt_future + adj_diff
            else:
                tt_newfuture = tt_future

# =DEBUG============================================================================
#                 margs_2016[name] = tt_2016.sum(axis = 1)
#                 margs_2040_base[name] = tt_future.sum(axis = 1)
#                 margs_2040_new[name] = tt_newfuture.sum(axis = 1)
# =============================================================================
            fout[name] = tt_newfuture

    mc_obj.pre_MC_trip_table = mtx.store_omx_as_dict(modified_trips)
# =DEBUG============================================================================
        #pd.DataFrame(margs_2016).to_csv(mc_obj.config.scen_path + r"base.csv")
        #pd.DataFrame(margs_2040_new).to_csv(mc_obj.config.scen_path + r"newfuture.csv")
        #pd.DataFrame(margs_2040_base).to_csv(mc_obj.config.scen_path + r"origfuture.csv")
# =============================================================================

def transit_ivt_modify_skim(mc_obj, transit_ivt_factor):
    ''' apply transit travel time factor to skims

    Factor transit invehicle time across all modes and time periods

    :param mc_obj: the mode choice object containing parameters and inputs
    :param transit_ivt_factor: factor to be applied to transit ivt
    '''
    mc_obj.DAT_B_skim_PK['Total_IVTT'] = mc_obj.DAT_B_skim_PK['Total_IVTT'] * transit_ivt_factor
    mc_obj.DAT_B_skim_OP['Total_IVTT'] = mc_obj.DAT_B_skim_OP['Total_IVTT'] * transit_ivt_factor
    mc_obj.DAT_CR_skim_PK['Total_IVTT'] = mc_obj.DAT_CR_skim_PK['Total_IVTT'] * transit_ivt_factor
    mc_obj.DAT_CR_skim_OP['Total_IVTT'] = mc_obj.DAT_CR_skim_OP['Total_IVTT'] * transit_ivt_factor
    mc_obj.DAT_RT_skim_PK['Total_IVTT'] = mc_obj.DAT_RT_skim_PK['Total_IVTT'] * transit_ivt_factor
    mc_obj.DAT_RT_skim_OP['Total_IVTT'] = mc_obj.DAT_RT_skim_OP['Total_IVTT'] * transit_ivt_factor
    mc_obj.DAT_LB_skim_PK['Total_IVTT'] = mc_obj.DAT_LB_skim_PK['Total_IVTT'] * transit_ivt_factor
    mc_obj.DAT_LB_skim_OP['Total_IVTT'] = mc_obj.DAT_LB_skim_OP['Total_IVTT'] * transit_ivt_factor
    mc_obj.WAT_skim_PK['Total_IVTT'] = mc_obj.WAT_skim_PK['Total_IVTT'] * transit_ivt_factor
    mc_obj.WAT_skim_OP['Total_IVTT'] = mc_obj.WAT_skim_OP['Total_IVTT'] * transit_ivt_factor

def __transit_time_reduction(skim, idx_list, time_saving_list, factor):
    # apply benefit of new transit projects - factor to account for aggregate zone effect

    for i in range(len(idx_list)):
        skim[tuple(idx_list[i])] -= min(time_saving_list[i],skim[tuple(idx_list[i])] * factor)
    return skim

def transit_improvement_projects(mc_obj, project_impact_file, transit_project_factor):
    ''' reduce transit OVTT and IVTT according to transit project impacts

    Adjust transit times based on transit projects

    :param mc_obj: the mode choice object containing parameters and inputs
    :param project_impact_file: file with transit projects
    :param transit_project_factor: factor to apply to zones with transit projects
    '''
    factor = transit_project_factor
    TAZ_savings = pd.read_csv(mc_obj.config.param_path + project_impact_file)
    ivtt_idx_list = ( TAZ_savings[TAZ_savings['IVTT difference']!=0][['TAZ_0_skim','TAZ_1_skim']].values.tolist()
                    +  TAZ_savings[TAZ_savings['IVTT difference']!=0][['TAZ_1_skim','TAZ_0_skim']].values.tolist())

    ivtt_list = TAZ_savings[TAZ_savings['IVTT difference']!=0]['IVTT difference'].values.tolist() * 2

    ovtt_idx_list = ( TAZ_savings[TAZ_savings['OVTT difference']!=0][['TAZ_0_skim','TAZ_1_skim']].values.tolist()
                    +  TAZ_savings[TAZ_savings['OVTT difference']!=0][['TAZ_1_skim','TAZ_0_skim']].values.tolist())

    ovtt_list = TAZ_savings[TAZ_savings['OVTT difference']!=0]['OVTT difference'].values.tolist() * 2

    mc_obj.DAT_B_skim_PK['Total_IVTT'] = __transit_time_reduction(mc_obj.DAT_B_skim_PK['Total_IVTT'], ivtt_idx_list, ivtt_list,factor)
    mc_obj.DAT_B_skim_OP['Total_IVTT'] = __transit_time_reduction(mc_obj.DAT_B_skim_OP['Total_IVTT'], ivtt_idx_list, ivtt_list,factor)
    mc_obj.DAT_CR_skim_PK['Total_IVTT'] = __transit_time_reduction(mc_obj.DAT_CR_skim_PK['Total_IVTT'], ivtt_idx_list, ivtt_list,factor)
    mc_obj.DAT_CR_skim_OP['Total_IVTT'] = __transit_time_reduction(mc_obj.DAT_CR_skim_OP['Total_IVTT'], ivtt_idx_list, ivtt_list,factor)
    mc_obj.DAT_RT_skim_PK['Total_IVTT'] = __transit_time_reduction(mc_obj.DAT_RT_skim_PK['Total_IVTT'], ivtt_idx_list, ivtt_list,factor)
    mc_obj.DAT_RT_skim_OP['Total_IVTT'] = __transit_time_reduction(mc_obj.DAT_RT_skim_OP['Total_IVTT'], ivtt_idx_list, ivtt_list,factor)
    mc_obj.DAT_LB_skim_PK['Total_IVTT'] = __transit_time_reduction(mc_obj.DAT_LB_skim_PK['Total_IVTT'], ivtt_idx_list, ivtt_list,factor)
    mc_obj.DAT_LB_skim_OP['Total_IVTT'] = __transit_time_reduction(mc_obj.DAT_LB_skim_OP['Total_IVTT'], ivtt_idx_list, ivtt_list,factor)
    mc_obj.WAT_skim_PK['Total_IVTT'] = __transit_time_reduction(mc_obj.WAT_skim_PK['Total_IVTT'], ivtt_idx_list, ivtt_list,factor)
    mc_obj.WAT_skim_OP['Total_IVTT'] = __transit_time_reduction(mc_obj.WAT_skim_PK['Total_IVTT'], ivtt_idx_list, ivtt_list,factor)

    mc_obj.DAT_B_skim_PK['Total_OVTT'] = __transit_time_reduction(mc_obj.DAT_B_skim_PK['Total_OVTT'], ovtt_idx_list, ovtt_list,factor)
    mc_obj.DAT_B_skim_OP['Total_OVTT'] = __transit_time_reduction(mc_obj.DAT_B_skim_OP['Total_OVTT'], ovtt_idx_list, ovtt_list,factor)
    mc_obj.DAT_CR_skim_PK['Total_OVTT'] = __transit_time_reduction(mc_obj.DAT_CR_skim_PK['Total_OVTT'], ovtt_idx_list, ovtt_list,factor)
    mc_obj.DAT_CR_skim_OP['Total_OVTT'] = __transit_time_reduction(mc_obj.DAT_CR_skim_OP['Total_OVTT'], ovtt_idx_list, ovtt_list,factor)
    mc_obj.DAT_RT_skim_PK['Total_OVTT'] = __transit_time_reduction(mc_obj.DAT_RT_skim_PK['Total_OVTT'], ovtt_idx_list, ovtt_list,factor)
    mc_obj.DAT_RT_skim_OP['Total_OVTT'] = __transit_time_reduction(mc_obj.DAT_RT_skim_OP['Total_OVTT'], ovtt_idx_list, ovtt_list,factor)
    mc_obj.DAT_LB_skim_PK['Total_OVTT'] = __transit_time_reduction(mc_obj.DAT_LB_skim_PK['Total_OVTT'], ovtt_idx_list, ovtt_list,factor)
    mc_obj.DAT_LB_skim_OP['Total_OVTT'] = __transit_time_reduction(mc_obj.DAT_LB_skim_OP['Total_OVTT'], ovtt_idx_list, ovtt_list,factor)
    mc_obj.WAT_skim_PK['Total_OVTT'] = __transit_time_reduction(mc_obj.WAT_skim_PK['Total_OVTT'], ovtt_idx_list, ovtt_list,factor)
    mc_obj.WAT_skim_OP['Total_OVTT'] = __transit_time_reduction(mc_obj.WAT_skim_PK['Total_OVTT'], ovtt_idx_list, ovtt_list,factor)

def __bike_time_distance_reduction(time_skim, dist_skim, idx_list, factor_list):
    # change bike time and distance skims by factor
    for i in range(len(idx_list)):
        time_skim[tuple(idx_list[i])] -= time_skim[tuple(idx_list[i])]*factor_list[i]
        dist_skim[tuple(idx_list[i])] -= dist_skim[tuple(idx_list[i])]*factor_list[i]
    return time_skim, dist_skim

def bike_advantage_modify_skim(mc_obj, bike_skim_factor):
    ''' cross the board factoring of bike time and distance

    :param mc_obj: the mode choice object containing parameters and inputs
    :param bike_skim_factor: factor to apply to bike time and distance
    '''
    mc_obj.bike_skim['BikeTime'] = mc_obj.bike_skim['BikeTime'] * bike_skim_factor
    mc_obj.bike_skim['Length (Skim)'] =  mc_obj.bike_skim['Length (Skim)'] * bike_skim_factor

def active_transportation_modify_skim(mc_obj,bike_project_file):
    ''' improve bike and walk skims according to new projects

    :param mc_obj: the mode choice object containing parameters and inputs
    :param bike_project_file: file with zones and factors to adjust bike time and distance
    '''
    TAZ_pair_df = pd.read_csv(mc_obj.config.param_path + bike_project_file)
    idx_list = TAZ_pair_df[['TAZ_0_skim','TAZ_1_skim']].values.tolist()
    factor_list = TAZ_pair_df['factor_avg'].values.tolist()

    mc_obj.bike_skim['BikeTime'], mc_obj.bike_skim['Length (Skim)'] = (
            __bike_time_distance_reduction(mc_obj.bike_skim['BikeTime'],
                                           mc_obj.bike_skim['Length (Skim)'],
                                           idx_list,
                                           factor_list))

    mc_obj.bike_skim['OneMileorLess'] = 1*(mc_obj.bike_skim['Length (Skim)']<=1)

def boston_pev(mc_obj, pev_factor):
    ''' reduce PEV variable by factor for zones within Boston

    :param mc_obj: the mode choice object containing parameters and inputs
    :param pev_factor: adjustments to PEV for all Boston zones
    '''
    mc_obj.AccPEV[0:447,0:447] *= pev_factor
    mc_obj.EgrPEV[0:447,0:447] *= pev_factor

def study_area_pev(mc_obj, pev_factor):
    ''' reduce PEV variable by factor for zones within Boston

    :param mc_obj: the mode choice object containing parameters and inputs
    :param pev_factor: adjustments to PEV for all Boston zones
    '''
    mc_obj.AccPEV[study_area,study_area] *= pev_factor
    mc_obj.EgrPEV[study_area,study_area] *= pev_factor

def parking_cost_change(mc_obj, parking_charge):
    ''' increase parking costs for trips ending in Boston

    :param mc_obj: the mode choice object containing parameters and inputs
    :param parking_charge: value to add to parking charge for all zones in Boston
    '''
    mc_obj.parking[0:md.max_zone,study_area] += parking_charge

def vmt_fee(mc_obj, per_mile_fee):
    ''' adjust vehicle cost per mile to represent VMT fee

    :param mc_obj: the mode choice object containing parameters and inputs
    :param per_mile_fee: value to add to vehicle operating cost
    '''
    mc_obj.cost_per_mile = mc_obj.cost_per_mile + per_mile_fee

def congestion_charge(mc_obj, cordon_area, charge):
    ''' apply a congestion charge for all trips to the cordon area

    :param mc_obj: the mode choice object containing parameters and inputs
    :param cordon_area: cordon area neighborhoods
    :param charge: charge to be applied for trips to the cordon area
    '''
    cong_zones = mc_obj.taz[mc_obj.taz['BOSTON_NB'].isin(cordon_area)]['ID'].values

    cong_charge_table = np.zeros((md.max_zone,md.max_zone))
    cong_charge_table[np.ix_(np.where(np.logical_not(mc_obj.taz['ID'].iloc[:md.max_zone].isin(cong_zones).values))[0],
      np.where(mc_obj.taz['ID'].iloc[:md.max_zone].isin(cong_zones).values)[0])] = charge

    mc_obj.drive_skim_PK['Auto_Toll (Skim)'] += cong_charge_table
    mc_obj.drive_skim_OP['Auto_Toll (Skim)'] += cong_charge_table

def transit_mode_fare_factor(mc_obj, lb_fact, rt_fact, cr_fact, boat_fact, walk_fact):
    ''' adjust transit fares by sub-mode

    :param mc_obj: the mode choice object containing parameters and inputs
    :param lb_fact: drive to local bus fare factor
    :param rt_fact: drive to rapid transit fare factor
    :param cr_fact: drive to commuter rail fare factor
    :param boat_fact: drive to ferry fare factor
    :param walk_fact: walk to any transit mode fare factor
    '''
    mc_obj.DAT_B_skim_PK['Total_Cost'] *= boat_fact
    mc_obj.DAT_CR_skim_PK['Total_Cost'] *= cr_fact
    mc_obj.DAT_RT_skim_PK['Total_Cost'] *= rt_fact
    mc_obj.DAT_LB_skim_PK['Total_Cost'] *= lb_fact
    mc_obj.WAT_skim_PK['Total_Cost'] *= walk_fact

    mc_obj.DAT_B_skim_OP['Total_Cost'] *= boat_fact
    mc_obj.DAT_CR_skim_OP['Total_Cost'] *= cr_fact
    mc_obj.DAT_RT_skim_OP['Total_Cost'] *= rt_fact
    mc_obj.DAT_LB_skim_OP['Total_Cost'] *= lb_fact
    mc_obj.WAT_skim_OP['Total_Cost'] *= walk_fact

def reduce_hbw_peak_trips(mc_obj, tdm_rate):
    ''' reduce HBW trips

    :param mc_obj: the mode choice object containing parameters and inputs
    :param tdm_rate: adjustment to hbw peak trips
    '''
    mc_obj.pre_MC_trip_table['HBW_PK_0Auto'] *= 1-tdm_rate
    mc_obj.pre_MC_trip_table['HBW_PK_wAuto'] *= 1-tdm_rate

def reduce_boston_veh_ownership(mc_obj, boston_veh_decrease):
    ''' shift HH from 1+ to zero vehicles for Boston only

    :param mc_obj: the mode choice object containing parameters and inputs
    :param boston_veh_decrease: factor to shift boston trips from with/vehicle
        to zero vehicle
    '''
    veh_segments = list(filter(re.compile('.*_0Auto').match, list(mc_obj.pre_MC_trip_table.keys())))
    for segment in veh_segments:
        mc_obj.pre_MC_trip_table[segment][0:447,0:md.max_zone] += (
                   mc_obj.pre_MC_trip_table[re.sub('0Auto','wAuto',segment)][0:447,0:md.max_zone]
                   * boston_veh_decrease)
        mc_obj.pre_MC_trip_table[re.sub('0Auto','wAuto',segment)][0:447,0:md.max_zone]  *= (1 - boston_veh_decrease)

def reduce_outside_boston_veh_ownership(mc_obj, veh_decrease):
    ''' shift HH from 1+ to zero vehicles for areas outside Boston

    :param mc_obj: the mode choice object containing parameters and inputs
    :param boston_veh_decrease: factor to shift trips from with/vehicle
        to zero vehicle
    '''
    veh_segments = list(filter(re.compile('.*_0Auto').match, list(mc_obj.pre_MC_trip_table.keys())))
    for segment in veh_segments:
        mc_obj.pre_MC_trip_table[segment][448:md.max_zone,0:md.max_zone] += (
                   mc_obj.pre_MC_trip_table[re.sub('0Auto','wAuto',segment)][448:md.max_zone,0:md.max_zone] * veh_decrease)
        mc_obj.pre_MC_trip_table[re.sub('0Auto','wAuto',segment)][448:md.max_zone,0:md.max_zone] *= (1 - veh_decrease)

def reduce_studyarea_veh_ownership(mc_obj, veh_decrease):
    ''' shift HH from 1+ to zero vehicles for study area only

    :param mc_obj: the mode choice object containing parameters and inputs
    :param bveh_decrease: factor to shift study area trips from with/vehicle
        to zero vehicle
    '''
    veh_segments = list(filter(re.compile('.*_0Auto').match, list(mc_obj.pre_MC_trip_table.keys())))
    for segment in veh_segments:
        mc_obj.pre_MC_trip_table[segment][study_area,0:md.max_zone] += (
                   mc_obj.pre_MC_trip_table[re.sub('0Auto','wAuto',segment)][study_area,0:md.max_zone]
                   * veh_decrease)
        mc_obj.pre_MC_trip_table[re.sub('0Auto','wAuto',segment)][study_area,0:md.max_zone]  *= (1 - veh_decrease)

def reduce_outside_studyarea_veh_ownership(mc_obj, veh_decrease):
    ''' shift HH from 1+ to zero vehicles for areas outside study area

    :param mc_obj: the mode choice object containing parameters and inputs
    :param boston_veh_decrease: factor to shift trips from with/vehicle
        to zero vehicle
    '''
    veh_segments = list(filter(re.compile('.*_0Auto').match, list(mc_obj.pre_MC_trip_table.keys())))
    for segment in veh_segments:
        mc_obj.pre_MC_trip_table[segment][(~study_area),0:md.max_zone] += (
                   mc_obj.pre_MC_trip_table[re.sub('0Auto','wAuto',segment)][(~study_area),0:md.max_zone] * veh_decrease)
        mc_obj.pre_MC_trip_table[re.sub('0Auto','wAuto',segment)][(~study_area),0:md.max_zone] *= (1 - veh_decrease)

def shared_mobility(mc_obj, param_file, coeffs):
    ''' setup shared mobility scenarios

    :param mc_obj: the mode choice object containing parameters and inputs
    :param param_file: parameter file with shared mobility alternatives
    :param coeffs: shared mobility fare, waiting time, and subsidy terms
    '''
    mc_obj.sm = mode_choice.Shared_Mobility(coeffs)
    mc_obj.param_file = mc_obj.config.param_path + param_file

def trip_pricing(mc_obj, trip_price):
    ''' set per-trip pricing

    create a new parameter file with the modified constant values
    to reflect additional trip pricing for auto modes

    :param mc_obj: the mode choice object containing parameters and inputs
    :param trip_price: trip price value
    '''
    new_param_file = mc_obj.config.scen_path + r"param_with_pricing.xlsx"
    writer1 = pd.ExcelWriter(new_param_file,engine = 'openpyxl')

    for purpose in md.purposes:
        param = pd.read_excel(mc_obj.param_file, sheet_name = purpose,index_col = 0)
        for mode in (md.drive_modes + md.smart_mobility_modes):
            if mode in param.index:
                for asc in ['ASC_0_PK','ASC_1_PK','ASC_0_OP','ASC_1_OP']:
                    cost_coeff = param.loc[mode, 'Cost']
                    param.loc[mode,asc] += (trip_price / md.AO_dict[mode]) * cost_coeff
        param.to_excel(writer1, sheet_name = purpose)
    writer1.save()
    mc_obj.param_file = new_param_file

# ====CAV Functionality - needs to be restructured=================================================
#
#     param_out, trip_table_conventional, trip_table_CAV = __CAV_input_generator(mc_obj, switches)
#     mc2 = mode_choice.Mode_Choice(config,run_now = False)
#     mc2.load_input()
#
#     mc_obj.pre_MC_trip_table = copy.deepcopy(trip_table_conventional)
#     mc2.pre_MC_trip_table = copy.deepcopy(trip_table_CAV)
#
#     if (switches.get('CAV_management_policy',False) == False):
#         mc2.param_file = mc_obj.config.param_file_av
#         mc2.cost_per_mile -= 0.05 # should only be applied to CAVs
#
#     print('Running CAV scenario for families with no vehicle or no access to CAV...')
#     mc_obj.run_model()
#     print('Running CAV scenario for families with access to CAV...')
#     mc2.run_model()
#
#     # combine post mode choice trip tables
#     for purpose in ['HBW','HBO','NHB', 'HBSc1','HBSc2','HBSc3']:
#         for peak in ['PK','OP']:
#             for veh_own in ['0','1']:
#                 for mode in mc_obj.table_container.get_table(purpose)[f'{veh_own}_{peak}']:
#                     mc_obj.table_container.get_table(purpose)[f'{veh_own}_{peak}'][mode] += mc2.table_container.get_table(purpose)[f'{veh_own}_{peak}'][mode]
#
# def __CAV_input_generator(mc_obj, switches):
#     # point to AV modified parameter file, split by adoption rate
#     AV_adoption = mc_obj.config.AV_adoption
#
#     # set param for CAV households, alternative baseline and mangement policies
#     if switches.get('smart_mobility',False) == True:
#         if switches.get('smart_mobility_management_policy',False):
#             param_file = mc_obj.config.SM_calib_param_file
#         else:
#             param_file = mc_obj.config.SM_calib_param_mngpol_file
#     else:
#         param_file = mc_obj.config.param_file
#
#     # create two sets of trip tables: one with all 0-veh HHs and 1-x% of the 1+ veh HHs, the other with x% of the 1+ veh HHs
#     trip_table_conventional = copy.deepcopy(mc_obj.pre_MC_trip_table)
#     trip_table_CAV = copy.deepcopy(mc_obj.pre_MC_trip_table)
#
#     hh_0_veh = list(filter(re.compile('.*_0Auto').match, list(trip_table_conventional.keys())))
#     hh_1_veh = list(filter(re.compile('.*_wAuto').match, list(trip_table_conventional.keys())))
#     for segment in hh_1_veh:
#         trip_table_conventional[segment] *= (1-AV_adoption)
#         trip_table_CAV[segment] *= AV_adoption
#
#     for segment in hh_0_veh:
#         trip_table_CAV[segment] *= 0
#
#     return param_file, trip_table_conventional, trip_table_CAV
#
# =============================================================================


text:
import matplotlib #plotting package
import scipy #science package
import numpy as np #this is an array manipulation package
print('Ready to follow the course!!')