text:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sb
import ppscore as pps
import sklearn
from sklearn.ensemble import AdaBoostRegressor, RandomForestRegressor
from sklearn.preprocessing import PolynomialFeatures, OneHotEncoder
from sklearn.linear_model import LinearRegression, RANSACRegressor, Lasso, ElasticNet, Ridge
from sklearn.metrics import mean_squared_error, r2_score
import statsmodels.api as sm
import statsmodels.formula.api as smf
from sklearn.tree import DecisionTreeRegressor
from sklearn.compose import make_column_transformer
from sklearn.model_selection import train_test_split, learning_curve, ShuffleSplit
plt.style.use('dark_background')

df = pd.read_csv('fifa.csv', sep = ';')
new_df = df.sort_values(by = 'overall', ascending = False)
pd.options.display.float_format = '{:,.4f}'.format

class Eda():
    def __init__():
        '''Initializing the class'''

    def nations(country):
        '''
            Function to return countries that most players represent in fifa
                Args:
                    pandas.Series or DataFrame column
                ___________
                Returns:
                    pandas.DataFrame object
        '''
        nat = df[country].value_counts().head(10)
        nat = nat.to_frame(name = 'Total Players')
        return nat

    def popular_clubs(x):
        '''
            Function to return popular clubs in fifa
                Args:
                    pandas.Series or DataFrame column
                ___________
                Returns:
                    pandas.DataFrame object
        '''
        clb = df[x].value_counts().head(40)
        clb = clb.to_frame(name = 'Total Players')
        return clb

    def young_player(name):
        '''
            Function to return youngest players with highest potential
                ________________
                Returns:
                    pandas.DataFrame
        '''
        return df.nsmallest(10, columns = name).sort_values(by = 'potential', ascending = False)

    def old_squad(x,y):
        '''
            Function to return the team with the average oldest squad in fifa
                ________________
                Returns:
                    pandas.DataFrame
        '''
        ols = df.groupby([x])[y].mean().sort_values(ascending = False).head(10)
        ols = ols.to_frame(name = 'Age')
        return ols

    def best_cdm(pos, x):
        '''
            Function to return the best central defensive midfielders
                _______________
                Returns:
                    pandas.DataFrame
        '''
        return df[df[pos] == x][['name','age','team','overall','potential']].head(10)

    def best_str(pos, x):
        '''
            Function to return the best strikers
                ____________
                Returns:
                    pandas.DataFrame
        '''
        return df[df[pos] == x][['name','age','team','overall','potential']].head(10)

    def best_lw(pos, x):
        '''
            Function to return best left-wing forwards
                __________
                Returns:
                    pandas.DataFrame
        '''
        return new_df[new_df[pos] == x][['name','age','team','overall','potential']].head(10)

    def best_rw(pos, x):
        '''
            Function to return best right-wing forwards\
                _________
                Returns:
                    pandas.DataFrame
        '''
        return new_df[new_df[pos] == x][['name','age','team','overall','potential']].head(10)

    def best_cf(pos, x):
        '''
            Function to return best center forwards
                ___________
                Returns:
                    pandas.DataFrame
        '''
        return df[df[pos] == x][['name','age','team','overall','potential','age']].head(10)

    def best_cam(pos, x):
        '''
            Function to return best centre attacking midfielders
                __________
                Returns:
                    pandas.DataFrame
        '''
        return df[df[pos] == x][['name','age','team','overall','potential']].head(10)

    def country_squad(country):
        '''
            Function to return country team sheet
                Args:
                    pandas.Series or DataFrame column
                _____________
                Returns:
                    pandas.DataFrame
        '''
        return df.groupby('nationality').get_group(country).drop(['player_id'], axis = 1).head(23)

    def player_info(name):
        '''
            Function to return player information
                Args:
                    Player Name
                ____________
                Returns:
                    pandas.core.frame.DataFrame
        '''
        return df.loc[df['name']==name]

    def player_rate(x):
        '''
            Function for return player information based on overall ratings
                Args:
                    int: Player rating
                ____________
                Returns:
                    pandas.core.frame.DataFrame
        '''
        return df[df['overall'] == x][['name','overall','potential','position','team']].head(10)

    def player_age(x):
        '''
            Function for return player information based on age
                Args:
                    int: Player age
                ____________
                Returns:
                    pandas.core.frame.DataFrame
        '''
        return df[df['age'] == x][['name','nationality','team','position','age','overall']].head(10)

    def high_pot(x):
        '''
            Function to return players with highest potential
                Args:
                    pandas.Series or pandas.DataFrame column
                ____________
                Returns:
                    pandas.DataFrame
        '''
        return df.sort_values(x, ascending = False)[['name','age','team','position','overall','potential']].head(20)

    def pps_corr(x):
        '''
            Function for calculating the Predictive Power Score
                Args:
                    pandas.Series or DataFrame
                ___________
                Returns:
                    pandas.DataFrame or list of PPS dicts
        '''
        return pps.matrix(x)

    def pps_heatmap(df):
        '''
            Function for calculating the Predictive Power Score and plotting a heatmap
                Args:
                    Pandas DataFrame or Series object
                __________
                Returns:
                    figure
        '''
        nuff = pps.matrix(df)
        nuff1 = nuff[['x', 'y', 'ppscore']].pivot(columns='x', index='y', values='ppscore')
        plt.figure(figsize = (15, 8))
        ax = sb.heatmap(nuff1, vmin=0, vmax=1, cmap="Oranges_r", linewidths=0.5, annot=True)
        ax.set_title("PPS matrix")
        ax.set_xlabel("feature")
        ax.set_ylabel("target")
        return ax

    def pairplot(x):
        '''
            Function for plotting a scatter matrix
                Args:
                    pandas.DataFrame
                ___________
                Returns:
                    figure
        '''
        plt.figure(figsize=(15,8))
        ax = sb.pairplot(x, height = 2.5)
        return ax

    def corr(x):
        '''
            Function for computing pairwise correlation of columns.
                Args:
                    pandas.DataFrame
                ____________
                Returns:
                    Correlation Matrix
        '''
        return x[['player_id','overall','name','team','age','hits','potential']].corr()

    def descriptive(x):
        '''
            Function for generating descriptive statistics.
                Args:
                    pandas.Series or pandas.DataFrame
                ________________
                Returns:
                    Summary statistics of the Series or Dataframe provided
        '''
        return x.describe()

    def heatmap(x):
        '''
            Function to return correlation matrix on a heatmap
                Args:
                    pandas.Series or DataFrame
                _____________
                Returns:
                    figure
        '''
        plt.figure(figsize = (12, 8))
        sb.heatmap(x.corr(), annot = True,cmap = 'Oranges')
        plt.show()

    def linear_regression(feature, target, z):
        '''
            Function for fitting a linear model with coefficients w = (w1, …, wp) to minimize the residual sum of squares between the
            observed targets in the dataset, and the targets predicted by the linear approximation

                Args:
                    feature variable, target variable, prediction value(int)
                ________________
                Returns:
                    array: Prediction using the linear model
                    figure:
        '''
        ## Instantiating the algorithm ##
        model = LinearRegression()
        ## preparing the data and reshaping ##
        x = df[feature].values.reshape(-1, 1)
        y = df[target].values
        ## fitting the linear model ##
        model.fit(x, y)
        ## plotting the model ##
        plt.figure(figsize = (15, 8))
        sb.regplot(x, y)
        plt.xlabel(feature)
        plt.ylabel(target)
        plt.title('Regression Plot')
        plt.show()

        return model.predict(np.array([z]).reshape(-1, 1))

    def regression_performance(feature, target, size, depth = None, model = str):
        '''
        Function to calculate the performance of a linear regression model.
            Args:
                feature variable, target variable, model type, test size, maximum depth of tree if model == tree
            _________________
            Returns:
                ❅ Mean squared error: (float)
                ❅ Coefficient of Determination R²: (float)
                ❅ Figure
        '''
        ## Arranging data and reshaping ##
        x = df[feature].values.reshape(-1, 1)
        y = df[target].values
        ## Splitting data into training and testing ##
        x_train,x_test,y_train,y_test = train_test_split(x,y,test_size = size, random_state = 0)
        ## Instantiating the model ##
        if model == 'linear':
            model = LinearRegression()
        elif model == 'ransac':
            model = RANSACRegressor()
        elif model == 'tree':
            model = DecisionTreeRegressor(max_depth=depth)
        ## fitting the model ##
        model.fit(x_train, y_train)
        ## prediction results ##
        y_train_pred = model.predict(x_train)
        y_test_pred = model.predict(x_test)
        ## Printing the MSE and the R² score ##
        print('Mean Squared Error: Training Data %.2f' %(mean_squared_error(y_train, y_train_pred)))
        print('Mean Squared Error: Test Data %.2f' %(mean_squared_error(y_test, y_test_pred)))
        print('Coefficient of Determination R²: Training Data %.2f' %(r2_score(y_train, y_train_pred)))
        print('Coefficient of Determination R²: Test Data %.2f' %(r2_score(y_test, y_test_pred)))
        ## Residual Analysis ##
        plt.figure(figsize = (10, 8))
        plt.scatter(y_train_pred, y_train_pred - y_train, c = 'lightgoldenrodyellow', label = 'Training data')
        plt.scatter(y_test_pred, y_test_pred - y_test, c = 'tomato', label = 'Test data')
        plt.hlines(y = 0, xmin = -10, xmax = 250, lw = 2, color = 'snow')
        plt.xlabel('Predicted Values')
        plt.ylabel('Residuals')
        plt.legend(loc = 'best')
        plt.show()

    def ransac_regression(feature, target):
        '''
            Function for fitting data to robust linear estimation using RANSAC algorithm
                Args:
                    feature variable, target variable
                _____________
                Returns:
                    ❅ figure
        '''
        ## Arranging data ##
        x = df[feature].values.reshape(-1, 1)
        y = df[target].values
        ## Instantiating the algorithm ##
        ransac = RANSACRegressor()
        ## Fitting the model ##
        ransac.fit(x, y)
        ## Adding inlier_mask ##
        inlier = ransac.inlier_mask_
        outlier = np.logical_not(inlier)
        ## Adjusting x-axis to plot lines ##
        line_x = np.arange(x.min(),x.max())[:,np.newaxis]
        line_y_ransac = ransac.predict(line_x.reshape(-1, 1))
        ## Plotting figure ##
        plt.figure(figsize = (10, 8))
        plt.scatter(x[inlier], y[inlier], c = 'darkorange', marker = 'o', label = 'Inliers')
        plt.scatter(x[outlier], y[outlier], c = 'limegreen', marker = "*", label = 'Outliers')
        plt.plot(line_x, line_y_ransac, color = 'navajowhite')
        plt.legend(loc = 'best')
        plt.xlabel(feature)
        plt.ylabel(target)
        plt.show()

    def robust_regression(feature, target):
        '''
            Function to robustly fit a linear model to faulty data using the RANSAC algorithm.
                Args:
                    feature variable, target variable
                ____________
                Returns:
                    ❅ float: Estimated Coefficients
                    ❅ figure
        '''
        ## Arranging data ##
        x = df[feature].values.reshape(-1, 1)
        y = df[target].values
        ## Instantiating the model ##
        model = LinearRegression()
        model.fit(x, y)
        ## Instantiating the RANSAC algorithm ##
        ransac = RANSACRegressor()
        ransac.fit(x, y)
        ## Adding inlier_mask ##
        inlier = ransac.inlier_mask_
        outlier = np.logical_not(inlier)
        ## Adjusting x-axis to plot lines ##
        line_x = np.arange(x.min(),x.max())[:,np.newaxis]
        line_y = model.predict(line_x)
        line_y_ransac = ransac.predict(line_x.reshape(-1, 1))
        ## Printing the coefficients ##
        print('Estimated Coefficients(Linear Regression, RANSAC Regression):')
        print(model.coef_, ransac.estimator_.coef_)
        ## Plotting the figure ##
        lw = 2
        plt.figure(figsize = (10, 8))
        plt.scatter(x[inlier], y[inlier], c = 'darkorange', marker = 'o', label = 'Inliers')
        plt.scatter(x[outlier], y[outlier], c = 'limegreen', marker = "*", label = 'Outliers')
        plt.plot(line_x, line_y, color='red', linewidth=lw, label='Linear regressor')
        plt.plot(line_x, line_y_ransac, color = 'navajowhite', linewidth = lw, label ='RANSAC Regressor')
        plt.legend(loc = 'best')
        plt.xlabel(feature)
        plt.ylabel(target)
        plt.show()

    def multiple_regression(features, target):
        '''
            Function for fitting model with Multiple Linear Regression using statsmodels.api module
                Args:
                    pd.Series or list
                _______________
                Returns:
                    pd.DataFrame: Ordinary Least Squares Regression Results
        '''
        ## Arranging data ##
        x = df[features]
        y = df[target]
        ## Adding the constant ##
        x_constant = sm.add_constant(x)
        pd.DataFrame(x_constant)
        ## Fit and summary ##
        model = sm.OLS(y, x_constant)
        lr = model.fit()
        lr.summary()
        ## Predicted Values ##
        print('Predicted Values:', lr.predict())
        return lr.summary()

    def polynomial_transformation(feature, target, degree):
        '''
            Function to generate polynomial and interaction features.
                Args:
                    feature variable, target variable, degree of transformation
                _______________
                Returns:
                    ❅ float: Coefficient of Determination
                    ❅ figure
        '''
        ## Arranging data ##
        x = df[feature].values
        y = df[target].values
        ## Instantiating the algorithms ##
        poly = PolynomialFeatures(degree=degree)
        lr = LinearRegression()
        ## Fitting and transforming the data ##
        poly_x1 = poly.fit_transform(x.reshape(-1, 1))
        lr.fit(poly_x1, y.reshape(-1, 1))
        ## Making predictions ##
        y_pred1 = lr.predict(poly_x1)
        ## Plotting the figure ##
        plt.figure(figsize = (15, 8))
        plt.scatter(x, y)
        plt.plot(x, y_pred1, color = 'gold')
        plt.title('Polynomial Transformation with Linear Regression')
        print('R² Score:{:.4f}'.format(r2_score(y, y_pred1)))
        plt.show()

    def decision_tree(feature, target, depth):
        '''
            Function for fitting data to a decision tree regressor
                Args:
                    feature variable, target variable, maximum depth of the tree
                ______________
                Returns:
                    ❅ figure: Decision Tree Regressor
        '''
        ## Instantiating the algorithm ##
        tree = DecisionTreeRegressor(max_depth = depth)
        ## Arranging the data ##
        x = df[[feature]].values
        y = df[[target]].values
        ## Fitting the data ##
        tree.fit(x, y)
        ## Sorting the data ##
        sort_idx = x.flatten().argsort()
        ## Plotting the data ##
        plt.figure(figsize = (10, 8))
        plt.scatter(x[sort_idx], y[sort_idx])
        plt.plot(x[sort_idx], tree.predict(x[sort_idx]), color = 'gold')
        plt.xlabel(feature)
        plt.ylabel(target)
        return tree.predict(x[sort_idx])

    def train_test(features, labels, target, size, model = str, ratio = None, depth = None):
        """
            Function for preprocessing the data, splitting data into training and test data and calculate the accuracy of the model

            Parameters:
                features: features for predicting target variable
                labels:   Categorical values that require preprocessing
                target:   target variables
                size:     Test size after splitting data
                model:    Machine learning algorithm
                ratio:    valid only if model is ElasticNet
                          The ElasticNet mixing parameter
                depth:    depth of the tree. Valid only if model = tree and ada
            ______________________

            Returns:
                Mean Squared Error of the training and the test data
                Coefficent of Determination R² of the training and the test data
                Residual Analysis figure
        """
        ## Data Preprocessing ##
        column_trans = make_column_transformer(
        (OneHotEncoder(), labels),remainder = 'passthrough')
        X = df[features]
        y = df[target]
        x = column_trans.fit_transform(X)
        ## Splitting data for training and testing
        x_train,x_test,y_train,y_test = train_test_split(x,y,test_size = size, random_state = 0)
        ## Instantiating the model ##
        if model == 'linear':
            model = LinearRegression()
        elif model == 'ransac':
            model = RANSACRegressor()
        elif model == 'tree':
            model = DecisionTreeRegressor(max_depth=depth)
        elif model == 'elastic':
            model = ElasticNet(l1_ratio=ratio, alpha = 1, max_iter=100)
        elif model == 'ridge':
            model = Ridge()
        elif model == 'ada':
            model = AdaBoostRegressor(DecisionTreeRegressor(max_depth = depth),
                           n_estimators = 500, random_state = 42)
        elif model == 'forest':
            model = RandomForestRegressor(n_estimators=500, criterion='mse',
                                   random_state=42, n_jobs=1)
        ## fitting the model ##
        model.fit(x_train, y_train)
        ## prediction results ##
        y_train_pred = model.predict(x_train)
        y_test_pred = model.predict(x_test)
        ## Printing the MSE and the R² score ##
        print('Mean Squared Error: Training Data %.2f' %(mean_squared_error(y_train, y_train_pred)))
        print('Mean Squared Error: Test Data %.2f' %(mean_squared_error(y_test, y_test_pred)))
        print('Coefficient of Determination R²: Training Data %.2f' %(r2_score(y_train, y_train_pred)))
        print('Coefficient of Determination R²: Test Data %.2f' %(r2_score(y_test, y_test_pred)))
        ## Residual Analysis ##
        plt.figure(figsize = (10, 8))
        plt.scatter(y_train_pred, y_train_pred - y_train, c = 'lightgoldenrodyellow', label = 'Training data')
        plt.scatter(y_test_pred, y_test_pred - y_test, c = 'tomato', label = 'Test data')
        plt.hlines(y = 0, xmin = -10, xmax = 250, lw = 2, color = 'snow')
        plt.xlabel('Predicted Values')
        plt.ylabel('Residuals')
        plt.title('Residual Analysis')
        plt.legend(loc = 'best')
        plt.show()

    def learning_curve(features, labels, target, estimator, title, ylim = None, cv = None, n_jobs = 1, train_sizes = np.linspace(.1, 1.0, 5)):
        '''
            Function to determine cross-validated training and test scores for different training set sizes
                Args:
                    features: list or columns
                    labels : list; columns with categorical values
                    target: pd.Series or str; target variable
                    estimator: model algorithm
                    cv: Cross validation set
            ______________________
            Returns:
                ❅ Train and Test mean scores
                ❅ Learning Curves Figure
        '''
        ## Data Preprocessing ##
        column_trans = make_column_transformer(
        (OneHotEncoder(), labels),remainder = 'passthrough')
        X = df[features]
        y = df[target]
        x = column_trans.fit_transform(X)
        ## plotting figure ##
        plt.figure(figsize = (10, 8))
        plt.title(title)
        if ylim is not None:
            plt.ylim(*ylim)
        plt.xlabel('Training Examples')
        plt.ylabel('Score')
        ## Model Validation with Learning curves
        train_sizes, train_scores, test_scores = learning_curve(
            estimator, x, y, cv = cv, n_jobs = n_jobs, train_sizes =train_sizes)
        train_scores_mean = np.mean(train_scores, axis = 1)
        train_scores_std = np.std(train_scores, axis = 1)
        test_scores_mean = np.mean(test_scores, axis = 1)
        test_scores_std = np.std(test_scores, axis = 1)
        plt.grid()
        plt.fill_between(train_sizes, train_scores_mean - train_scores_std,
                         train_scores_mean + train_scores_std, alpha = 0.1, color = 'r')
        plt.fill_between(train_sizes, test_scores_mean - test_scores_std,
                         test_scores_mean + test_scores_std, alpha = 0.1, color = 'g')
        plt.plot(train_sizes, train_scores_mean, 'o-', color = 'r', label = 'Training score')
        plt.plot(train_sizes, test_scores_mean, 'o-', color = 'g', label = 'Cross Validation Score')
        plt.legend(loc = 'best')
        plt.show()
        return train_scores_mean, test_scores_mean

