text:
'''
    File name: Exam.ipynb
    Author: Gabriel Moreira
    Date last modified: 14/10/2018
    Python Version: 3.6.5
'''
import pandas
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from scipy.stats import t
from matplotlib import gridspec
from sklearn import preprocessing
from sklearn.linear_model import LinearRegression, Lasso
from sklearn.model_selection import KFold, LeaveOneOut
from sklearn.linear_model import LogisticRegression

# Data extration

plt.rcParams['figure.dpi'] = 130
dataframe=pandas.read_csv(r"/data/MedicalData1.csv",sep=';',decimal=',')

listColNames=list(dataframe.columns)

XY=dataframe.values
ColNb_Y=listColNames.index('Disease progression')

Y=XY[:,ColNb_Y].reshape((XY.shape[0],1))   #reshape is to make sure that Y is a column vector
X = np.delete(XY, ColNb_Y, 1)

X_scaled = preprocessing.scale(X)

[n_lines, n_col] = X.shape

listColNames.pop(ColNb_Y)     #to make it contains the column names of X only

for Col in range(len(listColNames)):
    plt.plot(X[:,Col],Y[:],'.')
    plt.xlabel(listColNames[Col])
    plt.ylabel('Disease progression')
    plt.show()

"""
    Correlations between th 18 variables
    and Disease Progression
"""
R2 = []
lr = LinearRegression()
for Col in range(X.shape[1]):
    x = X[:,Col].reshape(XY.shape[0], 1)
    lr.fit(x,Y)
    R2.append(lr.score(x,Y))

markerline, stemlines, baseline = plt.stem(np.array([i for i in range(n_col)]), np.array(R2), '-.')
plt.xticks(np.arange(len(listColNames)), listColNames, rotation='vertical');
plt.ylabel('Correlation (R2)')
plt.title('Correlation with Disease Progression')
plt.setp(baseline, 'color', 'navy', 'linewidth', 1)
plt.setp(stemlines, 'color', 'navy', 'linewidth', 1, 'alpha', 0.5);
plt.setp(markerline, 'color', 'navy', 'markersize', 3, 'linewidth', 1);

Acid1DenIndex = listColNames.index('Acid 1 density')
Acid1Den = X[:,Acid1DenIndex].reshape(XY.shape[0], 1)
lr.fit(Acid1Den, Y[:])

print('> R2 DiseaseProgression(Acid1Density) = ' + str(round(lr.score(Acid1Den, Y[:]),4)))

"""
    Evaluates the prediction stability
"""
kf = KFold(n_splits=4)
eSplit = np.zeros(4)
b0 = np.zeros(4)
b1 = np.zeros(4)
i = 0

for train, test in kf.split(Acid1Den):
    lr.fit(Acid1Den[train], Y[train])
    eSplit[i] = np.mean(abs(np.subtract(Y, lr.predict(Acid1Den))))
    b1[i] = lr.coef_
    b0[i] = lr.intercept_
    i += 1

plt.figure()
plt.plot([i for i in range(1,5)], eSplit, color='navy')
plt.title('eSplit (Y~Acid 1 Density)')
plt.xlabel('i-th Fold')
plt.ylabel('eSplit')
plt.axhline(np.mean(eSplit), color='lightgray')
plt.text(2.7, np.mean(eSplit), 'mean eSplit', color='gray')

print('K-FOLDS (k=4) and linear regression for Disease Progression ~ Acid 1 Density\n')
print('> Mean eSplit = %.2f' % round(np.mean(eSplit),2))
print('> Std Variation (b1) = %.2f' % round(np.sqrt(np.var(b1)),2));

"""
    Uses K-folds with K=4 and a linear regression on
    Biomarker 8 and Disease Progression.
    Prints the R2, Mean eSplit and the std_var(b1).
    Plots the eSplit for each fold.
"""

Biomarker8 = X[:,listColNames.index('Biomarker 8')].reshape([len(Y),1])
lr.fit(Biomarker8[:], Y[:])
print('> R2 Progression(Biomarker 8) = ', str(round(lr.score(Biomarker8[:], Y[:]),4)))

i = 0
for train, test in kf.split(Biomarker8):
    lr.fit(Biomarker8[train], Y[train])
    b0[i] = lr.intercept_
    b1[i] = lr.coef_
    eSplit[i] = np.mean(abs(np.subtract(Y, lr.predict(Biomarker8))))
    i += 1

fig = plt.figure(figsize=(9, 3))
gs = gridspec.GridSpec(1, 2, width_ratios=[1, 1])
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[1])

ax0.plot([i for i in range(1,5)], eSplit, color='navy')
ax0.set_title('eSplit (Y~Biomarker8)')
ax0.set_xlabel('i-th Fold')
ax0.set_ylabel('eSplit')
ax0.axhline(np.mean(eSplit), color='lightgray')
ax0.text(2.7, np.mean(eSplit), 'mean eSplit', color='gray')

print('> Mean eSplit = %.2f' % round(np.mean(eSplit),4))
print('> Std Variation (b1) = %.4f\n' % round(np.sqrt(np.var(b1)),4))

"""
    Uses K-folds with K=4 and a linear regression on
    Pressure1 and Disease Progression.
    Prints the R2, Mean eSplit and the std_var(b1)
    Plots the eSplit for each fold
"""

Pressure1 = X[:,listColNames.index('Pressure 1')].reshape([len(Y),1])
lr.fit(Pressure1[:], Y[:])
print('> R2 Progression(Pressure 1) = ', str(round(lr.score(Pressure1[:], Y[:]),4)))

i = 0
for train, test in kf.split(Pressure1):
    lr.fit(Pressure1[train], Y[train])
    b0[i] = lr.intercept_
    b1[i] = lr.coef_
    eSplit[i] = np.mean(abs(np.subtract(Y, lr.predict(Pressure1))))
    i += 1

ax1.plot([i for i in range(1,5)], eSplit, color='navy')
ax1.set_title('eSplit (Y~Pressure1)')
ax1.set_label('i-th Fold')
ax1.set_ylabel('eSplit')
ax1.set_xlabel('i-th Fold')
ax1.axhline(np.mean(eSplit), color='lightgray')
ax1.text(2.7, np.mean(eSplit), 'mean eSplit', color='gray')

print('> Mean eSplit = %.2f' % round(np.mean(eSplit),4))
print('> Std Variation (b1) = %.2f' % round(np.sqrt(np.var(b1)),4))

"""
    This section has the functions used in the statistic t-test.
"""

def ss(x, y):
    """
        Calculates the residual variance.
    """
    n = len(x)
    lr.fit(x[:], y[:])
    x = x.reshape([n,1])
    y = y.reshape([n,1])
    s = np.sqrt(np.sum(np.power(np.subtract(y, lr.predict(x)), 2)) / (n-2))
    return s

def sx2(x):
    """
        Calculates the variance of the variable X, diving by (n-2)
        instead of n.
    """
    n = len(x)
    x = x.reshape([n,1])
    sx2 = (1.0/(n-2))*np.sum(np.power(np.subtract(x, np.array([np.mean(x)]*n).reshape([n,1])),2))
    return sx2

def testStat(x,y, alpha):
    """
        Computes a statistical test on the hypothesis that Beta-1 is zero.
        Type-I error being equal to alpha.
    """
    n = len(x)
    x = x.reshape([n,1])
    y = y.reshape([n,1])
    lr.fit(x[:], y[:])
    b1 = lr.coef_
    interval = ss(x, y)*np.sqrt((1.0/((n-1)*sx2(x))))*t.ppf(1-alpha/2.0, n-2)
    pvalue = t.sf(abs(b1 / (ss(x, y)*np.sqrt(1/((n-1)*sx2(x))))), n-2)
    print('> H0: Beta1 = 0. Type-I error (alpha) = %.2f.' % alpha)
    print('> p-value = %.4E.' % pvalue)
    print('> b1 = %.5f.' % b1[0])
    print('> Confidence interval (%.f percent) = [%.5f, +%.5f].' % (100-alpha*100,
          -interval, interval))
    if pvalue < alpha:
        print('> H0 rejected. We can assume Beta-1 is different than 0.')
    else:
        print('> We cannot reject H0.')

"""
    t-tests with confidence intervals and p-values to determine
    whether or not we can reject the null hypothesis (H0).
    H0: The slope (Beta1) is null - there is no significative relationship.
"""

print('Biomarker 8:')
testStat(Biomarker8, Y, 0.05)
print('\nPressure 1:')
testStat(Pressure1, Y, 0.05)
print('\nAcid 1 Density:')
testStat(Acid1Den, Y, 0.05)

def hii(x):
    """
        Calculates the for the diagonal H matrix.
        Ŷ = HY
    """
    n = len(x)
    x = x.reshape([n,1])
    x_mean = np.array([np.mean(x)]*n)
    num = np.matmul(np.subtract(x.reshape([n,1]), x_mean.reshape(n,1)),
                    np.subtract(x.reshape([1,n]), x_mean.reshape(1,n)))
    den = np.sum(np.square(np.subtract(x, x_mean)))
    h = np.array([1/n]*n**2).reshape(n,n) + np.divide(num,den)
    h_diag = np.diagonal(h)
    return h_diag

Biomarker5 = np.array(X[:,listColNames.index('Biomarker 5')])
n = len(Biomarker5)
Biomarker5 = Biomarker5.reshape([n ,1])

h_diag = hii(Biomarker5)
plt.figure()
plt.axhline(np.mean(h_diag), color='lightgray', linewidth=0.5)
plt.text(-10, np.mean(h_diag), ' average diag(h)', color='gray')
plt.scatter(Biomarker5, h_diag, alpha=0.6, color='navy')
plt.title('Outlier Identification using H-Matrix')
plt.xlabel('Biomarker 5')
plt.ylabel('diag(H) - Outlier Detection');

"""
    Functions used in calculation of the
    Cooks Distance and to find the outliers.
"""

def residue(x,y):
    """
        Calculates the residues vector resultant from
        a linear model being fitted to the pairs (x,y).
    """
    n = len(x)
    x = x.reshape(n, 1)
    y = y.reshape(n, 1)
    lr.fit(x,y)
    res = np.subtract(y, lr.predict(x))

    return res

def res_i(x,y):
    """
        Calculates the standardised residues vector resultant
        from a linear model being fitted to the pairs (x,y).
    """
    n = len(x)
    x = x.reshape([n,1])
    y = y.reshape([n,1])
    num = residue(x,y)
    den = np.multiply(ss(x,y), np.sqrt(np.subtract(np.ones(n), hii(x))))
    res_std = np.divide(num.reshape([n,1]), den.reshape([n,1]))

    return res_std

def cook(x, y):
    """
        Calculates the cook's distances.
    """
    n = len(x)
    x = x.reshape([n,1])
    y = y.reshape([n,1])
    num = hii(x)
    den = np.multiply(np.ones(n)*2, np.subtract(np.ones(n), hii(x)))
    cook = np.multiply(np.divide(num,den).reshape([n,1]), np.square(res_i(x,y)))

    return cook

def outlierFinder(x, y):
    """
        Determines the outliers in a dataset, using the cook's distance
        and a cutoff criterion of 4/n.
        Nominal_cook is the array containing the data not classified as
        outliers.
        <...>_ind are the indexes for the outliers and the nominal data.
    """
    n = len(x)
    cookd = cook(x,y)
    nominal_cook = cookd[np.where(cookd <= 4/n)]
    outliers_cook = cookd[np.where(cookd > 4/n)]
    nominal_ind = np.where(cookd <= 4/n)
    outlier_ind = np.where(cookd > 4/n)

    return [nominal_cook, outliers_cook, nominal_ind, outlier_ind]

y = cook(Biomarker5, Y)

plt.figure()
plt.axhline(y=0, color='gray', linewidth=0.2)
plt.scatter(Biomarker5, cook(Biomarker5, Y), color='navy', alpha=0.7)
plt.title('Cooks Distance for each observation of Biomarker 5')
plt.xlabel('Biomarker 5')
plt.ylabel('Cook''s Distance');

"""
    Represents, graphically, the outliers and the rest of the
    data in different colors and for their Cooks distance.
    Uses a cutoff of 4/n.
"""

[nominal_cook, outliers_cook, nominal_ind, outliers_ind] = outlierFinder(Biomarker5, Y)
plt.figure()
plt.plot(outliers_cook,np.zeros_like(outliers_cook), 'o', color='orange', alpha=0.8)
plt.plot(nominal_cook,np.zeros_like(nominal_cook), 'o', color='navy', alpha=0.3)
plt.title('Outlier identification using Cooks distance')
plt.xlabel('Cook''s distance')
frame1 = plt.gca()
frame1.axes.get_yaxis().set_visible(False)
plt.legend(['Outliers (4/n cutoff)', 'Rest of the data points', ])
plt.axvline(x=4/n, linewidth=0.2, color='black', linestyle='--')
plt.text(4/n,0, '  cutoff', color='gray', rotation='vertical');