reason: import numpy,import statsmodels


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

# In[1]:

import pandas as pd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
import sklearn

# The datas used in this notebook were created by the previous notebook (Part2_Data_Cleaning).

# ## Import the Cleaned Dataset

# In[2]:

pollution_o = pd.read_excel('pollution_outliers.xlsx')
pollution_No = pd.read_excel('pollution_NoOutliers.xlsx')

# In[3]:

pollution_o.head()

# In[4]:

pollution_No.head()

# Looking at both these tables, we can see that both the datasets have a column that reiterate the index. We will remove these columns.

# In[5]:

pollution_o.drop(columns = ['Unnamed: 0'], inplace = True)
pollution_No.drop(columns = ['Unnamed: 0'], inplace = True)

# In[6]:

pollution_o.head()

# In[7]:

pollution_No.head()

# ---

# ## Creating x and y parts

# ##### Dataset with Outliers

# In[8]:

x_o = pollution_o.iloc[:, :-1].values
x_o

# In[9]:

y_o = pollution_o.iloc[:, -1].values
y_o

# In[10]:

x_o.shape

# In[11]:

y_o.shape

# We can see that the x-set has 11 columns (all the features) and 41,757 rows, while the y-set has the same number of rows but only 1 columns (the target variable of PM 2.5 levels).