"""
    Linear regression between Disease Progression and Biomarker 5
    Uses Cook's distances to identify Outliers and color-mapping
    to represent them differently.
    Also plots the standardised residues.
"""

# Linear Regression with the identified outliers present
lr.fit(Biomarker5, Y)
R2 = round(lr.score(Biomarker5, Y), 4)
text = 'R2 = ' + str(R2)

# Calculates the residues array
res = res_i(Biomarker5, Y)

# Divides the figure two represent the two plots
fig = plt.figure(figsize=(8, 6))
gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1])
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[1])

# Plots stuff - different colors for outliers and normal data
ax0.plot(Biomarker5[outliers_ind[0]], Y[outliers_ind[0]], 'v', color='orange')
ax0.plot(Biomarker5[nominal_ind[0]], Y[nominal_ind[0]], 'o', color='navy', alpha=0.5)
ax0.plot(Biomarker5, lr.predict(Biomarker5), color='navy', alpha=0.8)
ax0.set_title('Disease Progression vs. Biomarker 5 (complete)')
ax0.set_ylabel('Disease Progression')
ax0.legend(['Outliers (4/n cutoff)','Rest of the data points'], loc='upper left')
ax0.text(x=-2.5, y=16, s=text, color='gray')
ax0.xaxis.set_major_locator(plt.NullLocator())
ax0.xaxis.set_major_formatter(plt.NullFormatter())
markerline1, stemlines1, baseline1 = plt.stem(Biomarker5[outliers_ind[0]],res[outliers_ind[0]], '-.')
markerline2, stemlines2, baseline2 = plt.stem(Biomarker5[nominal_ind[0]],res[nominal_ind[0]], '-.')
ax1.set_xlabel('Biomarker 5')
ax1.set_ylabel('Residues std (r_i)')
plt.setp(baseline1, 'color', 'orange', 'linewidth', 0)
plt.setp(stemlines1, 'color', 'orange', 'linewidth', 1);
plt.setp(markerline1, 'color', 'orange', 'markersize', 3, 'linewidth', 1, 'alpha', 0.5);
plt.setp(baseline2, 'color', 'navy', 'linewidth', 0)
plt.setp(stemlines2, 'color', 'navy', 'linewidth', 1, 'alpha', 0.5);
plt.setp(markerline2, 'color', 'navy', 'markersize', 3, 'linewidth', 1);
plt.axhline(y=0, color='navy')

"""
    Linear regression between Disease Progression and Biomarker 5 without outliers.
    Uses Cook's distances to identify Outliers and color-mapping
    to represent them differently.
    Also plots the standardised residues.
"""

# Remove the Outliers from the dataset
Biomarker5_clean = np.delete(Biomarker5, outliers_ind[0]).reshape([len(Biomarker5)-len(outliers_cook),1])
Y_clean = np.delete(Y, outliers_ind[0]).reshape([len(Y)-len(outliers_cook),1])

# Calculates the residues array
res = res_i(Biomarker5_clean, Y_clean)

# Linear regression performed on the clean dataset
lr.fit(Biomarker5_clean, Y_clean)
R2_clean = round(lr.score(Biomarker5_clean, Y_clean),4)
text = 'R2 = ' + str(R2_clean)

# Divides the figure two represent the two plots
fig = plt.figure(figsize=(8, 6))
gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1])
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[1])

# Plots stuff
ax0.set_title('Disease Progression vs. Biomarker 5 (without outliers)')
ax0.set_ylabel('Disease Progression')
ax0.text(x=-2, y=16, s=text, color='gray')
ax0.plot(Biomarker5_clean, lr.predict(Biomarker5_clean), color='navy', alpha=0.8)
ax0.plot(Biomarker5_clean, Y_clean, 'o', color='navy', alpha=0.5)
ax0.xaxis.set_major_locator(plt.NullLocator())
ax0.xaxis.set_major_formatter(plt.NullFormatter())
ax1.set_xlabel('Biomarker 5')
ax1.set_ylabel('Residues std (r_i)')
markerline, stemlines, baseline = plt.stem(Biomarker5_clean,res, '-.')
plt.setp(baseline, 'color', 'navy', 'linewidth', 0)
plt.setp(stemlines, 'color', 'navy', 'linewidth', 1, 'alpha', 0.5);
plt.setp(markerline, 'color', 'navy', 'markersize', 3, 'linewidth', 1);
plt.axhline(y=0, color='navy');

"""
    Preliminary Study on the influence of the alpha parameter (LASSO)
    has on the correlation coefficient R2 and the number of
    selected / unselected variables,
    i.e., number of null coefficients in the Beta vector.
"""

alpha = [0.05+0.01*i for i in range(1, 101)]
lassoR2 = []
NumEliminatedVars = []

for i in range(100):
    lassoReg = Lasso(alpha = alpha[i])
    lassoReg.fit(X_scaled, Y)
    lassoR2.append(round(lassoReg.score(X_scaled, Y), 2))
    NumEliminatedVars.append(np.size(np.where(lassoReg.coef_ == 0)))

# Plots more stuff
fig, ax1 = plt.subplots()
ax1.set_title('R2 and No. of eliminated variables vs. alpha')
ax1.set_xlabel('alpha (LASSO)')
ax1.set_ylabel('No. of eliminated variables', color='green')
ax1.plot(alpha, NumEliminatedVars, '-', alpha=0.8, color='green')
ax1.tick_params(axis='y', labelcolor='green')
ax2 = ax1.twinx()
ax2.set_ylabel('R2', color='navy')
ax2.plot(alpha, lassoR2, '-', alpha=0.8, color='navy')
ax2.tick_params(axis='y', labelcolor='navy')
fig.tight_layout()
plt.show()

def lassoVarSelect(X, Y, theta):
    """
        Selects the variables using Lasso Regression and for a correlation
        coefficient of theta*R2(Linear Regression).
        Varies the alpha automatically and stops iterating once it reaches the
        desired R2 (Without using cross-validation).
    """

    # First a mulitple linear regression is performed to find the R2 max
    lr = LinearRegression()
    lr.fit(X,Y)
    R2Goal = theta*lr.score(X,Y)
    alpha = [0.05+0.01*i for i in range(1, 101)]

    # Increasing alpha until the R2 goal is achieved
    for i in range(100):
        lassoReg = Lasso(alpha = alpha[i])
        lassoReg.fit(X_scaled, Y)
        lassoR2 = lassoReg.score(X_scaled, Y)
        if  lassoR2 <= R2Goal:
            break

    return [lassoR2, lassoReg.coef_]

# Selecting variables with a R2 cutoff of 95% of the value for linear multiple regression
theta = 0.85 # Adjust here! (85% of the R2 for alpha = 0)
[R2, betaVector] = lassoVarSelect(X_scaled,Y,theta)
numSelectVar = np.size(np.where(betaVector != 0))

print('> R2 = %.4f' % round(R2,4))
print('> Number of selected variables: %d' % numSelectVar)

def optimizeSelection(X,Y, alpha, method):
    """
        Method can be linear or log for logistic regression
        Uses Leave-one-out and Lasso regression while iterating over
        multiple alphas.
        Saves the Betas Coefs using in each alpha-iteration and in
        each LOO-iteration
        Outputs: eSplit, the number of variables selected and the beta
        vector in each iteration
    """

    loo = LeaveOneOut()
    num_itera = len(alpha)
    eSplit = np.empty((num_itera, loo.get_n_splits(X)))
    NumVarSelect = np.empty((num_itera, loo.get_n_splits(X)))
    Beta = np.empty((num_itera, loo.get_n_splits(X), np.shape(X)[1]))

    if method == 'linear':

        for i in range(num_itera):
            lassoReg = Lasso(alpha = alpha[i])
            j = 0
            for train_index, test_index in loo.split(X):
                X_train, X_test = X[train_index], X[test_index]
                Y_train, Y_test = Y[train_index], Y[test_index]
                lassoReg.fit(X_train, Y_train)
                Beta[i][j][:] = lassoReg.coef_
                NumVarSelect[i][j] = np.size(Beta[i][j])-np.size(np.where(Beta[i][j] == 0))
                Y_hat = lassoReg.predict(X_test)
                eSplit[i][j] = abs(Y_hat - Y_test)
                j+=1

    elif method == 'log':

        for i in range(num_itera):
            logReg = LogisticRegression(penalty='l1', C=alpha[i])
            logReg.fit(X,Y.reshape(n_lines,))
            j = 0
            for train_index, test_index in loo.split(X):
                X_train, X_test = X[train_index], X[test_index]
                Y_train, Y_test = Y[train_index], Y[test_index]
                logReg.fit(X_train, Y_train)
                Beta[i][j][:] = logReg.coef_
                NumVarSelect[i][j] = np.size(Beta[i][j])-np.size(np.where(Beta[i][j] == 0))
                Y_hat = logReg.predict(X_test)
                eSplit[i][j] = abs(Y_hat - Y_test)
                j+=1

    return [eSplit, NumVarSelect, Beta]

"""
    Calls the function optimizeSelection using the
    specified alpha vector for the iteration.
"""

alpha = [0.05+0.01*i for i in range(1, 101)]
[eSplit, NumVarSelect, Beta] = optimizeSelection(X_scaled,Y, alpha, method='linear')

"""
    Calculates the average eSplit for each alpha and the
    average number of selected variables for each alpha
"""

eSplitMean = np.mean(eSplit, axis=1)
NumVarSelectMean = np.round(np.mean(NumVarSelect, axis=1))

"""
    Finds the minimum eSplit that verifies number of selected variables <= 3,
    the corresponding alpha and the actual number of selected variables
"""
# find the minimum eSplit with max 3 vars selected
eSplit_min = np.min(eSplitMean[np.where(NumVarSelectMean<=3)])

# find the alpha that generates the lowest eSplit for max 3 vars selected
alpha_min = alpha[np.where(eSplitMean == eSplit_min)[0].tolist()[0]]

# find the number of variables satisfying: number of selected vars <= 3 and lowest eSplit
NumVarSelect_min = NumVarSelectMean[np.where(eSplitMean == eSplit_min)[0].tolist()[0]]

"""
    Shows the results of the optimization graphically
"""

fig, ax1 = plt.subplots()
plt.axvline(x=alpha_min, color='lightgray', linestyle='--')
ax1.axhline(y=eSplit_min, color='lightgray', linestyle='--')
ax1.set_title('eSplit and Num. selected variables vs. alpha (LOO)')
ax1.set_xlabel('alpha (LASSO)')
ax1.set_ylabel('Mean eSplit', color='green')
ax1.plot(alpha, eSplitMean, alpha=0.8, color='green')
ax1.plot(alpha_min, eSplit_min,'o', color='green')
ax1.tick_params(axis='y', labelcolor='green')
ax1.text(alpha_min, eSplit_min,'  eSplit min', color='green')
ax2 = ax1.twinx()
ax2.axhline(y=NumVarSelect_min, color='lightgray', linestyle='--')
ax2.set_ylabel('Num. selected variables (mean)', color='navy')
ax2.plot(alpha, NumVarSelectMean, alpha=0.8, color='navy')
ax2.plot(alpha_min, NumVarSelect_min,'o', color='navy')
ax2.text(alpha_min, NumVarSelect_min,'   3 vars', color='navy')
ax2.tick_params(axis='y', labelcolor='navy')
fig.tight_layout()
plt.show()

"""
    Analysis of the stability of the selected variables.
    Shows graphically which variables are selected in each iteration,
    for the alpha value that produces the minimum eSplit.
"""

plt.figure()

for i in range(len(listColNames)):
    plt.axhline(y=i, color='gray', linewidth=0.1)

for j in range(66):
    selVars = np.where(Beta[np.where(eSplitMean == eSplit_min)[0].tolist()[0]][j] != 0)
    plt.plot(np.ones(np.size(selVars))*j, selVars[0], 'o', color=(0.1, 0.2+0.01*j, 0.15+0.012*j), alpha=0.9-0.007*j)

plt.title('Selected Variables in each iteration of LOO (eSplit-min)')
plt.xlabel('i-th Leave-one-out iteration')
plt.yticks(np.arange(len(listColNames)), listColNames)
plt.text(40,2,'Mean eSplit = ' + str(round(eSplit_min,3)), color='gray');

"""
    Reads and prepares the new data frame
"""

dataframe=pandas.read_csv("./MedicalData2.csv",sep=',',decimal=b'.')

listColNames=list(dataframe.columns)

XY=dataframe.values
ColNb_Y=listColNames.index('Pathology type')

#reshape is to make sure that Y is a column vector
Y=XY[:,ColNb_Y].reshape((XY.shape[0],1))
X = np.delete(XY, ColNb_Y, 1)

X_scaled = preprocessing.scale(X)

#to make it contain the column names of X only
listColNames.pop(ColNb_Y);

# Find the probability of each pathology type for each observation from 1 ... to n=66
logReg = LogisticRegression(penalty='l2')
logReg.fit(X_scaled,Y.reshape(n_lines,))
prob = logReg.predict_proba(X_scaled)

# Plots the probability
width = 0.6
plt.bar(np.arange(0,66), prob[:,0], width, color='navy')
plt.bar(np.arange(0,66), prob[:,1], width, bottom=prob[:,0], color='orange')
plt.title('Probability of Pathology Type being 1 or 2')
orange_patch = mpatches.Patch(color='orange', label='Pathology type 2')
blue_patch = mpatches.Patch(color='navy', label='Pathology type 1')
plt.legend(handles=[orange_patch, blue_patch], loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlabel('Observations')
plt.ylabel('Probability')
plt.show()

"""
    Plots the predictions and calculates the prediction
    error to implement color-mapping (red and green) to
    distinguish between good and bad predictions.
"""
# calculates the error in each prediction
err = np.abs(np.subtract(logReg.predict(X_scaled).reshape(-1,1), Y))

# marker size with a radius proportional to the prediction probability
s = []
for i in range(n_lines):
    if prob[i,1] >= 0.5:
        s.append(60*prob[i,1])
    else:
        s.append(60*prob[i,0])

# defines a color vector, green when error = 0 and red when error = 1
color = ['green', 'red']
color = [color[int(i)] for i in err]

# Plots predictions with color mapping for correct and incorrect
plt.scatter(np.arange(66), logReg.predict(X_scaled), s=s, c=color, alpha=0.6)
plt.xlabel('Observations')
plt.ylabel('Pathology type prediction')
red_patch = mpatches.Patch(color='red', label='Incorrect prediction')
blue_patch = mpatches.Patch(color='green', label='Correct prediction')
plt.legend(handles=[red_patch, blue_patch]);

numErrors = np.sum(err == 1)
errorRate = (numErrors / n_lines)
print('> Logistic Regression prediction completed with %d errors.' % numErrors)
print('> Corresponds to an error rate of %.2f.' % errorRate)

"""
    Performs a logistic regression using L1 penalty
    and plots the predicted probability of the
    pathology being type 1 or 2.
"""

# Logistic Regression for multiple penalty factors. The lower the alpha, the greater the regularisation
alpha = [0.001+0.002*i for i in range(1, 100)]
[eSplit, NumVarSelect, Beta] = optimizeSelection(X_scaled,Y.reshape(n_lines,), alpha, method='log')

"""
    Calculates the average eSplit for each alpha and the
    average number of selected variables for each alpha.
"""

eSplitMean = np.mean(eSplit, axis=1)
NumVarSelectMean = np.round(np.mean(NumVarSelect, axis=1))

"""
    Finds the minimum eSplit that verifies number of selected variables = 1,
    the corresponding alpha and the actual number of selected variables.
"""

# find the minimum eSplit with max 1 vars selected
eSplit_min = np.min(eSplitMean[np.where(NumVarSelectMean==1)])

# find the alpha that generates the lowest eSplit for 1 var selected
alpha_min = alpha[np.where(eSplitMean == eSplit_min)[0].tolist()[0]]

# find the number of variables satisfying: number of selected vars = 1 and lowest eSplit
NumVarSelect_min = NumVarSelectMean[np.where(eSplitMean == eSplit_min)[0].tolist()[0]]

"""
    Shows the results of the optimization graphically.
"""

fig, ax1 = plt.subplots()
plt.axvline(x=alpha_min, color='lightgray', linestyle='--')
ax1.axhline(y=eSplit_min, color='lightgray', linestyle='--')
ax1.set_title('eSplit and Num. selected variables vs. alpha (LOO)')
ax1.set_xlabel('alpha (LogReg)')
ax1.set_ylabel('Mean eSplit', color='green')
ax1.plot(alpha, eSplitMean, alpha=0.8, color='green')
ax1.plot(alpha_min, eSplit_min,'o', color='green')
ax1.tick_params(axis='y', labelcolor='green')
ax1.text(alpha_min, eSplit_min,'  eSplit min', color='green')
ax2 = ax1.twinx()
ax2.axhline(y=NumVarSelect_min, color='lightgray', linestyle='--')
ax2.set_ylabel('Num. selected variables (mean)', color='navy')
ax2.plot(alpha, NumVarSelectMean, alpha=0.8, color='navy')
ax2.plot(alpha_min, NumVarSelect_min,'o', color='navy')
ax2.text(alpha_min, NumVarSelect_min,'   1 var', color='navy')
ax2.tick_params(axis='y', labelcolor='navy')
fig.tight_layout()
plt.show()

"""
    Analysis of the stability of the selected variables.
    Shows graphically which variables are selected in each iteration,
    for the alpha value that produces the minimum eSplit.
"""

plt.figure()

for i in range(len(listColNames)):
    plt.axhline(y=i, color='gray', linewidth=0.1)

for j in range(66):
    selVars = np.where(Beta[np.where(eSplitMean == eSplit_min)[0].tolist()[0]][j] != 0)
    plt.plot(np.ones(np.size(selVars))*j, selVars[0], 'o', color=(0.1, 0.2+0.01*j, 0.15+0.012*j), alpha=0.9-0.007*j)

plt.title('Selected Variables in each iteration of LOO (eSplit-min)')
plt.xlabel('i-th Leave-one-out iteration')
plt.yticks(np.arange(len(listColNames)), listColNames)
plt.text(40,2,'Mean eSplit = ' + str(round(eSplit_min,3)), color='gray');

"""
    Plots every variable 1...18 with color mapping to Pathology
    type 1 and Pathology type 2 to identify if Cells 2 density
    is in fact linked to the pathology type.
"""

f, ax = plt.subplots(nrows = 6, ncols = 3)
f.tight_layout()

green_patch = mpatches.Patch(color='green', label='Pathology Type 1')
blue_patch = mpatches.Patch(color='navy', label='Pathology Type 2')
plt.legend(handles=[green_patch, blue_patch], loc='center left', bbox_to_anchor=(1.2, 0.5));

for index in range(n_col):
    ColNb_Var = index
    c = []
    for i in range(n_lines):
        if Y[i] == 1:
            c.append('navy')
        else:
            c.append('green')
    j = int(ColNb_Var % 3)
    i = int(np.trunc(ColNb_Var / 3))
    ax[i][j].scatter(X[:,ColNb_Var].reshape(1,-1), np.zeros_like(X[:,ColNb_Var].reshape(1,-1)), c=c, alpha=0.6)
    ax[i][j].set_xlabel(listColNames[ColNb_Var])
    ax[i][j].get_yaxis().set_ticks([])

"""
    Uses just 'Cells 2 density' to
    perform the logistic regression.
"""

logReg = LogisticRegression()
logReg.fit(X_scaled[:,listColNames.index('Cells 2 density')].reshape(-1, 1),Y.reshape(n_lines,))
R = logReg.score(X_scaled[:,listColNames.index('Cells 2 density')].reshape(-1, 1),Y.reshape(n_lines,))
print('> Fit score = %.2f' % round(R,2))

"""
    Plots the predictions and calculates the prediction
    error to implement color-mapping (red and green) to
    distringuish between good and bad predictions.
"""
# calculates the error in each prediction
err = np.abs(np.subtract(logReg.predict(X_scaled[:,listColNames.index('Cells 2 density')].reshape(-1, 1)), Y[:,0]))
print(type(logReg.predict(X_scaled[:,listColNames.index('Cells 2 density')].reshape(-1, 1))))

# marker size with a radius proportional to the prediction probability
s = []
for i in range(n_lines):
    if prob[i,1] >= 0.5:
        s.append(60*prob[i,1])
    else:
        s.append(60*prob[i,0])

# defines a color vector, green when error = 0 and red when error = 1
color = ['green', 'red']
color = [color[int(i)] for i in err]

# Plots predictions with color mapping for correct and incorrect
plt.scatter(np.arange(66), logReg.predict(X_scaled[:,listColNames.index('Cells 2 density')].reshape(-1, 1)), s=s, c=color, alpha=0.6)
plt.xlabel('Observations')
plt.ylabel('Pathology type prediction')
plt.title('Pathology type prediction using Cells 2 density')
red_patch = mpatches.Patch(color='red', label='Incorrect prediction')
green_patch = mpatches.Patch(color='green', label='Correct prediction')
plt.legend(handles=[red_patch, green_patch]);

max_stars_repo_name: gabmoreira/datascience

max_issues_repo_name: gabmoreira/datascience

max_forks_repo_name: gabmoreira/datascience


text:
import argparse

from scipy.optimize import differential_evolution

from sklearn.ensemble import RandomForestClassifier
from imblearn.metrics import geometric_mean_score

import numpy as np
import pickle

with open('../X_train.pickle', 'rb') as f:
    X_train = pickle.load(f)

with open('../y_train.pickle', 'rb') as f:
    y_train = pickle.load(f)

with open('../X_test.pickle', 'rb') as f:
    X_test = pickle.load(f)

with open('../y_test.pickle', 'rb') as f:
    y_test = pickle.load(f)

__author__ = "Manolomon"
__license__ = "MIT"
__version__ = "1.0"

import pickle

def logger(xa, convergence):
    x_str = np.array_repr(xa).replace('\n', '')
    print(x_str)

def fitness_func(individual):  # Fitness Function
    global X_train
    global y_train
    classifier = RandomForestClassifier(
        n_estimators=int(round(individual[0])),
        max_depth=None if individual[1] < 2 else int(round(individual[1])),
        min_samples_split=individual[2],
        min_samples_leaf=individual[3],
        max_features="auto" if individual[4] < 0.5 else "sqrt",
        bootstrap=False if individual[5] < 0.5 else True,
        random_state=int(round(individual[6]))
    )
    classifier.fit(X_train, y_train)

    g_mean = geometric_mean_score(
        y_test, classifier.predict(X_test), average='weighted')
    del classifier
    return -1 * g_mean

bounds = [  # Parameters tuned in https://github.com/miguelfzafra/Latest-News-Classifier/blob/master/0.%20Latest%20News%20Classifier/09.%20Report/Latest%20News%20Classifier.pdf plus random_state
    (1, 100),  # n_estimators: int, default=100
    (1, 10),  # max_depth: int, default=None
    (0, 1),  # min_samples_split: int or float, default=2
    (0, 0.5),  # min_samples_leaf: int or float, default=1
    (0, 1),  # max_features: {"auto", "sqrt", "log2"}, int or float, default="auto"
    (0, 1),  # bootstrap: bool, default=True
    (0, 10)  # random_state: int or RandomState, default=None
]

if __name__ == "__main__":
    ap = argparse.ArgumentParser(
        description='Random Forest Hyperparameter tuning for software requirements categorization using Differential Evolution')
    ap.add_argument("-v", "--verbose",
                    help="increase output verbosity", action="store_true")
    ap.add_argument('--np', dest='np', type=int,
                    required=True, help='Population size')
    ap.add_argument('--max_gen', dest='max_gen', type=int,
                    required=True, help='Genarations')
    ap.add_argument('--f', dest='f', type=float,
                    required=True, help='Scale Factor')
    ap.add_argument('--cr', dest='cr', type=float,
                    required=True, help='Crossover percentage')
    ap.add_argument('--datfile', dest='datfile', type=str,
                    help='File where it will be save the score (result)')

    args = ap.parse_args()

    result = differential_evolution(fitness_func, bounds, disp=True, popsize=args.np, maxiter=args.max_gen,
                                    mutation=args.f, recombination=args.cr, strategy='rand1bin', callback=logger)

    # print("Best individual: [n_estimators=%s, max_depth=%s, min_samples_split=%s, min_samples_leaf=%s, max_features=%s, bootstrap=%s, random_state=%s] f(x)=%s" % (
    #    int(round(result.x[0])),
    #    None if result.x[1] < 2 else int(round(result.x[1])),
    #    result.x[2],
    #    result.x[3],
    #    "auto" if result.x[4] < 0.5 else "sqrt",
    #    False if result.x[5] < 0.5 else True,
    #    int(round(result.x[6])),
    #    result.fun*(-100)))

    if args.datfile:
        with open(args.datfile, 'w') as f:
            f.write(str(result.fun*(-100)))

ext: py

max_stars_repo_name: quality-attributes/autoML-models

max_issues_repo_name: quality-attributes/autoML-models

max_forks_repo_name: quality-attributes/autoML-models

reason: import numpy,from scipy


text:
#!/usr/bin/python
#
# Plots the results from the 2D pose graph optimization. It will draw a line
# between consecutive vertices.  The commandline expects two optional filenames:
#
#   ./plot_results.py --initial_poses optional --optimized_poses optional
#
# The files have the following format:
#   ID x y yaw_radians

import matplotlib.pyplot as plot
import numpy
import sys
from optparse import OptionParser

parser = OptionParser()
parser.add_option("--initial_poses", dest="initial_poses",
                  default="", help="The filename that contains the original poses.")
parser.add_option("--optimized_poses", dest="optimized_poses",
                  default="", help="The filename that contains the optimized poses.")
(options, args) = parser.parse_args()

# Read the original and optimized poses files.
poses_original = None
if options.initial_poses != '':
  poses_original = numpy.genfromtxt(options.initial_poses, usecols = (1, 2))

poses_optimized = None
if options.optimized_poses != '':
  poses_optimized = numpy.genfromtxt(options.optimized_poses, usecols = (1, 2))

# Plots the results for the specified poses.
plot.figure()
if poses_original is not None:
  plot.plot(poses_original[:, 0], poses_original[:, 1], '-', label="Original",
            alpha=0.5, color="green")

if poses_optimized is not None:
  plot.plot(poses_optimized[:, 0], poses_optimized[:, 1], '-', label="Optimized",
            alpha=0.5, color="blue")

plot.axis('equal')
plot.legend()
# Show the plot and wait for the user to close.
plot.show()

reason: import numpy


text:
import numpy as np
import pandas as pd
from sklearn.ensemble import ExtraTreesClassifier, RandomForestClassifier
from sklearn.feature_selection import RFE, SelectPercentile, f_classif
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline, make_union
from sklearn.preprocessing import MinMaxScaler, PolynomialFeatures
from sklearn.svm import LinearSVC
from tpot.builtins import StackingEstimator
from sklearn.preprocessing import FunctionTransformer
from copy import copy

# NOTE: Make sure that the outcome column is labeled 'target' in the data file
tpot_data = pd.read_csv('PATH/TO/DATA/FILE', sep='COLUMN_SEPARATOR', dtype=np.float64)
features = tpot_data.drop('target', axis=1)
training_features, testing_features, training_target, testing_target = \
            train_test_split(features, tpot_data['target'], random_state=None)

# Average CV score on the training set was: 0.8301746327045152
exported_pipeline = make_pipeline(
    make_union(
        make_pipeline(
            SelectPercentile(score_func=f_classif, percentile=65),
            PolynomialFeatures(degree=2, include_bias=False, interaction_only=False),
            StackingEstimator(estimator=RandomForestClassifier(bootstrap=True, criterion="gini", max_features=0.4, min_samples_leaf=17, min_samples_split=6, n_estimators=100)),
            MinMaxScaler()
        ),
        FunctionTransformer(copy)
    ),
    RFE(estimator=ExtraTreesClassifier(criterion="gini", max_features=0.7000000000000001, n_estimators=100), step=0.8),
    LinearSVC(C=0.1, dual=True, loss="squared_hinge", penalty="l2", tol=0.001)
)

exported_pipeline.fit(training_features, training_target)
results = exported_pipeline.predict(testing_features)


text:
# This is the script for cutting out clips with shots, no-shots and attacks and putting them into data_dir

from config import *

import mPyPl as mp
from mpyplx import *
import os
from pipe import Pipe
from moviepy.editor import *
import itertools
import numpy as np
import json
import argparse
from distutils import util

def convtime(s):
    t = s.split(':') | mp.select(float) | mp.as_list
    return t[2]+t[1]*60+t[0]*60*60

parser = argparse.ArgumentParser()
parser.add_argument('--debugMatchId', type=int, help='Produce clips for the specified match only', required=False)
parser.add_argument('--debugMaxClips', type=int, help='Produce up to N clips for the specified match/half time only', required=False)
parser.add_argument('--debug', help='Enable detailed traces', required=False, default='False')
parser.add_argument('--debugShotsOnly', help='Produce Shots clip only', required=False, default='False')
parser.add_argument('--goalClipEnabled', help='Toggle "Goal" Shot clips', required=False, default='True')
parser.add_argument('--overwriteExisting', help='If a clip exists, overwrite or not', required=False, default='False')
parser.add_argument('--clipTrimTime', help='Specify half-clip duration', required=False, default=2.5)
parser.add_argument('--clipTrimAttackTimeRange', help='Specify +/- X seconds centered on event time', required=False, default=10)
parser.add_argument('--clipTrimAttackStart', help='Specify when attack clip starts (before event time) in seconds', required=False, default=1)
parser.add_argument('--clipTrimAttackEnd', help='Specify when attack clip ends (after event time) in seconds', required=False, default=4)
parser.add_argument('--noShotClipDuration', help='Specify no-shot clip duration in seconds', required=False, default=5)
parser.add_argument('--noShotClipInterval', help='Specify how many no-shot to take between 2 shot events', required=False, default=10)
parser.add_argument('--resizeClips', help='Toggle clip resizing', required=False, default='False')
args = parser.parse_args()

shotFolder = "shot"
noShotFolder = "noshot"
attacksFolder = "attack"

debug = util.strtobool(args.debug)
shotsOnly = util.strtobool(args.debugShotsOnly)
overwriteExist = util.strtobool(args.overwriteExisting)
goalClipEnabled = util.strtobool(args.goalClipEnabled)
resizeClips = util.strtobool(args.resizeClips)

shotClips = []
noshotClips = []
attackClips = []

print(f'         Debug: {debug}')
print(f'     ShotsOnly: {shotsOnly}')
print(f'   GoalEnabled: {goalClipEnabled}')
print(f'     Overwrite: {overwriteExist}')
print(f'Clip Trim-time: {args.clipTrimTime}')

def shotFilter(tags):
    if goalClipEnabled:
        return "shot" in tags or "shot_goal" in tags
    else:
        return "shot" in tags

def attackFilter(tags):
    return "attack" in tags

def addClipToList(list, matchId, n, h, eventMatchPhase, x, mt, dt, satoff, startClip, endClip, isGoal, goalOff, shotTypes, prev_time, next_time):
    clipMeta = {"matchId": matchId,
                "n": n,
                "Half": h,
                "EventMatchPhase": eventMatchPhase,
                "MatchTime": x,
                "PhaseStartTime": mt,
                "Kickoff": dt,
                "SatOff": satoff,
                "GoalOff": goalOff,
                "ClipStart": startClip,
                "ClipEnd": endClip,
                "IsGoal": isGoal,
                "Type": shotTypes,
                "Previous_time": prev_time,
                "Next_time": next_time
        }
    list.append(clipMeta)

# Convert back seconds to HH:mm:ss, for easily jump in marks.jsonl
def GetTime(seconds):
    hh = seconds//(60*60)
    mm = (seconds-hh*60*60)//60
    ss = seconds-(hh*60*60)-(mm*60)
    return (f"{int(hh):0>2d}:{int(mm):0>2d}:{int(ss):0>2d}")

matches = from_json(os.path.join(source_dir,'matches.json'))

for match in matches:

    if args.debugMatchId != None:
        if match["Id"] != args.debugMatchId:
            #print(f"Skipping match {match['Id']} - Debug mode for match {args.debugMatchId}")
            continue
        else:
            print(f"\nDebug mode for match '{args.debugMatchId}'")

    dt = convtime(match["ClippingSettings"]["Kickoff"])
    satoff = match["ClippingSettings"]["SatOffset"]
    SecondHalfDelay = match["ClippingSettings"]["HalfDelay"]
    matchDay = match["MatchDay"]
    goaloff = match["ClippingSettings"]["ActionOverrides"]["Goal"]["ActionDelay"]

    competition = match["Competition"].lower()
    season = match["Season"].replace("/", "")

    h = max(0, int(match["Half"]))

    print(f"Processing {match['Video']} ({match['Id']})\n")
    videoFilePath = os.path.join(source_dir,match["Video"])
    marksFilePath = os.path.join(source_dir,"Marks.jsonl")
    video = VideoFileClip(videoFilePath)

    data = []
    with open(marksFilePath,'r') as f:
        for cnt, line in enumerate(f):
            data.append(json.loads(line))

    # Get list of times when match halves start
    mt = data | mp.filter('eventType', lambda x: "start" in x) | mp.select_field('matchTime') | mp.as_list
    mt = convtime(mt[h])

    # Shots are fine, added also Goal which are not marked as shots
    cuts = data | mp.filter('eventType', shotFilter) | mp.apply('matchTime', 'start', convtime)

    # Consider, for negative examples, also Attacks, which are Shots nearby the goal area (or in other places of the field)
    # They are also filtered (later), removing those overlapping with shots/goals
    attacks = data | mp.filter('eventType', attackFilter) | mp.apply('matchTime', 'start', convtime)

    clipTrimTime = float(args.clipTrimTime)                       # +/- X seconds centered on event time

    clipTrimAttackTimeRange = float(args.clipTrimAttackTimeRange) # +/- X seconds centered on event time
    clipTrimAttackStart = float(args.clipTrimAttackStart)         # when attack clip starts (before event time)
    clipTrimAttackEnd = float(args.clipTrimAttackEnd)             # when attack clip ends (after event time)

    noShotClipDuration = float(args.noShotClipDuration)
    noShotClipInterval = float(args.noShotClipInterval)

    n=0
    prev_time = mt
    for cutInfo in cuts:

        if args.debugMaxClips != None and n >= args.debugMaxClips:
            print(f"Debug mode - Reached max number of clips ({args.debugMaxClips})\n")
            break

        x = cutInfo['start']

        eventMatchPhase = cutInfo["phase"]

        # write out correct "shot" clip

        isGoal = 'shot_goal' in cutInfo['eventType']
        shotType = cutInfo['eventType']

        additionalShiftForGoal = 0
        shotKindForFileName = "Shot"
        if isGoal:
            shotKindForFileName = "Goal"
            additionalShiftForGoal = goaloff

        startClip = x + dt - mt - satoff - clipTrimTime - additionalShiftForGoal
        endClip =  x + dt - mt - satoff + clipTrimTime - additionalShiftForGoal

        shotTypeFolder = shotType
        shotTypesForFileName = shotType

        # Skip cuts that are in the wrong "half" of the video
        if eventMatchPhase != h:
            continue

        if eventMatchPhase==0 and h == eventMatchPhase and startClip > video.duration:
            print(f"Cut in the wrong 'half' of the video (h={h+1}). Skipping\n")
            break

        if eventMatchPhase==1 and h == eventMatchPhase and startClip <= 0:
            print(f"Cut in the wrong 'half' of the video (h={h+1}). Skipping\n")
            continue

        if debug:
            print("Shot\n====================================================")
            print(f"n:               {n}\n" +
                  f"Half:            {h+1}\n" +
                  f"EventMatchPhase: {eventMatchPhase+1}\n" +
                  f"MatchTime:       {GetTime(x)}\n" +
                  f"PhaseStartTime:  {GetTime(mt)}\n" +
                  f"Kickoff:         {GetTime(dt)}\n" +
                  f"SatOff:          {GetTime(satoff)}\n" +
                  f"ClipStart:       {GetTime(startClip)}\n" +
                  f"ClipEnd:         {GetTime(endClip)}\n" +
                  f"IsGoal:          {isGoal}\n" +
                  f"Type:            {shotTypesForFileName}\n")

        if startClip < 0:
            print(f"StartClip < 0. Forcing to 0.0\n")
            startClip = 0

        if endClip > video.duration:
            print(f"EndClip > video.duration. Forcing to video.duration\n")
            endClip = video.duration

        addClipToList(shotClips, match["Id"], n, h+1, eventMatchPhase+1, GetTime(x), GetTime(mt), GetTime(dt), GetTime(satoff), GetTime(startClip), GetTime(endClip), isGoal, goaloff, shotType, "", "")

        cut = video.subclip(startClip, endClip)
        outputDir = os.path.join(base_dir, "processed", shotFolder, shotTypeFolder)
        os.makedirs(outputDir, exist_ok=True)
        videoOutFileName = os.path.join(outputDir,f"{match['Id']}_{h+1}_{n}_{shotTypesForFileName}.full.mp4")
        if not os.path.isfile(videoOutFileName) or overwriteExist == True: # skip if exists
            cut.write_videofile(videoOutFileName)

        if resizeClips:
            resizedOutputFile = videoOutFileName.replace(".full.mp4",".resized.mp4")
            resized = cut.fx(vfx.resize, width=video_size[0])
            if not os.path.isfile(resizedOutputFile) or overwriteExist == True: # skip if exists
                resized.write_videofile(resizedOutputFile)

        # Skips No-Shots in DEBUG
        if shotsOnly == False:

            # Write out "noshot" clip (in between this and last shot)

            nx = (prev_time+x)/2
            if x - nx > noShotClipInterval:

                startClip = nx + dt - mt + satoff
                endClip = nx + dt - mt + satoff + noShotClipDuration

                # Skip cuts that are in the wrong "half" of the video
                if eventMatchPhase != h:
                    continue

                if eventMatchPhase==0 and h == eventMatchPhase and startClip > video.duration:
                    #print(f"Cut in the wrong 'half' of the video (h={h+1}). Skipping\n")
                    break
                if eventMatchPhase==1 and h == eventMatchPhase and startClip <= 0:
                    #print(f"Cut in the wrong 'half' of the video (h={h+1}). Skipping\n")
                    continue

                if debug:
                    print("No-shot\n====================================================")
                    print(f"n:               {n}\n" +
                          f"Half:            {h+1}\n" +
                          f"EventMatchPhase: {eventMatchPhase+1}\n" +
                          f"Prev:            {GetTime(prev_time)}\n" +
                          f"Next:            {GetTime(nx)}\n" +
                          f"PhaseStartTime:  {GetTime(mt)}\n" +
                          f"Kickoff:         {GetTime(dt)}\n" +
                          f"SatOff:          {GetTime(satoff)}\n" +
                          f"ClipStart:       {GetTime(startClip)}\n" +
                          f"ClipEnd:         {GetTime(endClip)}\n")

                if startClip < 0:
                    print(f"StartClip < 0. Forcing to 0.0\n")
                    startClip = 0

                if endClip > video.duration:
                    print(f"EndClip > video.duration. Forcing to video.duration\n")
                    endClip = video.duration

                addClipToList(noshotClips, match["Id"], n, h+1, eventMatchPhase+1, GetTime(x), GetTime(mt), GetTime(dt), GetTime(satoff), GetTime(startClip), GetTime(endClip), "false", goaloff, "", GetTime(prev_time), GetTime(nx))

                cut = video.subclip(startClip, endClip)

                outputDir = os.path.join(base_dir, "processed", noShotFolder)
                os.makedirs(outputDir, exist_ok=True)
                videoOutFileName = os.path.join(outputDir, f"{match['Id']}_{h+1}_{n}.full.mp4")
                if not os.path.isfile(videoOutFileName) or overwriteExist == True: # skip if exists
                    cut.write_videofile(videoOutFileName)

                if resizeClips:
                    resizedOutputFile = videoOutFileName.replace(".full.mp4",".resized.mp4")
                    resized = cut.fx(vfx.resize, width=video_size[0])
                    if not os.path.isfile(resizedOutputFile) or overwriteExist == True: # skip if exists
                        resized.write_videofile(resizedOutputFile)
            else:
                print("WARNING: Time between shots too small, cannot cut no-shot\n")

        n+=1
        prev_time = x

    # Dump lists
    if len(shotClips)>0:
        with open(os.path.join(base_dir, "processed", shotFolder, 'shots.txt'), 'w') as outfile:
            json.dump(shotClips, outfile, indent=4)

    if len(noshotClips)>0:
        with open(os.path.join(base_dir, "processed", noShotFolder, 'noshots.txt'), 'w') as outfile:
            json.dump(noshotClips, outfile, indent=4)

    # Skips Attacks Shots in DEBUG
    if shotsOnly == True:
        video.reader.close()
        video.audio.reader.close_proc()
        continue

    print("Attacks\n====================================================")

    n=0
    cutList = list(cuts)
    for attackInfo in attacks:

        if args.debugMaxClips != None and n >= args.debugMaxClips:
            print(f"Debug mode - Reached max number of clips ({args.debugMaxClips})\n")
            break

        x = attackInfo['start']

        doClip = True
        for cutInfo in cutList:
            attackClipOverlapStart = cutInfo['start'] - clipTrimAttackTimeRange
            attackClipOverlapEnd = cutInfo['start'] + clipTrimAttackTimeRange
            if x >= attackClipOverlapStart and x <= attackClipOverlapEnd:
                print(f'attackEvent: {x} shotEvent: {cutInfo["start"]} shotClipStart: {attackClipOverlapStart} shotClipEnd: {attackClipOverlapEnd}')
                print(f"WARNING: Attack {attackInfo['MarkGuid']} overlaps with Shot {cutInfo['MarkGuid']}")
                doClip = False
                break

        if doClip:

            startClip = x + dt - mt - satoff - clipTrimAttackStart
            endClip =  x + dt - mt - satoff + clipTrimAttackEnd

            # Skip cuts that are in the wrong "half" of the video
            if eventMatchPhase != h:
                continue

            if eventMatchPhase==0 and h == eventMatchPhase and startClip > video.duration:
                print(f"Cut in the wrong 'half' of the video (h={h+1}). Skipping\n")
                break
            if eventMatchPhase==1 and h == eventMatchPhase and startClip < 0:
                #print(f"Cut in the wrong 'half' of the video (h={h+1}). Skipping\n")
                continue

            if debug:
                print(f"n:               {n}\n" +
                    f"Half:            {h+1}\n" +
                    f"EventMatchPhase: {eventMatchPhase+1}\n" +
                    f"MatchTime:       {GetTime(x)}\n" +
                    f"PhaseStartTime:  {GetTime(mt)}\n" +
                    f"Kickoff:         {GetTime(dt)}\n" +
                    f"SatOff:          {GetTime(satoff)}\n" +
                    f"ClipStart:       {GetTime(startClip)}\n" +
                    f"ClipEnd:         {GetTime(endClip)}\n")

            startClip = x + dt - mt - satoff - clipTrimAttackStart
            endClip =  x + dt - mt - satoff + clipTrimAttackEnd

            # Skip cuts that are in the wrong "half" of the video
            if eventMatchPhase != h:
                continue

            if eventMatchPhase==0 and h == eventMatchPhase and startClip > video.duration:
                print(f"Cut in the wrong 'half' of the video (h={h+1}). Skipping\n")
                break
            if eventMatchPhase==1 and h == eventMatchPhase and startClip < 0:
                #print(f"Cut in the wrong 'half' of the video (h={h+1}). Skipping\n")
                continue

            if debug:
                print(f"n:               {n}\n" +
                    f"Half:            {h+1}\n" +
                    f"EventMatchPhase: {eventMatchPhase+1}\n" +
                    f"MatchTime:       {GetTime(x)}\n" +
                    f"PhaseStartTime:  {GetTime(mt)}\n" +
                    f"Kickoff:         {GetTime(dt)}\n" +
                    f"SatOff:          {GetTime(satoff)}\n" +
                    f"ClipStart:       {GetTime(startClip)}\n" +
                    f"ClipEnd:         {GetTime(endClip)}\n")

            if startClip < 0:
                print(f"StartClip < 0. Forcing to 0.0\n")
                startClip = 0

            if endClip > video.duration:
                print(f"EndClip > video.duration. Forcing to video.duration\n")
                endClip = video.duration

            addClipToList(attackClips, match["Id"], n, h+1, eventMatchPhase+1, GetTime(x), GetTime(mt), GetTime(dt), GetTime(satoff), GetTime(startClip), GetTime(endClip), "false", goaloff, "", "", "")

            cut = video.subclip(startClip, endClip)
            outputDir = os.path.join(base_dir, "processed", attacksFolder)
            os.makedirs(outputDir, exist_ok=True)
            videoOutFileName = os.path.join(outputDir, f"{match['Id']}_{h+1}_{n}.full.mp4")
            if not os.path.isfile(videoOutFileName) or overwriteExist == True: # skip if exists
                cut.write_videofile(videoOutFileName)

            if resizeClips:
                resizedOutputFile = videoOutFileName.replace(".full.mp4",".resized.mp4")
                resized = cut.fx(vfx.resize, width=video_size[0])
                if not os.path.isfile(resizedOutputFile) or overwriteExist == True: # skip if exists
                    resized.write_videofile(resizedOutputFile)
        n+=1

    video.reader.close()
    video.audio.reader.close_proc()

    if len(attackClips)>0:
        with open(os.path.join(base_dir, "processed", attacksFolder, 'attacks.txt'), 'w') as outfile:
            json.dump(attackClips, outfile, indent=4)


text:
# -*- coding: utf-8 -*-
"""Atkinson and Boore (2006, :cite:`atkinson06`) model."""

import numpy as np

from . import model

__author__ = 'Albert Kottke'

class AtkinsonBoore2006(model.GroundMotionModel):
    """Atkinson and Boore (2006, :cite:`atkinson06`) model.

    Developed for the Eastern North America with a reference velocity of 760
    or 2000 m/s.

    Parameters
    ----------
    scenario : :class:`pygmm.model.Scenario`
        earthquake scenario

    """

    NAME = 'Atkinson and Boore (2006)'
    ABBREV = 'AB06'

    # Load the coefficients for the model
    COEFF = dict(
        bc=model.load_data_file('atkinson_boore_2006-bc.csv', 2),
        rock=model.load_data_file('atkinson_boore_2006-rock.csv', 2), )

    PERIODS = COEFF['bc']['period']

    COEFF_SITE = model.load_data_file('atkinson_boore_2006-site.csv', 2)
    COEFF_SF = model.load_data_file('atkinson_boore_2006-sf.csv', 2)

    INDEX_PGD = 0
    INDEX_PGV = 1
    INDEX_PGA = 2
    INDICES_PSA = np.arange(3, 27)

    PARAMS = [
        model.NumericParameter('mag', True),
        model.NumericParameter('dist_rup', True),
        model.NumericParameter('v_s30', True)
    ]

    def __init__(self, scenario: model.Scenario):
        """Initialize the model."""
        super(AtkinsonBoore2006, self).__init__(scenario)
        self._ln_resp = self._calc_ln_resp()
        self._ln_std = self._calc_ln_std()

    def _calc_ln_resp(self) -> np.ndarray:
        """Calculate the natural logarithm of the response.

        Returns
        -------
        ln_resp : class:`np.array`:
            natural log of the response
        """
        s = self._scenario
        c = self.COEFF['bc'] if s.v_s30 else self.COEFF['rock']

        # Compute the response at the reference condition
        r0 = 10.0
        r1 = 70.0
        r2 = 140.0

        f0 = np.maximum(np.log10(r0 / s.dist_rup), 0)
        f1 = np.minimum(np.log10(s.dist_rup), np.log10(r1))
        f2 = np.maximum(np.log10(s.dist_rup / r2), 0)

        # Compute the log10 PSA in units of cm/sec/sec
        log10_resp = (c.c_1 + c.c_2 * s.mag + c.c_3 * s.mag ** 2 +
                      (c.c_4 + c.c_5 * s.mag) * f1 +
                      (c.c_6 + c.c_7 * s.mag) * f2 +
                      (c.c_8 + c.c_9 * s.mag) * f0 + c.c_10 * s.dist_rup)

        # Apply stress drop correction
        log10_resp += self._calc_stress_factor()

        if s.v_s30:
            # Compute the site amplification
            pga_bc = (10 ** log10_resp[self.INDEX_PGA])

            log10_site = self._calc_log10_site(pga_bc)

            log10_resp += log10_site

        # Convert from cm/sec/sec to gravity
        log10_resp -= np.log10(980.665)

        ln_resp = np.log(10 ** log10_resp)
        return ln_resp

    def _calc_ln_std(self) -> np.ndarray:
        """Calculate the logarithmic standard deviation.

        Returns
        -------
        ln_std : class:`np.array`:
            natural log standard deviation

        """
        ln_std = np.ones_like(self.PERIODS) * 0.30
        return ln_std

    def _calc_stress_factor(self) -> float:
        """Calculate the stress correction factor proposed by Atkinson and
        Boore (2011) :cite:`atkinson11`.

        Returns
        -------
        log10_stress_factor : class:`np.array`:
            log base 10 of the stress factor
        """
        s = self._scenario
        c = self.COEFF_SF

        stress_drop = 10. ** (3.45 - 0.2 * s.mag)
        v1 = c.delta + 0.05
        v2 = (0.05 + c.delta * np.maximum(s.mag - c.m_1, 0) / (c.m_h - c.m_1))

        log10_stress_factor = (np.minimum(2., stress_drop / 140.) *
                               np.minimum(v1, v2))

        return np.interp(self.PERIODS, c.period, log10_stress_factor)

    def _calc_log10_site(self, pga_bc: float) -> np.ndarray:
        """Calculate the log10 of the site amplification.

        Parameters
        ----------
        pga_bc : float
            peak ground acceleration (PGA, g) at the B/C boundary.

        Returns
        -------
        log_10_site : :class:`np.ndarray`
            log base 10 of the  site amplification.
        """
        s = self._scenario
        c = self.COEFF_SITE
        VS_1 = 180.
        VS_2 = 300.
        VS_REF = 760.

        if s.v_s30 <= VS_1:
            b_nl = c.b_1
        elif VS_1 < s.v_s30 <= VS_2:
            b_nl = (
                (c.b_1 - c.b_2) * np.log(s.v_s30 / VS_2) / np.log(VS_1 / VS_2))
        elif VS_2 < s.v_s30 <= VS_REF:
            b_nl = (c.b_2 * np.log(s.v_s30 / VS_REF) / np.log(VS_2 / VS_REF))
        else:
            # Vs30 > VS_REF
            b_nl = 0

        pga_bc = max(pga_bc, 60.)

        log10_site = np.log10(
            np.exp(c.b_lin * np.log(s.v_s30 / VS_REF) + b_nl * np.log(pga_bc /
                                                                      100.)))
        return np.interp(self.PERIODS, c.period, log10_site)

max_stars_repo_name: glavrentiadis/pygmm

max_issues_repo_name: glavrentiadis/pygmm

max_forks_repo_name: glavrentiadis/pygmm


text:
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 28 16:53:42 2019

@author: WT
"""

import os
import librosa
import soundfile as sf
import numpy as np
import re
from torch.utils.data import Dataset, DataLoader
import pandas as pd
from .utils import save_as_pickle, load_pickle, tokener, vocab
from sklearn.preprocessing import MinMaxScaler, StandardScaler
from tqdm import tqdm
import logging

logging.basicConfig(format='%(asctime)s [%(levelname)s]: %(message)s', \
                    datefmt='%m/%d/%Y %I:%M:%S %p', level=logging.INFO)
logger = logging.getLogger(__file__)

class args():
    def __init__(self):
        self.n_mfcc = 13
        self.n_fft = 2048
        self.hop_length = 512
        self.batch_size = 5
        self.use_lg_mels = 1
        self.n_mels = 80
        self.max_frame_len = 3171

def separate_info(x):
    id_ = re.search("[\d+-]+", x)[0]
    speaker = re.search("[\d+]+-", x)[0][:-1]
    text = re.search("([A-Z]+.+)+", x)[0]
    return id_, speaker, text

def pad_sos_eos(x):
    return [0] + x + [2]

def normal_scaler(x, mu, std):
    return (x-mu)/std

def min_max_scaler(x, min_, max_, feature_range=[0.3,2.3]):
    x_std = (x - min_)/(max_- min_)
    x_scaled = x_std*(feature_range[1] - feature_range[0]) + feature_range[0]
    return x_scaled

def stack_downsampler(feature, args, stack=3):
    '''stack: number of frames to stack
    shape = 3 X n_features X length '''
    no_stacks = feature.shape[-1]//stack
    feature = feature.reshape(3, -1, feature.shape[2]//3)
    stacked_feature = np.zeros((feature.shape[0], stack*feature.shape[1], feature.shape[-1]//stack))
    for i in range(no_stacks):
        pass
    return feature

def extract_feature(filepath, args):
    samples, sample_rate = sf.read(filepath)
    if args.use_lg_mels == 1:
        S = librosa.feature.melspectrogram(samples.astype(float), sr=sample_rate, n_mels=args.n_mels,\
                                           n_fft=int(sample_rate*(args.n_fft/1000)), \
                                           hop_length=int(sample_rate*(args.hop_length/1000)),\
                                           power=1.0)
        original_len = S.shape[-1]; #print(S.shape)
        if (args.max_frame_len - original_len) >= 0:
            S = np.append(S, (1e-9)*np.ones((args.n_mels, args.max_frame_len-original_len)), axis=1)
            S = librosa.amplitude_to_db(S) # convert to log mel energies
            delta1_S = librosa.feature.delta(S, order=1)
            delta2_S = librosa.feature.delta(S, order=2)
            #S[:,:original_len] = scaler.fit_transform(S[:,:original_len].T).T
            #delta1_S[:,:original_len] = scaler.fit_transform(delta1_S[:,:original_len].T).T
            #delta2_S[:,:original_len] = scaler.fit_transform(delta2_S[:,:original_len].T).T
            features = np.stack((S, delta1_S, delta2_S), axis=0); #break; break;break
            #features = librosa.power_to_db(S, ref=np.max)
            #features = scaler.fit_transform(S.T).T
        else:
            features = None
    else:
        mfcc = librosa.feature.mfcc(y=samples.astype(float), sr=sample_rate, n_mfcc=args.n_mfcc,\
                                                        n_fft=int(sample_rate*(args.n_fft/1000)), \
                                                        hop_length=int(sample_rate*(args.hop_length/1000)))
        original_len = mfcc.shape[-1]
        if (args.max_frame_len - original_len) >= 0:
            mfcc = np.append(mfcc, np.zeros((args.n_mfcc, args.max_frame_len-original_len)), axis=1)
            delta1_mfcc = librosa.feature.delta(mfcc, order=1)
            delta2_mfcc = librosa.feature.delta(mfcc, order=2)
            #mfcc[:,:original_len] = scaler.fit_transform(mfcc[:,:original_len].T).T
            #delta1_mfcc[:,:original_len] = scaler.fit_transform(delta1_mfcc[:,:original_len].T).T
            #delta2_mfcc[:,:original_len] = scaler.fit_transform(delta2_mfcc[:,:original_len].T).T
            features = np.stack((mfcc, delta1_mfcc, delta2_mfcc), axis=0)
        else:
            features = None
    return features, original_len

def get_mfcc_data(args):
    ''' Extracts MFCC 0th, 1st, 2nd order coefficients, tokenizes text transcript, build vocab, convert text tokens to ids
    and saves results in pickle file
    use_lg_mels: if True, use log of Mel spectrogram instead of MFCC coefficients as features
    '''
    logging.info("Extracting MFCC/Mel features from data...")
    folder = args.folder
    df = pd.DataFrame(columns=['id', 'speaker', 'text','features', 'features_len'])
    tokenizer_en = tokener()
    #scaler = MinMaxScaler(feature_range=(0.0001, 2), copy=True)
    #scaler = StandardScaler()
    def dict_data(data, id_, features=True):
        if id_ in data.keys():
            if features:
                return data[id_][0]
            else:
                return data[id_][1]

    for speaker in tqdm(os.listdir("./data/%s/" % folder)):
        for chapter in os.listdir("./data/%s/%s/" % (folder, speaker)):
            data = {}
            for file in os.listdir("./data/%s/%s/%s/" % (folder, speaker, chapter)):
                if '.txt' in file:
                    s_path = "./data/%s/%s/%s/%s" % (folder, speaker, chapter, file)
                    df_dum = pd.read_csv(s_path, names=['id', 'text'])
                    df_dum.loc[:, 'text'] = df_dum.apply(lambda x: separate_info(x['id'])[2], \
                                                          axis=1).apply(lambda x: tokenizer_en.tokenize(x))
                    df_dum.loc[:, 'speaker'] = df_dum.apply(lambda x: separate_info(x['id'])[1], axis=1)
                    df_dum.loc[:, 'id'] = df_dum.apply(lambda x: separate_info(x['id'])[0], axis=1)

                elif '.flac' in file:
                    s_path = "./data/%s/%s/%s/%s" % (folder, speaker, chapter, file)
                    features, original_len = extract_feature(s_path, args)
                    if features is None:
                        continue
                    id_ = re.search("[\d+-]+", file)[0]
                    if (id_ is not None):
                        data[id_] = features, original_len

            df_dum["features"] = df_dum.apply(lambda x: dict_data(data, x['id'], features=True), axis=1)
            df_dum["features_len"] = df_dum.apply(lambda x: dict_data(data, x['id'], features=False), axis=1)
            df_dum.dropna(inplace=True)
            df = df.append(df_dum, ignore_index=True)

    logger.info("Building vocab and converting transcript to tokens...")
    v = vocab(level=args.level)
    v.build_vocab(df["text"])
    df.loc[:, "text"] = df.apply(lambda x: v.convert_w2idx(x["text"]), axis=1)
    df.loc[:, "text"] = df.apply(lambda x: pad_sos_eos(x["text"]), axis=1)

    logger.info("Speaker normalization...")
    logger.info("Getting speakers features mean and std...")
    speaker_stats = {}
    for speaker in tqdm(df["speaker"].unique()):
        df_dum = df[df["speaker"] == speaker]
        channel_mu = []; channel_std = []
        for channel in range(3):
            f_list = []
            for row, l in zip(df_dum["features"], df_dum["features_len"]):
                row = row[:,:,:int(l)]
                f_list.extend(list(row[channel].reshape(-1)))
            f_list = np.array(f_list)
            channel_mu.append(f_list.mean())
            channel_std.append(f_list.std())
        speaker_stats[speaker] = (channel_mu, channel_std)
    save_as_pickle("speaker_stats.pkl", speaker_stats)

    logger.info("Normalizing...")
    def speaker_norm(feature, stats):
        channel_mu, channel_std = stats
        for idx, (mu, std) in enumerate(zip(channel_mu, channel_std)):
            feature[idx, :, :] = (feature[idx, :, :] - mu)/std
        return feature

    df["features"] = df.apply(lambda x: speaker_norm(x["features"], speaker_stats[x["speaker"]]), axis=1)

    logging.info("Saving...")
    save_as_pickle("df_%s.pkl" % folder, df)
    save_as_pickle("vocab.pkl", v)
    logging.info("Saved!")

class padded_dataset(Dataset):
    def __init__(self, df, args, labels=True):
        def x_padder(x, max_len):
            if x.shape[-1] < max_len:
                if args.use_lg_mels == 1:
                    x = np.append(x, np.zeros((3, args.n_mels, max_len-x.shape[-1])), axis=2)
                else:
                    x = np.append(x, np.zeros((3, args.n_mfcc, max_len-x.shape[-1])), axis=2)
            return x

        def y_padder(x, max_len):
            x = np.array(x, dtype=int)
            if len(x) < max_len:
                x = np.append(x, np.ones((max_len-len(x)), dtype=int)) # 1 is the idx for <pad> token
            return x

        self.labels = labels
        self.x_max_len = int(max(df['features_len']))
        self.X = df["features"]
        self.X_len = df["features_len"]
        if labels == True:
            self.y_max_len = max(df['text'].apply(lambda x: len(x)))
            self.y = df["text"].apply(lambda x: y_padder(x, self.y_max_len))
        #self.X = df["features"].apply(lambda x: x_padder(x, self.x_max_len))
        else:
            self.y = y_padder(np.zeros(1, dtype=int), args.max_seq_len)

    def __len__(self):
        return len(self.X)

    def __getitem__(self, idx):
        if self.labels == True:
            # output X shape = 3 X n_features X length
            return self.X.iloc[idx], self.y.iloc[idx], self.X_len.iloc[idx]
        else:
            return self.X.iloc[idx], self.y, self.X_len.iloc[idx] # return <sos> if labels=False

def load_dataloaders(args):
    ''' loads mfcc data, convert words into ids and returns corresponding data_loaders '''
    logger.info("Loading dataloaders...")
    folder = args.folder
    train_path = "./data/" + "df_%s.pkl" % folder
    if not os.path.isfile(train_path):
        get_mfcc_data(args)
    df = load_pickle("df_%s.pkl" % folder)

    trainset = padded_dataset(df, args)
    train_loader = DataLoader(trainset, batch_size=args.batch_size, shuffle=True, \
                              num_workers=0, pin_memory=False)
    train_length = len(trainset)
    max_features_length = trainset.x_max_len
    max_seq_len = trainset.y_max_len
    return train_loader, train_length, max_features_length, max_seq_len


text:
#!/usr/bin/env python

''' ---------------- About the script ----------------

Assignment 4: Classification benchmarks: Logistic Regression

This script takes a MNIST digits dataset, trains a Logistic Regression Classifier, prints the evaluation metrics to the terminal and saves classification report and a confusion matrix to a directory. Optionally, user can define both training data and test data sizes, name of a classification report file. There are no positional arguments, only default values of optional arguments.

Optional argument:
    -trs, --train_size:       The size of the training data as a percentage, where the default = 0.8 (80%)
    -tes, --test_size:        The size of the test data as a percentage, where the default = 0.2 (20%)
    -n,   --name:             Name of the classification report to be saved as .csv file, where the default = logReg_report

Example:
    customized optional arguments:
        $ python Logistic_Regression.py -trs 0.9 -tes 0.1 -n lr_cm.csv

    default arguments:
        $ python Logistic_Regression.py

'''

"""---------------- Importing libraries ----------------
"""

# System tools
import os
import sys
sys.path.append(os.path.join(".."))

# Import teaching utils
import numpy as np
import utils.classifier_utils as clf_util

# Import pandas for working with dataframes
import pandas as pd

# visualization tool
import matplotlib.pyplot as plt

# Import sklearn metrics
from sklearn import metrics
from sklearn.datasets import fetch_openml
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score

# Command-line interface
import argparse

"""---------------- Main script ----------------
"""

def main():

    """------ Argparse parameters ------
    """
    # Instantiating the ArgumentParser  object as parser
    parser = argparse.ArgumentParser(description = "[INFO] Classify MNIST data and print out performance accuracy report")

    # Adding optional arguments
    parser.add_argument("-trs", "--train_size", required = False, default = 0.8, type = float, help = "The size of the training data as a percentage. Default = 0.8 (80%)")
    parser.add_argument("-tes", "--test_size", required = False, default = 0.2, type = float, help = "The size of the testing data as a percentage. Default = 0.2 (20%)")
    parser.add_argument("-n", "--name", required=False, default = "logReg_report", help="Name of the classification report to be saved as .csv file. Default = logReg_report")

    # Parsing the arguments
    args = vars(parser.parse_args())

    """------ Loading full data and preprocessing ------
    """

    print("[INFO] Loading and preprocessing data...")

    # Fetching data
    X, y = fetch_openml("mnist_784", version=1, return_X_y=True)

    # Converting to numpy arrays
    X = np.array(X) # Data
    y = np.array(y) # Classes (0-9)

    # Creating training data and test dataset
    X_train, X_test, y_train, y_test = train_test_split(X,
                                                        y,
                                                        train_size = args["train_size"],
                                                        test_size= args["test_size"])

    # Scaling the features from 0-255 to between 0 and 1
    X_train_scaled = X_train/255.0
    X_test_scaled = X_test/255.0

    """------ Loading sample data and preprocessing ------
    """

    # Load sample data
    ##digits = datasets.load_digits()
    # Convert to floats
    ##data = digits.data.astype("float")
    # MinMax regularization
    ##data = (data - data.min())/(data.max() - data.min())
    # Creating train and test datasets
    ##X_train, X_test, y_train, y_test = train_test_split(data,
                                                  ##digits.target,
                                                  ##train_size = args["train_size"],
                                                  ##test_size= args["test_size"])

    """------ Training the model ------
    """

    # Defining a model
    print("[INFO] training Logistic Regression model...")
    clf = LogisticRegression(penalty='none',
                         tol=0.1,
                         solver='saga',
                         multi_class='multinomial').fit(X_train_scaled, y_train) #Fitting classifier to test data

    """------ Evaluating the model: accuracy measures and classifier report ------
    """

    print("[INFO] evaluating Logistic Regression model...")
    # Predicted labels
    y_pred = clf.predict(X_test_scaled)
    # Calculating  accuracy of the model based on comparing the predicted labels with the actual labels
    accuracy = accuracy_score(y_test, y_pred)
    # Logistic Regression classifier report
    cm = metrics.classification_report(y_test, y_pred)

    """------ Saving classification report as .csv file ------
    """

    # Create ouput folder for saving the classification report if it doesn´t exist already
    if not os.path.exists("../out"):
        os.makedirs("../out")

    # Turning classification report into a dataframe
    report_df = pd.DataFrame(metrics.classification_report(y_test, y_pred, output_dict = True)).transpose()
    # Defining full filepath to save csv file
    outfile = os.path.join("..", "out", args["name"])
    # Saving a dataframe as .csv
    report_df.to_csv(outfile)
    # Printing that .csv file has been saved
    print(f"\n[INFO] classification report is saved in directory {outfile}")

    """------ Saving confusion matrix as .png file ------
    """

    # Defining full filepath to save .png file
    path_png = os.path.join("..", "out", "logReg_confusion_matrix.png")
    # Creating confusion matrix
    clf_util.plot_cm(y_test, y_pred, normalized=True)
    # Saving as .png file
    plt.savefig(path_png)
    # Printing that .png file has been saved
    print(f"\n[INFO] confusion matrix is saved in directory {path_png}")

    """------ Final messages and printing results ------
    """

    print("\nSUCCESS! Logistic regression model has been trained and evaluated. Full evaluation of the model is displayed below")
    print(f"\nWeighted average accuracy of the current logistic regression model is = {accuracy}")
    print(cm) # full classification report

# Define behaviour when called from command line
if __name__=="__main__":
    main()


text:
import logging
import os

import astropy.io.fits
import galsim

import test_case_definitions

_logger = logging.getLogger(__name__)

RANDOM_SEED = 1234

FITS_COMMENT_PREPEND = 'Celeste: '

def add_header_to_hdu(hdu, header_dict):
    header = galsim.fits.FitsHeader(hdu.header)
    for name, (value, comment) in header_dict.iteritems():
        header[name] = (value, FITS_COMMENT_PREPEND + comment)

def save_multi_extension_fits(hdu_list, filename):
    if not os.path.exists(os.path.dirname(filename)):
        os.mkdir(os.path.dirname(filename))
    if os.path.exists(filename):
        os.remove(filename)
    hdu_list.writeto(filename)
    _logger.info('Wrote multi-extension FITS image to %r', filename)

def main():
    logging.basicConfig(format="%(message)s", level=logging.INFO)

    _logger.info('Generating %d test cases', len(test_case_definitions.TEST_CASE_FNS))
    fits_hdus = astropy.io.fits.HDUList()
    for case_index, test_case_fn in enumerate(test_case_definitions.TEST_CASE_FNS):
        _logger.info('  Generating case %s', test_case_fn.__name__)
        uniform_deviate = galsim.UniformDeviate(RANDOM_SEED + case_index)
        test_case = test_case_definitions.GalSimTestCase()
        test_case_fn(test_case)
        for band_index in xrange(5):
            image = test_case.construct_image(band_index, uniform_deviate)
            galsim.fits.write(image, hdu_list=fits_hdus)
            add_header_to_hdu(fits_hdus[-1], test_case.get_fits_header(case_index, band_index))

    image_file_name = os.path.join('output', 'galsim_test_images.fits')
    save_multi_extension_fits(fits_hdus, image_file_name)

if __name__ == "__main__":
    main()

reason: import astropy


text:
from __future__ import print_function, division
import argparse
import os
import torch
import torch.nn as nn
import torch.nn.parallel
import torch.backends.cudnn as cudnn
import torch.optim as optim
import torch.utils.data
from torch.autograd import Variable
import torchvision.utils as vutils
import torch.nn.functional as F
import numpy as np
import time
from tensorboardX import SummaryWriter
from datasets import __datasets__
from models import __models__, model_loss
from utils import *
from torch.utils.data import DataLoader
from math import isnan
import torchvision
import gc
import cv2 as cv
from os.path import join, split, isdir, isfile, splitext, split, abspath, dirname

import matplotlib.pyplot as plt
from pylab import *
import matplotlib
matplotlib.use('pdf')

cudnn.benchmark = True

parser = argparse.ArgumentParser(description='Group-wise Correlation Stereo Network (GwcNet)')
parser.add_argument('--model', default='gwcnet-gc', help='select a model structure', choices=__models__.keys())
parser.add_argument('--maxdisp', type=int, default=384, help='maximum disparity')

parser.add_argument('--dataset', required=True, help='dataset name', choices=__datasets__.keys())
parser.add_argument('--datapath', required=True, help='data path')
parser.add_argument('--trainlist', required=True, help='training list')
parser.add_argument('--testlist', required=True, help='testing list')

parser.add_argument('--lr', type=float, default=0.001, help='base learning rate')
parser.add_argument('--batch_size', type=int, default=4, help='training batch size')
parser.add_argument('--test_batch_size', type=int, default=1, help='testing batch size')
parser.add_argument('--epochs', type=int, required=True, help='number of epochs to train')
parser.add_argument('--lrepochs', type=str, required=True, help='the epochs to decay lr: the downscale rate')

parser.add_argument('--logdir', required=True, help='the directory to save logs and checkpoints')
parser.add_argument('--loadckpt', help='load the weights from a specific checkpoint')
#parser.add_argument('--resume', action='store_true', help='continue training the model')
parser.add_argument('--resume', default=False, help='continue training the model')
parser.add_argument('--seed', type=int, default=1, metavar='S', help='random seed (default: 1)')

parser.add_argument('--summary_freq', type=int, default=20, help='the frequency of saving summary')
parser.add_argument('--save_freq', type=int, default=1, help='the frequency of saving checkpoint')

# parse arguments, set seeds
args = parser.parse_args()
torch.manual_seed(args.seed)
torch.cuda.manual_seed(args.seed)
os.makedirs(args.logdir, exist_ok=True)

# test_gt_path='/data/yyx/GwcNet-master/checkpoints/kitti/ft_from0/kitti_test_gt/'
# test_pred1_path='/data/yyx/GwcNet-master/checkpoints/kitti/ft_from0/kitti_test_pred/'
# os.makedirs(test_gt_path exist_ok=True)
# os.makedirs(test_pred1_path, exist_ok=True)

# create summary logger
print("creating new summary file")
logger = SummaryWriter(args.logdir)

# dataset, dataloader
StereoDataset = __datasets__[args.dataset]
train_dataset = StereoDataset(args.datapath, args.trainlist, True)
test_dataset = StereoDataset(args.datapath, args.testlist, False)
TrainImgLoader = DataLoader(train_dataset, args.batch_size, shuffle=True, num_workers=8, drop_last=True)
TestImgLoader = DataLoader(test_dataset, args.test_batch_size, shuffle=False, num_workers=4, drop_last=False)

# model, optimizer
model = __models__[args.model](args.maxdisp)
model = nn.DataParallel(model)
model.cuda()
optimizer = optim.Adam(model.parameters(), lr=args.lr, betas=(0.9, 0.999))

# for index,(name,value) in enumerate(model.named_parameters()):
#     #value.requires_grad = (index < last)
#     #value.requires_grad = False
#     print(index+1, name," : ",value.requires_grad)

# load parameters
start_epoch = 0
if args.resume:
    # find all checkpoints file and sort according to epoch id
    all_saved_ckpts = [fn for fn in os.listdir(args.logdir) if fn.endswith(".ckpt")]
    all_saved_ckpts = sorted(all_saved_ckpts, key=lambda x: int(x.split('_')[-1].split('.')[0]))
    # use the latest checkpoint file
    loadckpt = os.path.join(args.logdir, all_saved_ckpts[-1])
    print("loading the lastest model in logdir: {}".format(loadckpt))
    state_dict = torch.load(loadckpt)
    model.load_state_dict(state_dict['model'])
    optimizer.load_state_dict(state_dict['optimizer'])
    start_epoch = state_dict['epoch'] + 1
elif args.loadckpt:
    # load the checkpoint file specified by args.loadckpt
    print("loading model {}".format(args.loadckpt))
    state_dict = torch.load(args.loadckpt)
    model.load_state_dict(state_dict['model'])
    start_epoch = state_dict['epoch'] + 1
print("start at epoch {}".format(start_epoch))

def train():
    for epoch_idx in range(start_epoch, args.epochs):
        adjust_learning_rate(optimizer, epoch_idx, args.lr, args.lrepochs)
        print("current rate is ",args.lr)
        print(args.maxdisp)
        #avg_train_scalars = AverageMeter()
        all_loss=0.0
        #print(all_loss)
        # training
        # nannum=0
        # for batch_idx, sample in enumerate(TrainImgLoader):
        #     #print(len(sample))
        #     #print("!!!!!!!!!!!!!")
        #     global_step = len(TrainImgLoader) * epoch_idx + batch_idx
        #     start_time = time.time()
        #     do_summary = global_step % args.summary_freq == 0
        #     loss, scalar_outputs, image_outputs = train_sample(sample,compute_metrics=do_summary)
        #     #loss_all=loss_all+loss
        #     if math.isnan(loss):
        #         print(batch_idx,'loss is nan')
        #         nannum=nannum+1
        #     else:
        #         all_loss=loss+all_loss
        #     if do_summary:
        #         save_scalars(logger, 'train', scalar_outputs, global_step)
        #         save_images(logger, 'train', image_outputs, global_step)

        #     #avg_train_scalars.update(loss)
        #     #print("!!!",avg_train_scalars)
        #     del scalar_outputs, image_outputs
        #     if (batch_idx%100==0)or(batch_idx==len(TrainImgLoader)-1):
        #         print('Epoch {}/{}, Iter {}/{}, train loss = {:.3f}, time = {:.3f},losses averange={:.3f}'.format(epoch_idx, args.epochs,
        #                                                                                batch_idx,
        #                                                                                len(TrainImgLoader), loss,
        #                                                                                time.time() - start_time,all_loss/(batch_idx-nannum+1)))

        # #loss_all=loss_all.avg()
        # #avg_train_scalars=avg_train_scalars.mean()
        # #print("avg_train_loss",all_loss/(len(TestImgLoader)-nannum))
        #         #print("losses",loss_all/(batch_idx+1))
        # #print("loss averange=",loss_all/12537.0)
        # #saving checkpoints
        # if (epoch_idx + 1) % args.save_freq == 0:
        #     checkpoint_data = {'epoch': epoch_idx, 'model': model.state_dict(), 'optimizer': optimizer.state_dict()}
        #     torch.save(checkpoint_data, "{}/checkpoint_{:0>6}.ckpt".format(args.logdir, epoch_idx))
        # gc.collect()
        # print("avg_train_loss",all_loss/(len(TrainImgLoader)-nannum))

        # testing
        avg_test_scalars = AverageMeterDict()
        for batch_idx, sample in enumerate(TestImgLoader):
            #print("sample",sample)
            global_step = len(TestImgLoader) * epoch_idx + batch_idx
            start_time = time.time()
            do_summary = global_step % args.summary_freq == 0
            loss, scalar_outputs, image_outputs = test_sample(batch_idx,sample, compute_metrics=do_summary)
            if do_summary:
                save_scalars(logger, 'test', scalar_outputs, global_step)
                save_images(logger, 'test', image_outputs, global_step)
            avg_test_scalars.update(scalar_outputs)
            del scalar_outputs, image_outputs
            if (batch_idx%100==0)or(batch_idx==len(TestImgLoader)-1):
                print('Epoch {}/{}, Iter {}/{}, test loss = {:.3f}, time = {:3f}'.format(epoch_idx, args.epochs,
                                                                                batch_idx,
                                                                                len(TestImgLoader), loss,
                                                                                time.time() - start_time))
            #print(avg_test_scalars["EPE"])
        avg_test_scalars = avg_test_scalars.mean()
        save_scalars(logger, 'fulltest', avg_test_scalars, len(TrainImgLoader) * (epoch_idx + 1))
        print("avg_test_scalars", avg_test_scalars)
        gc.collect()

# train one sample
def train_sample(sample, compute_metrics=False):
    model.train()

    imgL, imgR, disp_gt = sample['left'], sample['right'], sample['disparity']
    imgL = imgL.cuda()
    imgR = imgR.cuda()
    disp_gt = disp_gt.cuda()

    optimizer.zero_grad()

    disp_ests = model(imgL, imgR)

    #print(len(disp_ests))
    #print(disp_ests[0].shape)
    mask = (disp_gt < args.maxdisp) & (disp_gt > 0)
    #print("##########",mask.shape)
    loss = model_loss(disp_ests, disp_gt, mask)
    #loss = model_loss(disp_ests, disp_gt, mask,imgL, imgR)

    scalar_outputs = {"loss": loss}
    image_outputs = {"disp_est": disp_ests, "disp_gt": disp_gt, "imgL": imgL, "imgR": imgR}
    if compute_metrics:
        with torch.no_grad():
            image_outputs["errormap"] = [disp_error_image_func()(disp_est, disp_gt) for disp_est in disp_ests]
            scalar_outputs["EPE"] = [EPE_metric(disp_est, disp_gt, mask) for disp_est in disp_ests]
            scalar_outputs["D1"] = [D1_metric(disp_est, disp_gt, mask) for disp_est in disp_ests]
            scalar_outputs["Thres1"] = [Thres_metric(disp_est, disp_gt, mask, 1.0) for disp_est in disp_ests]
            scalar_outputs["Thres2"] = [Thres_metric(disp_est, disp_gt, mask, 2.0) for disp_est in disp_ests]
            scalar_outputs["Thres3"] = [Thres_metric(disp_est, disp_gt, mask, 3.0) for disp_est in disp_ests]
    loss.backward()
    optimizer.step()

    return tensor2float(loss), tensor2float(scalar_outputs), image_outputs

# test one sample
@make_nograd_func
def test_sample(batch_idx,sample, compute_metrics=True):
    model.eval()
    #model.train()

    imgL, imgR, disp_gt ,mask= sample['left'], sample['right'], sample['disparity'] ,sample["mask"]
    imgL = imgL.cuda()
    imgL = imgL.cuda()
    imgR = imgR.cuda()
    mask=mask.cuda()

    disp_gt = disp_gt.cuda()
    disp_ests = model(imgL, imgR)
    #disp_ests = model(imgL, imgR)
    # print("left_fea",left_fea[0].shape)
    # print("right_fea",right_fea[0].shape)
    # print("cost",cost.shape)
    #mask = (disp_gt < args.maxdisp) & (disp_gt > 0)

    #select_fea_scale=left_fea[0]
    #svisualize_feature_map(select_fea_scale)
    #disp_ests[0][mask==0]=0

    #fea store
    # iml_2=F.upsample(imgL, [select_fea_scale.size()[2], select_fea_scale.size()[3]], mode='bilinear')
    # left_feature_path='/data/yyx/GwcNet-master/checkpoints/sceneflow_monkaa/gwcnet-gc/left_feature/'
    # _, H, W = select_fea_scale[0].shape
    # im_fea = torch.zeros(( H, W))
    # for f in range(select_fea_scale.size()[1]):
    #     im_fea = np.array(select_fea_scale[0][f,:, :].cpu()*255, dtype=np.uint16)
    #     cv.imwrite(join(left_feature_path, "left_feature-%d.png" % (f+1)),im_fea)
    # #im = np.array(iml_2[0,:,:,:].permute(1,2,0).cpu()*255, dtype=np.uint8)
    # im = np.array(iml_2[0,:,:,:].permute(1,2,0).cpu()*255, dtype=np.uint8)
    # cv.imwrite(join(left_feature_path, "itercolor-%d.jpg" % batch_idx),im)

    #driving test
    # test_gt_path='/data/yyx/GwcNet-master/checkpoints/sceneflow_monkaa/gwcnet-gc-gray/sceneflow_test_gt/'
    # test_pred1_path='/data/yyx/GwcNet-master/checkpoints/sceneflow_monkaa/gwcnet-gc-gray/sceneflow_test_pred/'

    # _, H, W = disp_ests[0].shape
    # gt=torch.zeros((H,W))
    # im_pred1 = torch.zeros(( H, W))

    # gt = np.array(disp_gt[0,:, :].cpu()*2, dtype=np.uint16)
    # cv.imwrite(join(test_gt_path, "sceneflow-gt-%d.png" % batch_idx),gt)
    # im_pred1 = np.array(disp_ests[3][0,:, :].cpu()*2, dtype=np.uint16)
    # cv.imwrite(join(test_pred1_path, "sceneflow-%d.png" % batch_idx),im_pred1)

    #kitti test
    # test_gt_path='/data/yyx/GwcNet-master/checkpoints/kitti/pri_192/kitti_test_gt/'
    # test_pred1_path='/data/yyx/GwcNet-master/checkpoints/kitti/pri_192/kitti_test_pred/'

    # # # os.makedirs(test_gt_path exist_ok=True)
    # # # os.makedirs(test_pred1_path, exist_ok=True)
    # _, H, W = disp_ests[0].shape
    # gt=torch.zeros((H,W))
    # im_pred1 = torch.zeros(( H, W))

    # gt = np.array(disp_gt[0,:, :].cpu()*256 +1, dtype=np.uint16)
    # cv.imwrite(join(test_gt_path,  "kitti-gt-%d.png" % batch_idx),gt)
    # im_pred1 = np.array(disp_ests[0][0,:, :].cpu()*256, dtype=np.uint16)
    # cv.imwrite(join(test_pred1_path,  "kitti-pred-%d.png" % batch_idx),im_pred1)

    loss = model_loss(disp_ests, disp_gt,mask)
    #loss = model_loss(disp_ests, disp_gt,mask,imgL, imgR)

    scalar_outputs = {"loss": loss}
    image_outputs = {"disp_est": disp_ests, "disp_gt": disp_gt, "imgL": imgL, "imgR": imgR}
    image_outputs["errormap"] = [disp_error_image_func()(disp_ests[-1], disp_gt)]
    scalar_outputs["EPE"] = [EPE_metric(disp_ests[-1], disp_gt, mask)]
    scalar_outputs["D1"] = [D1_metric(disp_ests[-1], disp_gt, mask)]
    scalar_outputs["D3"] = [D3_metric(disp_ests[-1], disp_gt, mask)]
    scalar_outputs["Thres1"] = [Thres_metric(disp_ests[-1], disp_gt, mask, 1.0)]
    scalar_outputs["Thres2"] = [Thres_metric(disp_ests[-1], disp_gt, mask, 2.0)]
    scalar_outputs["Thres3"] = [Thres_metric(disp_ests[-1], disp_gt, mask, 3.0)]
    # scalar_outputs["D1"] = [D1_metric(disp_est, disp_gt, mask) for disp_est in disp_ests]
    # scalar_outputs["D3"] = [D3_metric(disp_est, disp_gt, mask) for disp_est in disp_ests]
    # scalar_outputs["EPE"] = [EPE_metric(disp_est, disp_gt, mask) for disp_est in disp_ests]
    # scalar_outputs["Thres1"] = [Thres_metric(disp_est, disp_gt, mask, 1.0) for disp_est in disp_ests]
    # scalar_outputs["Thres2"] = [Thres_metric(disp_est, disp_gt, mask, 2.0) for disp_est in disp_ests]
    # scalar_outputs["Thres3"] = [Thres_metric(disp_est, disp_gt, mask, 3.0) for disp_est in disp_ests]

    # scalar_outputs = {"loss": loss}
    # image_outputs = {"disp_est": disp_ests, "disp_gt": disp_gt, "imgL": imgL, "imgR": imgR}

    # scalar_outputs["D1"] = [D1_metric(disp_est, disp_gt, mask) for disp_est in disp_ests]
    # scalar_outputs["D3"] = [D3_metric(disp_est, disp_gt, mask) for disp_est in disp_ests]
    # scalar_outputs["EPE"] = [EPE_metric(disp_est, disp_gt, mask) for disp_est in disp_ests]
    # scalar_outputs["Thres1"] = [Thres_metric(disp_est, disp_gt, mask, 1.0) for disp_est in disp_ests]
    # scalar_outputs["Thres2"] = [Thres_metric(disp_est, disp_gt, mask, 2.0) for disp_est in disp_ests]
    # scalar_outputs["Thres3"] = [Thres_metric(disp_est, disp_gt, mask, 3.0) for disp_est in disp_ests]

    if compute_metrics:
        image_outputs["errormap"] = [disp_error_image_func()(disp_est, disp_gt) for disp_est in disp_ests]

    return tensor2float(loss), tensor2float(scalar_outputs), image_outputs

def get_row_col(num_pic):
    squr = num_pic ** 0.5
    row = round(squr)
    col = row + 1 if squr - row > 0 else row
    return row, col

def visualize_feature_map(img_batch):
    feature_map = img_batch.cpu()
    feature_map=feature_map.data.numpy()
    print("fea-heat",feature_map.shape)

    feature_map_combination = []
    plt.figure()

    num_pic = feature_map.shape[1]
    print("num,pic",num_pic)
    row, col = get_row_col(num_pic)
    print(row,col)

    for i in range(0, num_pic):
        feature_map_split = feature_map[0,i, :, :]
        #print(feature_map_split.shape)
        feature_map_combination.append(feature_map_split)
        #plt.subplot(row, col, i + 1)
        plt.imshow(feature_map_split)
        axis('off')
        filename='fea_map/kitti/0/'+str(i+1)+'.png'
        plt.savefig(filename)
    #plt.show()

    # 各个特征图按1：1 叠加
    feature_map_sum=feature_map_split[0]
    for i in range(len(feature_map_combination)-1):
        feature_map_sum = feature_map_combination[i+1]+feature_map_sum
    print("!!!!",feature_map_sum.shape)
    plt.imshow(feature_map_sum)
    plt.savefig('fea_map/kitti/0/'+'feature_map_sum.png')

if __name__ == '__main__':
    train()

max_stars_repo_name: Danny-proton/gwc-bts-deeplab

max_issues_repo_name: Danny-proton/gwc-bts-deeplab

max_forks_repo_name: Danny-proton/gwc-bts-deeplab