# ---

# ##### Dataset without Outliers

# In[12]:

x_no = pollution_No.iloc[:, :-1].values
x_no

# In[13]:

y_no = pollution_No.iloc[:, -1].values
y_no

# In[14]:

x_no.shape

# In[15]:

y_no.shape

# This is the dataset without outliers. The x-set has 11 columns for all the features and 39,061 rows (less than the previous dataset because there are no outliers in this data). The y-set has the same number of observations and 1 column for the target variable of pollution levels.

# ---

# ## Splitting Dataset

# ##### Dataset with Outliers

# In[16]:

from sklearn.model_selection import train_test_split

x_o_train, x_o_test, y_o_train, y_o_test = train_test_split(x_o, y_o, test_size = 0.2, random_state = 4)
print('Train set:', x_o_train.shape, y_o_train.shape)
print('Test set:', x_o_test.shape, y_o_test.shape)

# This is the dataset with outliers.
#
# The training set has 33,405 observations and the test set has 8,352 rows. This is a 80:20 split of the full dataset with outliers.

# ---

# ##### Dataset without Outliers

# In[17]:

x_no_train, x_no_test, y_no_train, y_no_test = train_test_split(x_no, y_no, test_size = 0.2, random_state = 4)
print('Train set:', x_no_train.shape, y_no_train.shape)
print('Test set:', x_no_test.shape, y_no_test.shape)

# This is the dataset without outliers.
#
# The 80:20 split creates fewer observations for the x and y sets because this dataset is smaller (since outliers have been removed). The training set has 31,903 observations and the test set has 7,976 rows.

# ---

# ## Linear Regression

# ##### Data with Outliers

# In[18]:

from sklearn.linear_model import LinearRegression

linear = LinearRegression()
linear.fit(x_o_train, y_o_train)

# In[19]:

from sklearn.metrics import r2_score, mean_squared_error
from mlxtend.evaluate import bias_variance_decomp

y_o_train_pred = linear.predict(x_o_train)
mse = mean_squared_error(y_o_train, y_o_train_pred)
r2 = r2_score(y_o_train, y_o_train_pred)
mse, bias, var = bias_variance_decomp(linear, x_o_train, y_o_train, x_o_test, y_o_test, loss='mse', num_rounds=200, random_seed=1)

print("The model performance for training set")
print('MSE is {}'.format(mse))
print('R2 score is {}'.format(r2))
print("\n")

y_o_test_pred = linear.predict(x_o_test)
mse_test = mean_squared_error(y_o_test, y_o_test_pred)
r2_test = r2_score(y_o_test, y_o_test_pred)

print("The model performance for test set")
print('MSE is {}'.format(mse_test))
print('R2 score is {}'.format(r2_test))
print("\n")
print('Model Bias: %.3f' % bias)
print('Model Variance: %.3f' % var)

# The regression model using data with outliers, for the training set the r-square score is 0.261, which means that this model explains 26.1% of the total variation in the response variable of the pollution levels.
#
# As for test set, it appears to be slightly higher, the r-square is 0.264, the model explains 26.4% of the total variation, only 26.4% of the data fit the regression model.
#
# These r-squared values are pretty low. We have also displayed the bias and variance of this model. We can see that this model has extremely high bias and low variance. A high bias means that the accuracy of the predictions is not good (as seeen by the low r-squared value). Variance is the difference between many model's predictions. A high variance may lead to overfitting.
# Since we have low variance and high bias, our model is underfitting. The predictions are close or similar to each other but they are inaccurate. When the mdoel is underfitting, we can try polynomial regression to see if that improves the fit.

# In[20]:

model_o_lin = linear.fit(x_o_train, y_o_train)
print('intercept', model_o_lin.intercept_)
print('coefficients', model_o_lin.coef_)

# Looking at the above result, we can form a linear regression equation:
#
# ${pm2.5} = 97 + (-3.70) * {month} + 60.14 * {dewTemp} + (-76.53) * {temp} + (-15.69) * {pres} + (-9.59) * {windSpeed} + (-2.80) * {cumSnow} + (-8.72) * {cumRain} + (-13.14) * {NE} + (-16.15) * {NW} + 16.51 * {SE} + 12.78 * {SW} $

# This shows that a one-month increase leads to a decrease of pm2.5 levels by 3.70 micrograms per cubic meter of air. An increase of one degree celcius in the dew point temperature leads to a decrease of 60.14 micrograms per cubic meter of air. The rest of the coefficients can be interpreted similarly.
#
# However, it can be seen that the temperature has the largest influence on the level of PM 2.5 where a one degree celcius increase in the temperature reduces the PM 2.5 levels by 77 micrograms per cubic meter of air.
#
# Interesting finding: It appears that the general temporal change is that as the month increases, the pollution level decreases by 3.70 micrograms per cubic meter of air in Beijing. It would be increasing to see if there are similar trends in other countries. This finding also has real-life influences because it could advice policymakers on the way that they can/should create pollution laws that dictate people's behaviours. For instance, policymakers can create polciies that require people to use less fuel or diesel-consuming vehicles during the earlier months of the year.

# ##### Data without outliers

# In[21]:

from sklearn.linear_model import LinearRegression
linear = LinearRegression()
linear.fit(x_no_train, y_no_train)

# In[22]:

from sklearn.metrics import r2_score, mean_squared_error

y_no_train_pred = linear.predict(x_no_train)
mse = mean_squared_error(y_no_train, y_no_train_pred)
r2 = r2_score(y_no_train, y_no_train_pred)
mse, bias, var = bias_variance_decomp(linear, x_no_train, y_no_train, x_no_test, y_no_test, loss='mse', num_rounds=200, random_seed=1)

print("The model performance for training set")
print('MSE is {}'.format(mse))
print('R2 score is {}'.format(r2))
print("\n")

y_no_test_pred = linear.predict(x_no_test)
mse_test = mean_squared_error(y_no_test, y_no_test_pred)
r2_test = r2_score(y_no_test, y_no_test_pred)

print("The model performance for tet set")
print('MSE is {}'.format(mse_test))
print('R2 score is {}'.format(r2_test))
print("\n")

print('Model Bias: %.3f' % bias)
print('Model Variance: %.3f' % var)

# The regression model above is using data without outliers, for the training set the r-square score is 0.269, which means that this model explains 26.9% of the total variation in the response variable.
#
# As for test set, it appears to be slightly higher, the r-square is 0.260, the model explains 26% of the total variation in the response variable in pollution level.
#
# Overall, the difference in the r-square value in the first model and the second is close to 1% which is not a big difference. Hence, we cannot make a conclusion about which dataset would make a better predictive model.
#
# Since the r-squared for both models are still very small, which is around 25% to 27% we will try plotting with polynomial model to see if it can create a better prediction than the regression.
#
# Looking at the bias and variance values, we see that the bias is smaller for the model without outliers indicating better accuracy of the predictions. However, there is still clear high bias and low variance showing that the model is underfitting for this data. Hence, we should try a polynomial regression.

# In[23]:

model_no_lin = linear.fit(x_no_train, y_no_train)
print('intercept', model_no_lin.intercept_)
print('coefficients', model_no_lin.coef_)

# Looking at the above result, we can form a linear regression equation:
#
# ${pm2.5} = 92 + (-1.78) * {month} + 51.7 * {dewTemp} + (-62.74) * {temp} + (-12.49) * {pres} + (-9.06) * {windSpeed} + (-2.04) * {cumSnow} + (-7.73) * {cumRain} + (-13.41) * {NE} + (-15.06) * {NW} + 18.08 * {SE} + 10.38 * {SW} $

# This shows that a one-month increase leads to a decrease of pm2.5 levels by 1.78 micrograms per cubic meter of air. An increase of one degree celcius in the dew point temperature leads to a decrease of 51.7 micrograms per cubic meter of air. The rest of the coefficients can be interpreted similarly.
#
# However, it can be seen that the temperature has the largest influence on the level of PM 2.5 where a one degree celcius increase in the temperature reduces the PM 2.5 levels by 62.74 micrograms per cubic meter of air. The dew point temperature also has a large effect on the pm 2.5 levels where a one degree celcius increase in the dew point temperature leads to a 51.7 micrograms per cubic meter of air increase in pm 2.5 levels.
#
# Interesting finding: Similar to the previous equation, it appears that the general temporal change is that as the month increases, the pollution level decreases by 1.78 micrograms per cubic meter of air in Beijing. It would be increasing to see if there are similar trends in other countries. This finding also has real-life influences because it could advice policymakers on the way that they can/should create pollution laws that dictate people's behaviours. For instance, policymakers can create policies that require people to use less fuel or diesel-consuming vehicles during the earlier months of the year. Also, it can be seen that the wind directions also have an influence on pollution levels. We would prefer the ind to be coming from the north east or north west because the trend is for pollution levels to decrease in these cases in Beijing. This also gives us an idea of what areas are most polluted (the south east and south west of Beijing). This information gives policymakers an idea of the places they should creatte pollution laws.

# ---

# ## Polynomial Regression
#
# ##### Data with Outliers

# In[24]:

from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_squared_error, r2_score

poly_features_o = PolynomialFeatures(degree = 3)
x_o_train_poly = poly_features_o.fit_transform(x_o_train)
poly_model_o = LinearRegression()
poly_model_o.fit(x_o_train_poly, y_o_train)

# In[25]:

y_o_train_predict = poly_model_o.predict(x_o_train_poly)
rmse_train_o = np.sqrt(mean_squared_error(y_o_train, y_o_train_predict))
r2_train_o = r2_score(y_o_train, y_o_train_predict)

print("The model performance for training set")
print('MSE is {}'.format(rmse_train_o))
print('R2 score is {}'.format(r2_train_o))
print("\n")

y_o_test_predict = poly_model_o.predict(poly_features_o.fit_transform(x_o_test))
rmse_test_o = np.sqrt(mean_squared_error(y_o_test, y_o_test_predict))
r2_test_o = r2_score(y_o_test, y_o_test_predict)
mse, bias, var = bias_variance_decomp(poly_model_o, x_o_train_poly, y_o_train, poly_features_o.fit_transform(x_o_test), y_o_test, loss='mse', num_rounds=200, random_seed=1)

print("The model performance for test set")
print('MSE is {}'.format(rmse_test_o))
print('R2 score is {}'.format(r2_test_o))
print("\n")

print('Model Bias: %.3f' % bias)
print('Model Variance: %.3f' % var)

# The polynomial regression model is using data with outliers, for the training set the r-square score is 0.488, which means that this model explains 48.8% of the total variation in the target variable.
# As for test set, it appears to be almost the same, the r-square is 0.485, the model explains 48.5% of the total variation in the target variable.
#
# Compared with the linear regression model, the r-square in the polynomial model have a higher r-square value. It increases by 20% for both train and test set. Looking at the values of bias and variance, we can see that this model (like the others) has high bias and low variance. However, we can see that the bias has decreased whereas the variance has increased from the linear regression model. This shows that the accuracy has increased and, while the variance has increased from the linear regresison model, it is still low. Clearly, the model is still underfitting but it is a better fit for this dataset than the linear regression model.
#
# Now will see how the results are using the data without outliers.

# ##### Data without outlier

# In[26]:

from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_squared_error, r2_score

poly_features_no = PolynomialFeatures(degree = 3)
x_no_train_poly = poly_features_no.fit_transform(x_no_train)
poly_model_no = LinearRegression()
poly_model_no.fit(x_no_train_poly, y_no_train)

# In[27]:

y_no_train_predict = poly_model_no.predict(x_no_train_poly)
rmse_train_no = np.sqrt(mean_squared_error(y_no_train, y_no_train_predict))
r2_train_no = r2_score(y_no_train, y_no_train_predict)

print("The model performance for training set")
print('MSE is {}'.format(rmse_train_no))
print('R2 score is {}'.format(r2_train_no))
print("\n")

y_no_test_predict = poly_model_no.predict(poly_features_no.fit_transform(x_no_test))
rmse_test_no = np.sqrt(mean_squared_error(y_no_test, y_no_test_predict))
r2_test_no = r2_score(y_no_test, y_no_test_predict)
mse, bias, var = bias_variance_decomp(poly_model_no, x_no_train_poly, y_no_train, poly_features_no.fit_transform(x_no_test), y_no_test, loss='mse', num_rounds=200, random_seed=1)

print("The model performance for test set")
print('MSE is {}'.format(rmse_test_no))
print('R2 score is {}'.format(r2_test_no))
print("\n")

print('Model Bias: %.3f' % bias)
print('Model Variance: %.3f' % var)

# The model above is polynomial regression without outliers. The r-squared value is 0.46 for the test set and for the training set it is 0.48. This means that in the training set, the model explains 48% of the variation in the response variable. On the other hand, for the test set, the model explains 46% of the variation in the response variable.
#
#
# Based on the r-squared only, the model of polynomial regression with outliers is better because it has a higher r-squared. For this case, it does shows that outliers are necessary to keep, there might be other factors affecting the air quality features and the pollution model, which might cause some values to be higher than others. With the current climate change situation, it is possible that factors that are not recorded in the dataset can impact overall the dataset.
#
# *However*, it is important to look at other measures. Let us make a conclusion using bias and variance. Comparing the two polynomial regressions, the model without outliers has lower bias (higher accuracy) and higher variance (better predictions) compared to the dataset with outliers. Since the bias and variance are 'relatively' closer to each oher in this dataset than the dataset with outliers, we would conclude that this model, in overall, appears to be a better fit because of these values. However, there is still a large difference between the bias and variance: high bias and variance which indicates that this model is still underfitting the dataset.

# The coefficients of the polynomial regression are difficult to interpret due to squared and interaction effects created. Hence, we will not be writing or interepreting the regression equation for polynomial regression. However, we will be visualizing them.

# ---

# ## Visualizations

# Let us compare the linear and polynomial fits for each of the two datasets.

# In[28]:

import warnings
warnings.filterwarnings('ignore')

# ##### Dataset with Outliers

# In[51]:

fig, axes = plt.subplots(1, 2, figsize = (15,6))

sns.regplot(y_o_test, y_o_test_pred, data = pollution_o, color = 'red', scatter_kws = {'color':'black'}, line_kws = {'color':'green'}, ax = axes[0]).set_title('Linear Regression: True Values v/s Predicted Values')
axes[0].set_xlabel('True Values')
axes[0].set_ylabel('Predicted Values')

sns.regplot(y_o_test, y_o_test_predict, data = pollution_o, color = 'red', scatter_kws = {'color':'black'}, line_kws = {'color':'green'}, ax = axes[1])
plt.title('Polynomial Regression: True Values v/s Predicted Values')
plt.xlabel('True Values')
plt.ylabel('Predicted Values')
plt.show()

# These two regression plots show the comparison between the linear regression fitted values and the ploynomial regression fitted values. The difference is not large (as seen by a 20% increase in R-squared). However, it can be seen that the line is better fitted for the points in the polynomial regression since it goes through the middle of the cluster of points till the end (where the true values are close to 600 micrograms per cubic meter of air). On the other hand, the linear regression line is not as good a fit because it does not recorgnize the cluster of points at the curved end. and at the bottom-left of the graph.
#
# Hence, the polynomial regression is betteere in this case because the points at the curved end are recognizeed by this regression making it s a slightly better fit.

# ---

# ##### Dataset without Outliers

# In[52]:

fig, axes = plt.subplots(1, 2, figsize = (15, 6))

sns.regplot(y_no_test, y_no_test_pred, data = pollution_o, color = 'red', scatter_kws = {'color':'black'}, line_kws = {'color':'green'}, ax = axes[0]).set_title('Linear Regression: True Values v/s Predicted Values')
plt.xlabel('True Values')
plt.ylabel('Predicted Values')

sns.regplot(y_no_test, y_no_test_predict, data = pollution_o, color = 'red', scatter_kws = {'color':'black'}, line_kws = {'color':'green'})
plt.title('Polynomial Regression: True Values v/s Predicted Values')
plt.xlabel('True Values')
plt.ylabel('Predicted Values')
plt.show()

# These two regression plots are similar to the regression plots created with the datasets that had outliers. This is because the polynomial regression, in comparison to the linear regression, accounts for most of the clusters of points whereas, the linear regression does not account for the points at the curveed end of the graph.
#
# However, the polynomial regression for this dataset is not accounting for all the curvature at the end as well as the polynomial regression of the dataset with outliers. This can be seen because the line passes over the bottom-right cluster of points (and not through it). Hence, we would recommend using the dataset with outliers. There is also more variability in the polynomial regression at the end of the graph which makes fitting a line over these points dificult. Hence, just looking at this graph, we can say that the polynomial regression with outliers is a better predictive model; however, as determined earlier, the model created without outliers has better accuracy as determined by the bias and variance.

# ---

# ## Conclusion

# After running different regression models, using both with outliers and without outliers, the results of the different regression models shows that the dataset with outliers is a better fit for the prediction of the level of pollution based on different attributes of air quality. The polynomial model with outliers has better r-squared values than the ones without outliers. The final model has an r-square value of 0.48 for the training set and r-square value of 0.46 for the test set. These r-squared values between the training and test set are not much different. This means that the final model explains about 47% of the variation on the level of pollution is explained by the model.
#
# However, this conclusion is solely based on looking at the R-squared value, which is not a good standard. When we conducted an analysis of bias and variance for these models, the bias is higher for the model with outliers in comparison to the model without outliers indicating that the model has lower accuracy. Further, the variance is higher in the model without outliers which indicates that there is a better bias and variance tradeoff in this polynomial model in comparison to the model with outleirs. Hence, looking at all statistics and standards, we would recommend using the polynomial model without outleirs because it has a better accuracy and predictive power and the difference in r-squared between the two models is too small to be significant.
#
# Additionally, there might be more variables in the dataset that is not explained, that could possibly cause the generally low r-squared
#
# ##### The used of the model in real-life situations
#
# This model can be used in different real-life situations to create a better living environment for the people, such as, be used in the purpose of creating a better and healthier environment. Some attributes of air quality are dangerous for health, therefore reducing pollution level means reducing hazard for humans and for the environment.
#
# ##### Air Pollution and Economics
#
# Economics are also related to air pollution level, infrastructure and sector contributes a lot to the economic growth of a particular country. These sectors also contribute to air pollution. Through creating a model predicting air pollution, this should be used in assessing cost and benefits of economic activity in relation to air pollution. Thus, creating policy that allows the economic activities to keep working, while keeping the air pollution low. If both still have a positive correlation, then we must see whether the cost is higher than the benefit. The higher the cost, the more consequences the earth have to face

ext: py

reason: import numpy


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

# __author__="Danqing Wang"

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

"""This file contains code to read the train/eval/test data from file and process it, and read the vocab data from file and process it"""

import os
import re
import glob
import copy
import random
import json
import collections
from itertools import combinations
import numpy as np
from random import shuffle

import torch.utils.data
import time
import pickle

from nltk.tokenize import sent_tokenize

import tools.utils
from tools.logger import *

# <s> and </s> are used in the data files to segment the abstracts into sentences. They don't receive vocab ids.
SENTENCE_START = '<s>'
SENTENCE_END = '</s>'

PAD_TOKEN = '[PAD]' # This has a vocab id, which is used to pad the encoder input, decoder input and target sequence
UNKNOWN_TOKEN = '[UNK]' # This has a vocab id, which is used to represent out-of-vocabulary words
START_DECODING = '[START]' # This has a vocab id, which is used at the start of every decoder input sequence
STOP_DECODING = '[STOP]' # This has a vocab id, which is used at the end of untruncated target sequences

# Note: none of <s>, </s>, [PAD], [UNK], [START], [STOP] should appear in the vocab file.

class Vocab(object):
    """Vocabulary class for mapping between words and ids (integers)"""

    def __init__(self, vocab_file, max_size):
        """
        Creates a vocab of up to max_size words, reading from the vocab_file. If max_size is 0, reads the entire vocab file.
        :param vocab_file: string; path to the vocab file, which is assumed to contain "<word> <frequency>" on each line, sorted with most frequent word first. This code doesn't actually use the frequencies, though.
        :param max_size: int; The maximum size of the resulting Vocabulary.
        """
        self._word_to_id = {}
        self._id_to_word = {}
        self._count = 0 # keeps track of total number of words in the Vocab

        # [UNK], [PAD], [START] and [STOP] get the ids 0,1,2,3.
        for w in [PAD_TOKEN, UNKNOWN_TOKEN,  START_DECODING, STOP_DECODING]:
            self._word_to_id[w] = self._count
            self._id_to_word[self._count] = w
            self._count += 1

        # Read the vocab file and add words up to max_size
        with open(vocab_file, 'r', encoding='utf8') as vocab_f: #New : add the utf8 encoding to prevent error
            cnt = 0
            for line in vocab_f:
                cnt += 1
                pieces = line.split("\t")
                # pieces = line.split()
                w = pieces[0]
                # print(w)
                if w in [SENTENCE_START, SENTENCE_END, UNKNOWN_TOKEN, PAD_TOKEN, START_DECODING, STOP_DECODING]:
                    raise Exception('<s>, </s>, [UNK], [PAD], [START] and [STOP] shouldn\'t be in the vocab file, but %s is' % w)
                if w in self._word_to_id:
                    logger.error('Duplicated word in vocabulary file Line %d : %s' % (cnt, w))
                    continue
                self._word_to_id[w] = self._count
                self._id_to_word[self._count] = w
                self._count += 1
                if max_size != 0 and self._count >= max_size:
                    logger.info("[INFO] max_size of vocab was specified as %i; we now have %i words. Stopping reading." % (max_size, self._count))
                    break
        logger.info("[INFO] Finished constructing vocabulary of %i total words. Last word added: %s", self._count, self._id_to_word[self._count-1])

    def word2id(self, word):
        """Returns the id (integer) of a word (string). Returns [UNK] id if word is OOV."""
        if word not in self._word_to_id:
            return self._word_to_id[UNKNOWN_TOKEN]
        return self._word_to_id[word]

    def id2word(self, word_id):
        """Returns the word (string) corresponding to an id (integer)."""
        if word_id not in self._id_to_word:
            raise ValueError('Id not found in vocab: %d' % word_id)
        return self._id_to_word[word_id]

    def size(self):
        """Returns the total size of the vocabulary"""
        return self._count

    def word_list(self):
        """Return the word list of the vocabulary"""
        return self._word_to_id.keys()

class Word_Embedding(object):
    def __init__(self, path, vocab):
        """
        :param path: string; the path of word embedding
        :param vocab: object;
        """
        logger.info("[INFO] Loading external word embedding...")
        self._path = path
        self._vocablist = vocab.word_list()
        self._vocab = vocab

    def load_my_vecs(self, k=200):
        """Load word embedding"""
        word_vecs = {}
        with open(self._path, encoding="utf-8") as f:
            count = 0
            lines = f.readlines()[1:]
            for line in lines:
                values = line.split(" ")
                word = values[0]
                count += 1
                if word in self._vocablist:  # whether to judge if in vocab
                    vector = []
                    for count, val in enumerate(values):
                        if count == 0:
                            continue
                        if count <= k:
                            vector.append(float(val))
                    word_vecs[word] = vector
        return word_vecs

    def add_unknown_words_by_zero(self, word_vecs, k=200):
        """Solve unknown by zeros"""
        zero = [0.0] * k
        list_word2vec = []
        oov = 0
        iov = 0
        for i in range(self._vocab.size()):
            word = self._vocab.id2word(i)
            if word not in word_vecs:
                oov += 1
                word_vecs[word] = zero
                list_word2vec.append(word_vecs[word])
            else:
                iov += 1
                list_word2vec.append(word_vecs[word])
        logger.info("[INFO] oov count %d, iov count %d", oov, iov)
        return list_word2vec

    def add_unknown_words_by_avg(self, word_vecs, k=200):
        """Solve unknown by avg word embedding"""
        # solve unknown words inplaced by zero list
        word_vecs_numpy = []
        for word in self._vocablist:
            if word in word_vecs:
                word_vecs_numpy.append(word_vecs[word])
        col = []
        for i in range(k):
            sum = 0.0
            for j in range(int(len(word_vecs_numpy))):
                sum += word_vecs_numpy[j][i]
                sum = round(sum, 6)
            col.append(sum)
        zero = []
        for m in range(k):
            avg = col[m] / int(len(word_vecs_numpy))
            avg = round(avg, 6)
            zero.append(float(avg))

        list_word2vec = []
        oov = 0
        iov = 0
        for i in range(self._vocab.size()):
            word = self._vocab.id2word(i)
            if word not in word_vecs:
                oov += 1
                word_vecs[word] = zero
                list_word2vec.append(word_vecs[word])
            else:
                iov += 1
                list_word2vec.append(word_vecs[word])
        logger.info("[INFO] External Word Embedding iov count: %d, oov count: %d", iov, oov)
        return list_word2vec

    def add_unknown_words_by_uniform(self, word_vecs, uniform=0.25, k=200):
        """Solve unknown word by uniform(-0.25,0.25)"""
        list_word2vec = []
        oov = 0
        iov = 0
        for i in range(self._vocab.size()):
            word = self._vocab.id2word(i)
            if word not in word_vecs:
                oov += 1
                word_vecs[word] = np.random.uniform(-1 * uniform, uniform, k).round(6).tolist()
                list_word2vec.append(word_vecs[word])
            else:
                iov += 1
                list_word2vec.append(word_vecs[word])
        logger.info("[INFO] oov count %d, iov count %d", oov, iov)
        return list_word2vec

    # load word embedding
    def load_my_vecs_freq1(self, freqs, pro):
        word_vecs = {}
        with open(self._path, encoding="utf-8") as f:
            freq = 0
            lines = f.readlines()[1:]
            for line in lines:
                values = line.split(" ")
                word = values[0]
                if word in self._vocablist:  # whehter to judge if in vocab
                    if freqs[word] == 1:
                        a = np.random.uniform(0, 1, 1).round(2)
                        if pro < a:
                            continue
                    vector = []
                    for count, val in enumerate(values):
                        if count == 0:
                            continue
                        vector.append(float(val))
                    word_vecs[word] = vector
        return word_vecs

class DomainDict(object):
    """Domain embedding for Newsroom"""
    def __init__(self, path):
        self.domain_list = self.readDomainlist(path)
        # self.domain_list = ["foxnews.com", "cnn.com", "mashable.com", "nytimes.com", "washingtonpost.com"]
        self.domain_number = len(self.domain_list)
        self._domain_to_id = {}
        self._id_to_domain = {}
        self._cnt = 0

        self._domain_to_id["X"] = self._cnt
        self._id_to_domain[self._cnt] = "X"
        self._cnt += 1

        for i in range(self.domain_number):
            domain = self.domain_list[i]
            self._domain_to_id[domain] = self._cnt
            self._id_to_domain[self._cnt] = domain
            self._cnt += 1

    def readDomainlist(self, path):
        domain_list = []
        with open(path) as f:
            for line in f:
                domain_list.append(line.split("\t")[0].strip())
        logger.info(domain_list)
        return domain_list

    def domain2id(self, domain):
        """ Returns the id (integer) of a domain (string). Returns "X" for unknow domain.
        :param domain: string
        :return: id; int
        """
        if domain in self.domain_list:
            return self._domain_to_id[domain]
        else:
            logger.info(domain)
        return self._domain_to_id["X"]

    def id2domain(self, domain_id):
        """ Returns the domain (string) corresponding to an id (integer).
        :param id: int;
        :return: domain: string
        """
        if domain_id not in self._id_to_domain:
            raise ValueError('Id not found in DomainDict: %d' % domain_id)
        return self._id_to_domain[id]

    def size(self):
        return self._cnt

class Example(object):
    """Class representing a train/val/test example for text summarization."""
    def __init__(self, article_sents, abstract_sents, vocab, sent_max_len, label, domainid=None):
        """ Initializes the Example, performing tokenization and truncation to produce the encoder, decoder and target sequences, which are stored in self.

        :param article_sents: list of strings; one per article sentence. each token is separated by a single space.
        :param abstract_sents: list of strings; one per abstract sentence. In each sentence, each token is separated by a single space.
        :param domainid: int; publication of the example
        :param vocab: Vocabulary object
        :param sent_max_len: int; the maximum length of each sentence, padding all sentences to this length
        :param label: list of int; the index of selected sentences
        """

        self.sent_max_len = sent_max_len
        self.enc_sent_len = []
        self.enc_sent_input = []
        self.enc_sent_input_pad = []

        # origin_cnt = len(article_sents)
        # article_sents = [re.sub(r"\n+\t+", " ", sent) for sent in article_sents]
        # assert origin_cnt == len(article_sents)

        # Process the article
        for sent in article_sents:
            article_words = sent.split()
            self.enc_sent_len.append(len(article_words)) # store the length after truncation but before padding
            # self.enc_sent_input.append([vocab.word2id(w) for w in article_words]) # list of word ids; OOVs are represented by the id for UNK token
            self.enc_sent_input.append([vocab.word2id(w.lower()) for w in article_words]) # list of word ids; OOVs are represented by the id for UNK token
        self._pad_encoder_input(vocab.word2id('[PAD]'))

        # Store the original strings
        self.original_article = " ".join(article_sents)
        self.original_article_sents = article_sents

        if isinstance(abstract_sents[0], list):
            logger.debug("[INFO] Multi Reference summaries!")
            self.original_abstract_sents = []
            self.original_abstract = []
            for summary in abstract_sents:
                self.original_abstract_sents.append([sent.strip() for sent in summary])
                self.original_abstract.append("\n".join([sent.replace("\n", "") for sent in summary]))
        else:
            self.original_abstract_sents = [sent.replace("\n", "") for sent in abstract_sents]
            self.original_abstract = "\n".join(self.original_abstract_sents)

        # Store the label
        self.label = np.zeros(len(article_sents), dtype=int)
        if label != []:
            self.label[np.array(label)] = 1
        self.label = list(self.label)

        # Store the publication
        if domainid != None:
            if domainid == 0:
                logger.debug("domain id = 0!")
            self.domain = domainid

    def _pad_encoder_input(self, pad_id):
        """
        :param pad_id: int; token pad id
        :return:
        """
        max_len = self.sent_max_len
        for i in range(len(self.enc_sent_input)):
            article_words = self.enc_sent_input[i]
            if len(article_words) > max_len:
                article_words = article_words[:max_len]
            while len(article_words) < max_len:
                article_words.append(pad_id)
            self.enc_sent_input_pad.append(article_words)

class ExampleSet(torch.utils.data.Dataset):
    """ Constructor: Dataset of example(object) """
    def __init__(self, data_path, vocab, doc_max_timesteps, sent_max_len, domaindict=None, randomX=False, usetag=False):
        """ Initializes the ExampleSet with the path of data

        :param data_path: string; the path of data
        :param vocab: object;
        :param doc_max_timesteps: int; the maximum sentence number of a document, each example should pad sentences to this length
        :param sent_max_len: int; the maximum token number of a sentence, each sentence should pad tokens to this length
        :param domaindict: object; the domain dict to embed domain
        """
        self.domaindict = domaindict
        if domaindict:
            logger.info("[INFO] Use domain information in the dateset!")
            if randomX==True:
                logger.info("[INFO] Random some example to unknow domain X!")
                self.randomP = 0.1
        logger.info("[INFO] Start reading ExampleSet")
        start = time.time()
        self.example_list = []
        self.doc_max_timesteps = doc_max_timesteps
        cnt = 0
        with open(data_path, 'r') as reader:
            for line in reader:
                try:
                    e = json.loads(line)
                    article_sent = e['text']
                    tag = e["tag"][0] if usetag else e['publication']
                    # logger.info(tag)
                    if "duc" in data_path:
                        abstract_sent = e["summaryList"] if "summaryList" in e.keys() else [e['summary']]
                    else:
                        abstract_sent = e['summary']
                    if domaindict:
                        if randomX == True:
                            p = np.random.rand()
                            if p <= self.randomP:
                                domainid = domaindict.domain2id("X")
                            else:
                                domainid = domaindict.domain2id(tag)
                        else:
                            domainid = domaindict.domain2id(tag)
                    else:
                        domainid = None
                    logger.debug((tag, domainid))
                except (ValueError,EOFError) as e :
                    logger.debug(e)
                    break
                else:
                    example = Example(article_sent, abstract_sent, vocab, sent_max_len, e["label"], domainid) # Process into an Example.
                    self.example_list.append(example)
                cnt += 1
                # print(cnt)
            logger.info("[INFO] Finish reading ExampleSet. Total time is %f, Total size is %d", time.time() - start, len(self.example_list))
        self.size = len(self.example_list)

        # self.example_list.sort(key=lambda ex: ex.domain)

    def get_example(self, index):
        return self.example_list[index]

    def __getitem__(self, index):
        """
        :param index: int; the index of the example
        :return
            input_pad: [N, seq_len]
            label: [N]
            input_mask: [N]
            domain: [1]
        """
        item = self.example_list[index]
        input = np.array(item.enc_sent_input_pad)
        label = np.array(item.label, dtype=int)
        # pad input to doc_max_timesteps
        if len(input) < self.doc_max_timesteps:
            pad_number = self.doc_max_timesteps - len(input)
            pad_matrix = np.zeros((pad_number, len(input[0])))
            input_pad = np.vstack((input, pad_matrix))
            label = np.append(label, np.zeros(pad_number, dtype=int))
            input_mask = np.append(np.ones(len(input)), np.zeros(pad_number))
        else:
            input_pad = input[:self.doc_max_timesteps]
            label = label[:self.doc_max_timesteps]
            input_mask = np.ones(self.doc_max_timesteps)
        if self.domaindict:
            return torch.from_numpy(input_pad).long(), torch.from_numpy(label).long(), torch.from_numpy(input_mask).long(), item.domain
        return torch.from_numpy(input_pad).long(), torch.from_numpy(label).long(), torch.from_numpy(input_mask).long()

    def __len__(self):
        return self.size

class MultiExampleSet():
    def __init__(self, data_dir, vocab, doc_max_timesteps, sent_max_len, domaindict=None, randomX=False, usetag=False):
        self.datasets = [None] * (domaindict.size() - 1)
        data_path_list = [os.path.join(data_dir, s) for s in os.listdir(data_dir) if s.endswith("label.jsonl")]
        for data_path in data_path_list:
            fname = data_path.split("/")[-1]    # cnn.com.label.json
            dataname = ".".join(fname.split(".")[:-2])
            domainid = domaindict.domain2id(dataname)
            logger.info("[INFO] domain name: %s, domain id: %d" % (dataname, domainid))
            self.datasets[domainid - 1] = ExampleSet(data_path, vocab, doc_max_timesteps, sent_max_len, domaindict, randomX, usetag)

    def get(self, id):
        return self.datasets[id]

from torch.utils.data.dataloader import default_collate
def my_collate_fn(batch):
    '''
    :param batch: (input_pad, label, input_mask, domain)
    :return:
    '''
    start_domain = batch[0][-1]
    # for i in range(len(batch)):
    #     print(batch[i][-1], end=',')
    batch = list(filter(lambda x: x[-1] == start_domain, batch))
    print("start_domain %d" % start_domain)
    print("batch_len %d" % len(batch))
    if len(batch) == 0: return torch.Tensor()
    return default_collate(batch) # 用默认方式拼接过滤后的batch数据

max_stars_repo_name: Chillee/benchmark

max_issues_repo_name: stratoes/fastNLP

max_forks_repo_name: stratoes/fastNLP


text:
from pathlib import Path
import glob
import random

import cv2
import numpy as np
from .dlib_face_landmarks import get_face_landmarks

class FacemaskAugmenter:
    def __init__(self, facemask_paths=None):
        if facemask_paths is None:
            facemask_glob_pattern = str(Path(__file__).parent / "facemask-images/*.png")
            self.facemask_paths = glob.glob(facemask_glob_pattern)
        else:
            self.facemask_paths = facemask_paths

    def augment_faces(self, image: np.ndarray) -> np.ndarray:
        face_landmarks_list = get_face_landmarks(image)

        for face_landmarks in face_landmarks_list:
            image = self._augment(image, face_landmarks)

        return image

    def _augment(self, image: np.ndarray, face_landmarks: dict) -> np.ndarray:
        img_face = cv2.cvtColor(image, cv2.COLOR_RGB2BGR)
        left_upper_chin = np.array(face_landmarks['chin'][0])
        left_lower_chin = np.array(face_landmarks['chin'][5])
        right_upper_chin = np.array(face_landmarks['chin'][16])
        right_lower_chin = np.array(face_landmarks['chin'][11])
        nose_upper = np.array(face_landmarks['nose_bridge'][1])
        center_chin = np.array(face_landmarks['chin'][8])

        facemask = cv2.imread(random.choice(self.facemask_paths), -1)
        w_split = int(facemask.shape[1] / 2)

        facemask_left = facemask[:, :w_split]
        facemask_right = facemask[:, w_split:]

        left_h, left_w = facemask_left.shape[:2]
        size_h, size_w = img_face.shape[:2]

        # Warp facemask halves to desired landmark points
        input_pts = np.float32([[0, 0], [left_w, 0],
                                [left_w, left_h],[0, left_h]])

        out_left_pts = np.float32([left_upper_chin, nose_upper,
                                   center_chin, left_lower_chin])
        LM = cv2.getPerspectiveTransform(input_pts, out_left_pts)
        left_out = cv2.warpPerspective(facemask_left, LM, (size_w, size_h),
                                       flags=cv2.INTER_LINEAR)
        left_out_no_alpha = cv2.warpPerspective(facemask_left[:, :, 0:3],
                                                LM, (size_w, size_h),
                                                flags=cv2.INTER_LINEAR)

        out_right_pts = np.float32([nose_upper, right_upper_chin,
                                    right_lower_chin, center_chin])
        RM = cv2.getPerspectiveTransform(input_pts, out_right_pts)
        right_out = cv2.warpPerspective(facemask_right, RM, (size_w, size_h),
                                        flags=cv2.INTER_LINEAR)
        right_out_no_alpha = cv2.warpPerspective(facemask_right[:, :, 0:3],
                                                 RM, (size_w, size_h),
                                                 flags=cv2.INTER_LINEAR)

        left_msk = left_out[:, :, 3]
        left_msk_inv = cv2.bitwise_not(left_msk)

        right_msk = right_out[:, :, 3]
        right_msk_inv = cv2.bitwise_not(right_msk)

        msk_inv = cv2.bitwise_and(left_msk_inv, right_msk_inv)
        out = cv2.bitwise_or(left_out_no_alpha, right_out_no_alpha)

        bg = cv2.bitwise_and(img_face, img_face, mask=msk_inv)

        result = cv2.add(bg, out)
        result = cv2.cvtColor(result, cv2.COLOR_BGR2RGB)

        return result


text:
# -*- coding: utf-8 -*-
"""
Created on Wed May 12 23:38:34 2021

@author: Bruno M. Breggia
"""

from typing import List, Tuple
import numpy as np
import matplotlib.pyplot as plt

class LagrangeBaricentrico:
    """
    Clase que implementa el metodo de Interpolacion Baricentrica
    de Lagrange, tal como se describe en el paper "Barycentric Lagrange
    Interpolation", de Jean-Paul Berrut y Lloyd Trefethen (Society for
    Industrial and Applied Mathematics, 2004).
    DOI: 10.1137/S0036144502417715

    Este paper propone un metodo de interpolacion basado en el
    interpolador de Lagrange, capaz de funcionar con complejidad O(n),
    con previo calculo de "coeficientes baricentricos". El calculo de
    dichos coeficientes puede ser de complejidad O(n^2) o menor, segun
    la distribucion de los puntos en el eje de las abscisas. El agregado
    de puntos implica una actualizacion de los coeficientes con
    complejidad O(n).
    """

    def __init__(self, puntos:List=[]):
        """
        Inicializador de la clase, se le puede pasar una lista de tuplas
        de pares de numeros, para que sean los puntos iniciales para la
        interpolacion.

        Parameters
        ----------
        puntos : List[Tuple[float,float]], optional
            Lista de puntos iniciales para el interpolador. Por defecto
            la lista inicia vacia.

        Returns
        -------
        None.

        """

        self.puntos = puntos

        # Orden del polinomio
        self.orden = len(self.puntos)-1 if self.puntos else None

        # Calculo coeficientes baricentricos
        self.w = None
        self.precalculo()

        pass

    def __call__(self, px):
        return self.interpolar(px)

    def precalculo(self):
        """
        Calculo de los coeficientes baricentricos.
        Algoritmo con complejidad O(n^2)

        """
        X = 0

        self.w = np.array([0]*len(self.puntos), np.float64)
        for i in range(len(self.w)):
            prod = 1
            for j in range(len(self.w)):
                if j!= i:
                    prod *= (self.puntos[i][X] - self.puntos[j][X])
            self.w[i] = 1/prod
        pass

    def interpolar(self, px:float) -> float:
        """
        Evalua la interpolacion en el punto px.
        Algoritmo con complejidad O(n)

        Parameters
        ----------
        px : float
            Abscisa a evaluar.

        Returns
        -------
        py : float
            Ordenada correspondiente a la abscisa provista.

        """
        X, Y = 0, 1

        L = np.prod( [ px-self.puntos[i][X]
                      for i in range(len(self.puntos)) ], 0 )

        suma = 0
        for i in range(len(self.puntos)):
            suma += self.w[i]*self.puntos[i][Y]/(px-self.puntos[i][X])

        return L*suma

    def agregar(self, px, py):
        """
        Agrega un nuevo punto de referencia al interpolador.
        Se actualizan los coeficientes baricentricos.
        La actualizacion es de complejidad O(n)

        Parameters
        ----------
        px : float
            Abscisa del punto agregado.
        py : float
            Ordenada del punto agregado.

        """
        X = 0

        # Agrego a la lista de puntos
        self.puntos.append( (px,py) )

        # Actualizo orden
        self.orden += 1

        # Actualizo coeficientes baricentricos
        prod = 1
        for i in range(len(self.w)):
            self.w[i] /= (self.puntos[i][X] - px)
            prod *= (px - self.puntos[i][X])
        self.w = np.append(self.w, [1/prod], 0)

    def pop(self):
        """
        Remueve el ultimo punto de referencia del interpolador

        """
        X = 0

        # Agrego a la lista de puntos
        px, _ = self.puntos.pop()

        # Actualizo orden
        self.orden -= 1

        # Actualizo coeficientes baricentricos
        self.w = np.delete(self.w, -1)
        for i in range(len(self.w)):
            self.w[i] *= (self.puntos[i][X] - px)

    pass

# %% main
def main():

    puntos = [(-3, 7),
              (0, -5),
              (5, 3.3),
              # (2.3, 5)
              ]

    interpolador = LagrangeBaricentrico(puntos)
    testeo(interpolador)

    interpolador.agregar(2.3, 5)
    testeo(interpolador)

    interpolador.agregar(-4, -10)
    testeo(interpolador)

    interpolador.pop()
    testeo(interpolador)

    pass

# %% testeo
def testeo(interpolador):

    # Graficacion
    fig, ax = plt.subplots(1, 1, figsize=(15, 7), sharex=True)
    fig.suptitle("Funcion interpolante", fontsize=18)
    abscisas = np.arange(-6.5,6.5,0.02)
    ordenadas = interpolador(abscisas)
    ax.plot(abscisas,
            ordenadas,
            label='Interpolacion de Lagrange',
            color='green')
    ax.legend(loc="upper right", fontsize=15)
    ax.grid()
    #plt.show()

# %% if __name__ == '__main__'
if __name__ == '__main__':
    plt.close('all')
    main()


text:
# The plot server must be running
# Go to http://localhost:5006/bokeh to view this plot

import numpy as np
from bokeh.plotting import *

def grid_example():
    N = 80

    x = np.linspace(0, 4*np.pi, N)
    y = np.sin(x)

    output_server("grid")
    l = line(x,y, color="#0000FF", tools="pan,wheel_zoom,box_zoom,reset,previewsave",
         title="line")

    N = 9
    x2 = np.linspace(-2, 2, N)
    y2 = x2**2
    sizes = np.linspace(4, 10, N)
    xpts = np.array([-.09, -.12, .0, .12, .09])
    ypts = np.array([-.1, .02, .1, .02, -.1])

    aw = annular_wedge(
        x, y, 10, 20, 0.6, 4.1,
        inner_radius_units="screen", outer_radius_units = "screen",
        color="#8888ee", tools="pan,wheel_zoom,box_zoom,reset,previewsave", title="annular_wedge",
        name="glyphs_example")
    bez = bezier(x, y, x+0.2, y, x+0.1, y+0.1, x-0.1, y-0.1,
           color="#8888ee", tools="pan,wheel_zoom,box_zoom,reset,previewsave", title="bezier")
    q = quad(x, x-0.5, y, y-0.5,
         color="#8888ee", tools="pan,wheel_zoom,box_zoom,reset,previewsave", title="quad")
    gridplot([[l,aw],[bez,q]])
    return curplot()

if __name__ == "__main__":
    grid_example()
    # open a browser
    show()

max_stars_repo_name: tswicegood/bokeh

max_issues_repo_name: csaid/bokeh

max_forks_repo_name: csaid/bokeh


text:
import h5py
import numpy as np
import os.path
from PIL import Image
from glob import glob
from skimage.transform import resize

raw_training_x_path = './data/ALL/training/images/*.png'
raw_training_y_path= './data/ALL/training/av/*.png'

raw_data_path = [raw_training_x_path, raw_training_y_path]

HDF5_data_path = './data/HDF5/'

DESIRED_DATA_SHAPE = [576, 576]

def isHDF5exists(raw_data_path, HDF5_data_path):
    for raw in raw_data_path:
        if not raw:
            continue

        raw_splited = raw.split('/')
        HDF5 = ''.join([HDF5_data_path, '/'.join(raw_splited[2:-1]), '/*.hdf5'])

        if len(glob(HDF5)) == 0:
            return False

    return True

def read_input(path):
    x = np.array(Image.open(path)) / 255.
    if x.shape[-1] == 3:
        return x
    else:
        return x[..., np.newaxis]

def read_input_av_label(path, need_one_hot=False):
    global DESIRED_DATA_SHAPE
    x = Image.open(path)
    x = np.array(x) / 255.
    x = resize(x, DESIRED_DATA_SHAPE)

    new_x_whole = np.zeros((x.shape[0], x.shape[1]))
    new_x_a = np.zeros((x.shape[0], x.shape[1]))
    new_x_v = np.zeros((x.shape[0], x.shape[1]))
    new_x_cross= np.zeros((x.shape[0], x.shape[1]))
    new_x_unknown = np.zeros((x.shape[0], x.shape[1]))
    new_x_onehot = np.zeros((x.shape[0], x.shape[1], 3)) #bg, a, v

    for row_id, row in enumerate(x):
        for col_id, elem in enumerate(row):
            if elem[0] > 0.5 or elem[1] > 0.5 or elem[2] > 0.5:
                new_x_whole[row_id, col_id] = 1.0
            else:
                new_x_onehot[row_id, col_id] = np.array((1.0, 0.0, 0.0))
            if elem[0] > 0.5 and elem[1] > 0.5 and elem[2] > 0.5:
                new_x_unknown[row_id, col_id] = 1.0
                new_x_a[row_id, col_id] = 1.0
                new_x_onehot[row_id, col_id] = np.array((0.0, 1.0, 0.0))
            elif elem[0] > 0.5:
                new_x_a[row_id, col_id] = 1.0
                new_x_onehot[row_id, col_id] = np.array((0.0, 1.0, 0.0))
            elif elem[1] > 0.5:
                new_x_cross[row_id, col_id] = 1.0
                new_x_a[row_id, col_id] = 1.0
                # new_x_v[row_id, col_id] = 1.0
                new_x_onehot[row_id, col_id] = np.array((0.0, 1.0, 0.0))
            elif elem[2] > 0.5:
                new_x_v[row_id, col_id] = 1.0
                new_x_onehot[row_id, col_id] = np.array((0.0, 0.0, 1.0))

    new_x_whole = resize(new_x_whole[..., np.newaxis], DESIRED_DATA_SHAPE)
    new_x_a = resize(new_x_a[..., np.newaxis], DESIRED_DATA_SHAPE)
    new_x_v = resize(new_x_v[..., np.newaxis], DESIRED_DATA_SHAPE)
    new_x_cross = resize(new_x_cross[..., np.newaxis], DESIRED_DATA_SHAPE)
    new_x_unknown = resize(new_x_unknown[..., np.newaxis], DESIRED_DATA_SHAPE)
    new_x_onehot = resize(new_x_onehot, DESIRED_DATA_SHAPE)

    if not need_one_hot:
        return [new_x_whole, new_x_cross, new_x_a, new_x_v,  new_x_unknown]
    else:
        return [new_x_onehot]

def preprocessData(data_path, dataset, need_one_hot=False):
    global DESIRED_DATA_SHAPE

    data_path = list(sorted(glob(data_path)))

    if data_path[0].find('mask') > 0:
        return np.array([read_input(image_path) for image_path in data_path])
    elif data_path[0].find('/av/') > 0 or data_path[0].find('/arteries-and-veins/') > 0 :
        return np.array([read_input_av_label(image_path, need_one_hot) for image_path in data_path])
    else:
        return np.array([resize(read_input(image_path), DESIRED_DATA_SHAPE) for image_path in data_path])

def createHDF5(data, HDF5_data_path, one_hot=False):
    try:
        os.makedirs(HDF5_data_path, exist_ok=True)
    except:
        pass
    if not one_hot:
        f = h5py.File(HDF5_data_path + 'data.hdf5', 'w')
    else:
        f = h5py.File(HDF5_data_path + 'data_onehot.hdf5', 'w')
    f.create_dataset('data', data=data)
    return

def prepareDataset(dataset):
    global raw_data_path, HDF5_data_path
    global DESIRED_DATA_SHAPE

    if isHDF5exists(raw_data_path, HDF5_data_path):
        return

    for raw in raw_data_path:
        if not raw:
            continue

        raw_splited = raw.split('/')
        HDF5 = ''.join([HDF5_data_path, '/'.join(raw_splited[2:-1]), '/'])

        preprocessed = preprocessData(raw, dataset)
        createHDF5(preprocessed, HDF5)

        if raw.find('/av/') > 0 or raw.find('/arteries-and-veins/') > 0:
            raw_splited = raw.split('/')
            HDF5 = ''.join([HDF5_data_path, '/'.join(raw_splited[2:-1]), '/'])

            preprocessed = preprocessData(raw, dataset, need_one_hot=True)
            createHDF5(preprocessed, HDF5, one_hot=True)

def getTrainingData(XorY, dataset, need_one_hot=False):
    global HDF5_data_path

    raw_training_x_path, raw_training_y_path = raw_data_path[:2]

    if XorY == 0:
        raw_splited = raw_training_x_path.split('/')
    else:
        raw_splited = raw_training_y_path.split('/')

    if not need_one_hot:
        data_path = ''.join([HDF5_data_path, dataset, '/', '/'.join(raw_splited[3:-1]), '/data.hdf5'])
    else:
        data_path = ''.join([HDF5_data_path, dataset, '/', '/'.join(raw_splited[3:-1]), '/data_onehot.hdf5'])
    f = h5py.File(data_path, 'r')
    data = f['data']

    return data

def getTestData(XorYorMask, dataset, need_one_hot=False):
    global HDF5_data_path

    raw_test_x_path, raw_test_y_path = raw_data_path[2:]

    if XorYorMask == 0:
        raw_splited = raw_test_x_path.split('/')
    elif XorYorMask == 1:
        raw_splited = raw_test_y_path.split('/')
    else:
        if not raw_test_mask_path:
            return None
        raw_splited = raw_test_mask_path.split('/')

    if not need_one_hot:
        data_path = ''.join([HDF5_data_path, dataset, '/', '/'.join(raw_splited[3:-1]), '/data.hdf5'])
    else:
        data_path = ''.join([HDF5_data_path, dataset, '/', '/'.join(raw_splited[3:-1]), '/data_onehot.hdf5'])
    f = h5py.File(data_path, 'r')
    data = f['data']

    return data


text:
# Utils used for the data,
import numpy as np, matplotlib.pyplot as plt, sys
import SimPEG as simpeg
import simpegMT as simpegmt
import numpy.lib.recfunctions as recFunc
from scipy.constants import mu_0
from scipy import interpolate as sciint

def getAppRes(MTdata):
    # Make impedance
    zList = []
    for src in MTdata.survey.srcList:
        zc = [src.freq]
        for rx in src.rxList:
            if 'i' in rx.rxType:
                m=1j
            else:
                m = 1
            zc.append(m*MTdata[src,rx])
        zList.append(zc)
    return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]

def rotateData(MTdata,rotAngle):
    '''
    Function that rotates clockwist by rotAngle (- negative for a counter-clockwise rotation)
    '''
    recData = MTdata.toRecArray('Complex')
    impData = rec2ndarr(recData[['zxx','zxy','zyx','zyy']],complex)
    # Make the rotation matrix
    # c,s,zxx,zxy,zyx,zyy = sympy.symbols('c,s,zxx,zxy,zyx,zyy')
    # rotM = sympy.Matrix([[c,-s],[s, c]])
    # zM = sympy.Matrix([[zxx,zxy],[zyx,zyy]])
    # rotM*zM*rotM.T
    # [c*(c*zxx - s*zyx) - s*(c*zxy - s*zyy), c*(c*zxy - s*zyy) + s*(c*zxx - s*zyx)],
    # [c*(c*zyx + s*zxx) - s*(c*zyy + s*zxy), c*(c*zyy + s*zxy) + s*(c*zyx + s*zxx)]])
    s = np.sin(-np.deg2rad(rotAngle))
    c = np.cos(-np.deg2rad(rotAngle))
    rotMat = np.array([[c,-s],[s,c]])
    rotData = (rotMat.dot(impData.reshape(-1,2,2).dot(rotMat.T))).transpose(1,0,2).reshape(-1,4)
    outRec = recData.copy()
    for nr,comp in enumerate(['zxx','zxy','zyx','zyy']):
        outRec[comp] = rotData[:,nr]
    return simpegmt.DataMT.DataMT.fromRecArray(outRec)

def appResPhs(freq,z):
    app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
    app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
    return app_res, app_phs

def skindepth(rho,freq):
    ''' Function to calculate the skindepth of EM waves'''
    return np.sqrt( (rho*((1/(freq * mu_0 * np.pi )))))

def rec2ndarr(x,dt=float):
    return x.view((dt, len(x.dtype.names)))

def makeAnalyticSolution(mesh,model,elev,freqs):
  data1D = []
  for freq in freqs:
      anaEd, anaEu, anaHd, anaHu = simpegmt.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)
      anaE = anaEd+anaEu
      anaH = anaHd+anaHu

      anaZ = anaE/anaH
      # Add to the list
      data1D.append((freq,0,0,elev,anaZ[0]))
  dataRec = np.array(data1D,dtype=[('freq',float),('x',float),('y',float),('z',float),('zyx',complex)])
  return dataRec

def plotMT1DModelData(problem,models,symList=None):

    # Setup the figure
    fontSize = 15

    fig = plt.figure(figsize=[9,7])
    axM = fig.add_axes([0.075,.1,.25,.875])
    axM.set_xlabel('Resistivity [Ohm*m]',fontsize=fontSize)
    axM.set_xlim(1e-1,1e5)
    axM.set_ylim(-10000,5000)
    axM.set_ylabel('Depth [km]',fontsize=fontSize)
    axR = fig.add_axes([0.42,.575,.5,.4])
    axR.set_xscale('log')
    axR.set_yscale('log')
    axR.invert_xaxis()
    # axR.set_xlabel('Frequency [Hz]')
    axR.set_ylabel('Apparent resistivity [Ohm m]',fontsize=fontSize)

    axP = fig.add_axes([0.42,.1,.5,.4])
    axP.set_xscale('log')
    axP.invert_xaxis()
    axP.set_ylim(0,90)
    axP.set_xlabel('Frequency [Hz]',fontsize=fontSize)
    axP.set_ylabel('Apparent phase [deg]',fontsize=fontSize)

    # if not symList:
    #   symList = ['x']*len(models)
    import plotDataTypes as pDt
    # Loop through the models.
    modelList = [problem.survey.mtrue]
    modelList.extend(models)
    if False:
        modelList = [problem.mapping.sigmaMap*mod for mod in modelList]
    for nr, model in enumerate(modelList):
        # Calculate the data
        if nr==0:
            data1D = problem.dataPair(problem.survey,problem.survey.dobs).toRecArray('Complex')
        else:
            data1D = problem.dataPair(problem.survey,problem.survey.dpred(model)).toRecArray('Complex')
        # Plot the data and the model
        colRat = nr/((len(modelList)-1.999)*1.)
        if colRat > 1.:
            col = 'k'
        else:
            col = plt.cm.seismic(1-colRat)
        # The model - make the pts to plot
        meshPts = np.concatenate((problem.mesh.gridN[0:1],np.kron(problem.mesh.gridN[1::],np.ones(2))[:-1]))
        modelPts = np.kron(1./(problem.mapping.sigmaMap*model),np.ones(2,))
        axM.semilogx(modelPts,meshPts,color=col)

        ## Data
        # Appres
        pDt.plotIsoStaImpedance(axR,np.array([0,0]),data1D,'zyx','res',pColor=col)
        # Appphs
        pDt.plotIsoStaImpedance(axP,np.array([0,0]),data1D,'zyx','phs',pColor=col)
        try:
            allData = np.concatenate((allData,simpeg.mkvc(data1D['zyx'],2)),1)
        except:
            allData = simpeg.mkvc(data1D['zyx'],2)
    freq = simpeg.mkvc(data1D['freq'],2)
    res, phs = appResPhs(freq,allData)

    stdCol = 'gray'
    axRtw = axR.twinx()
    axRtw.set_ylabel('Std of log10',color=stdCol)
    [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
    axPtw = axP.twinx()
    axPtw.set_ylabel('Std ',color=stdCol)
    [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
    axRtw.plot(freq, np.std(np.log10(res),1),'--',color=stdCol)
    axPtw.plot(freq, np.std(phs,1),'--',color=stdCol)

    # Fix labels and ticks

    yMtick = [l/1000 for l in axM.get_yticks().tolist()]
    axM.set_yticklabels(yMtick)
    [ l.set_rotation(90) for l in axM.get_yticklabels()]
    [ l.set_rotation(90) for l in axR.get_yticklabels()]
    [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
    [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
    for ax in [axM,axR,axP]:
        ax.xaxis.set_tick_params(labelsize=fontSize)
        ax.yaxis.set_tick_params(labelsize=fontSize)
    return fig

def printTime():
    import time
    print time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime())

def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
    # Find the unique locations
    # Need to find the locations
    recDataTemp = MTdata.toRecArray()
    # Check if survey.std has been assigned.
    ## NEED TO: write this...
    # Calculte and add the DET of the tensor to the recArray
    if 'det' in rxType3D:
        Zon = (recDataTemp['zxxr']+1j*recDataTemp['zxxi'])*(recDataTemp['zyyr']+1j*recDataTemp['zyyi'])
        Zoff = (recDataTemp['zxyr']+1j*recDataTemp['zxyi'])*(recDataTemp['zyxr']+1j*recDataTemp['zyxi'])
        det = np.sqrt(Zon.data - Zoff.data)
        recData = recFunc.append_fields(recDataTemp,['zdetr','zdeti'],[det.real,det.imag] )
    else:
        recData = recDataTemp

    uniLocs = rec2ndarr(np.unique(recData[['x','y','z']])).data
    mtData1DList = []
    if 'zxy' in rxType3D:
        corr = -1 # Shift the data to comply with the quadtrature of the 1d problem
    else:
        corr = 1
    for loc in uniLocs:
        # Make the receiver list
        rx1DList = []
        for rxType in ['z1dr','z1di']:
            rx1DList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(loc,2).T,rxType))
        # Source list
        locrecData = recData[np.sqrt(np.sum( (rec2ndarr(recData[['x','y','z']]).data - loc )**2,axis=1)) < 1e-5]
        dat1DList = []
        src1DList = []
        for freq in locrecData['freq']:
            src1DList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rx1DList,freq))
            for comp  in ['r','i']:
                dat1DList.append( corr * locrecData[rxType3D+comp][locrecData['freq']== freq].data )

        # Make the survey
        sur1D = simpegmt.SurveyMT.SurveyMT(src1DList)

        # Make the data
        dataVec = np.hstack(dat1DList)
        dat1D = simpegmt.DataMT.DataMT(sur1D,dataVec)
        sur1D.dobs = dataVec
        # Need to take MTdata.survey.std and split it as well.
        std=0.05
        sur1D.std =  np.abs(sur1D.dobs*std) #+ 0.01*np.linalg.norm(sur1D.dobs)
        mtData1DList.append(dat1D)

    # Return the the list of data.
    return mtData1DList

def resampleMTdataAtFreq(MTdata,freqs):
    '''
    Function to resample MTdata at set of frequencies

    '''
    # Make a rec array
    MTrec = MTdata.toRecArray().data

    # Find unique locations
    uniLoc = np.unique(MTrec[['x','y','z']])
    uniFreq = MTdata.survey.freqs
    # Get the comps
    dNames = MTrec.dtype

    # Loop over all the locations and interpolate
    for loc in uniLoc:
        # Find the index of the station
        ind = np.sqrt(np.sum((rec2ndarr(MTrec[['x','y','z']]) - rec2ndarr(loc))**2,axis=1)) < 1. # Find dist of 1 m accuracy
        # Make a temporary recArray and interpolate all the components
        tArrRec = np.concatenate((simpeg.mkvc(freqs,2),np.ones((len(freqs),1))*rec2ndarr(loc),np.nan*np.ones((len(freqs),12))),axis=1).view(dNames)
        for comp in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
            int1d = sciint.interp1d(MTrec[ind]['freq'],MTrec[ind][comp],bounds_error=False)
            tArrRec[comp] = simpeg.mkvc(int1d(freqs),2)

        # Join together
        try:
            outRecArr = recFunc.stack_arrays((outRecArr,tArrRec))
        except NameError as e:
            outRecArr = tArrRec

    # Make the MTdata and return
    return simpegmt.DataMT.DataMT.fromRecArray(outRecArr)

max_forks_repo_name: simpeg/simpegmt

reason: import numpy,from scipy