<jupyter_start><jupyter_text>First we say what query and project we are interested inThis is using the NACSOS platform<jupyter_code>pid = 234
qid = 8318
qids = [8318, 8421, 8503]
merged_query = Query.objects.get(pk=8510)
docids = set(Doc.objects.filter(query__in=qids).values_list('id',flat=True))
docs = Doc.objects.filter(id__in=docids)
p = Project.objects.get(pk=pid)

#print(q.doc_set.count())
docs.count()<jupyter_output><empty_output><jupyter_text>And we retrieve those documents that have been screened in that query, and get the majority rating<jupyter_code>import statistics

def mode_na(x):
    try:
        return statistics.mode([y for y in x if y is not np.NaN])
    except:
        return np.NaN
    
columns = [
    'id',
    'content',
    'title',
    'PY',
    'wosarticle__de',
    'wosarticle__wc'
]

# First get only the ratings from updated tags

update_tags = Tag.objects.filter(pk__in=[5598, 5597, 5577, 5562, 5551, 5550, 5632, 5630, 5629]) 
updated_dos = DocOwnership.objects.filter(
    tag__in=update_tags, 
    relevant__gt=0,
    doc__content__iregex='\w'
)
dodf = pd.DataFrame.from_dict(updated_dos.order_by('finish').values(
    'id','doc_id','user__username','relevant'
))
dodf.loc[dodf['relevant']==2,'relevant'] = 0
dodf.loc[dodf['relevant']==3,'relevant'] = np.NaN

gdf = dodf.groupby('doc_id').agg(lambda x: mode_na(x)).reset_index()[['doc_id','relevant']].rename(columns={"doc_id":"id"})
gdf['seen'] = 1
gdf.loc[pd.isna(gdf['relevant']),'seen'] = np.NaN

updated_df = gdf

ddf = pd.DataFrame.from_dict(
    Doc.objects.filter(pk__in=updated_df['id']).values(*columns)
)
updated_df = pd.merge(ddf,updated_df)
# Now get majority ratings from the rest of documents
# Majority ratings
dos = DocOwnership.objects.filter(
    query__in=qids, 
    relevant__gt=0,
    doc__content__iregex='\w'
).exclude(doc__id__in=updated_df["id"])
dodf = pd.DataFrame.from_dict(dos.order_by('finish').values(
    'id','doc_id','user__username','relevant'
))
dodf.loc[dodf['relevant']==2,'relevant'] = 0
dodf.loc[dodf['relevant']==3,'relevant'] = np.NaN

gdf = dodf.groupby('doc_id').agg(lambda x: mode_na(x)).reset_index()[['doc_id','relevant']].rename(columns={"doc_id":"id"})
gdf['seen'] = 1
gdf.loc[pd.isna(gdf['relevant']),'seen'] = np.NaN

seen_df = gdf
ddf = pd.DataFrame.from_dict(
    Doc.objects.filter(pk__in=seen_df['id']).values(*columns)
)
seen_df = pd.merge(ddf,seen_df)

seen_df = pd.concat([seen_df,updated_df])
print(seen_df.shape)
seen_df.head()
#ADDED FOR TEST
testing=False
if testing:
    seen_df = seen_df.head(1000)
    print(seen_df.shape)<jupyter_output><empty_output><jupyter_text>We have a separate list of documents in the query that have not been screened<jupyter_code>t0 = time.time()
#docs = Doc.objects.filter(pk__lt=1000,query=merged_query)
docs = Doc.objects.filter(id__in=docids)
unseen_docs = docs.filter(
    content__iregex='\w',
    wosarticle__dt__iregex='Article|Review'
).exclude(
    wosarticle__dt__iregex='book|proceedings|comment|retracted|editorial'
).exclude(
    id__in=seen_df['id']
)

unseen_df = pd.DataFrame.from_dict(list(unseen_docs.values(*columns)))
unseen_df = unseen_df[columns]
unseen_df['relevant'] = 0
unseen_df['seen'] = 0

print(unseen_df.shape)
print(time.time()-t0)
unseen_df.head()
#ADDED FOR TEST
if testing:
    unseen_df = unseen_df.head(1000)
    print(unseen_df.shape)<jupyter_output><empty_output><jupyter_text>And we merge the two, with a column "seen" recording if a document has been seen or not<jupyter_code>django.db.connection.close()
df = (pd.concat([seen_df,unseen_df])
      .sort_values('id')
      .sample(frac=1, random_state=1)
      .reset_index(drop=True)
)
y = df['relevant']

seen_index = df[df['seen']==1].index
unseen_index = df[df['seen']==0].index

print(df.shape)
df.head(10)<jupyter_output>(286144, 8)<jupyter_text>We also get the assignments of these documents to categories, and add these as further columns in our dataframeagain we get the majority ratings<jupyter_code># first get only those ratings from the updated tags
dudf = pd.DataFrame.from_dict(
    DocUserCat.objects.filter(
        doc__pk__in=updated_df['id'],category__project=p
    ).values('user__username','doc__id','category__level','category__name')
)
updated_docusers = updated_dos.values_list('doc__id',"user__username")
dudf = dudf[dudf[['doc__id', 'user__username']].apply(tuple, axis=1).isin(updated_docusers)]
dudf['v'] = 1
dudf['cname'] =  dudf['category__level'].astype(str) + " - " + dudf['category__name']#
dudf = dudf.drop(columns=['category__level','category__name'])
dudf_wide = dudf.pivot_table(index=['doc__id','user__username'],columns="cname", values="v").fillna(0)
updated_gudf = dudf_wide.groupby(['doc__id']).agg(lambda x: st.mode(x)[0][0]).reset_index()


dudf = pd.DataFrame.from_dict(
    DocUserCat.objects.filter(
        doc__pk__in=set(seen_df['id']) - set(updated_df['id']),category__project=p
    ).values('user__username','doc__id','category__level','category__name')
)
dudf['v'] = 1
dudf['cname'] =  dudf['category__level'].astype(str) + " - " + dudf['category__name']#
dudf = dudf.drop(columns=['category__level','category__name'])
dudf_wide = dudf.pivot_table(index=['doc__id','user__username'],columns="cname", values="v").fillna(0)
gudf = dudf_wide.groupby(['doc__id']).agg(lambda x: st.mode(x)[0][0]).reset_index()

#merge updated and old
gudf = pd.concat([gudf,updated_gudf])

print(gudf.shape)
gudf.tail()<jupyter_output>(1038, 4)<jupyter_text>And merge<jupyter_code>df = pd.merge(df, gudf.rename(columns={"doc__id":"id"}), how="left").fillna(0)
print(df.shape)
df.head()<jupyter_output>(286144, 11)<jupyter_text>We have a look at the distribution of relevant and not relevant documents<jupyter_code>fig, ax = plt.subplots()
n = seen_df.shape[0]
y = df['relevant']
for i, (name, group) in enumerate(seen_df.groupby('relevant')):
    ax.bar(i, group.shape[0])
    ax.text(i, group.shape[0]+n*0.02, f'{group.shape[0]/n:.0%}',ha="center")

ax.set_xticks([0,1])
ax.set_ylim(ymax=ax.get_ylim()[1]*1.05)
ax.set_xlabel('Relevance')
ax.set_ylabel('Number of documents')<jupyter_output><empty_output><jupyter_text>Now we "vectorize" the data, representing each document as a weighted count of its words and phrases (bag of words and bigrams). Or if we already did this previously, then we reload the vectorized data<jupyter_code>import pickle
import scipy
revectorize = True
y = df['relevant']
if revectorize:
    vec = TfidfVectorizer(
        ngram_range=(1,2),
        min_df=4, max_df=0.8, strip_accents='unicode', 
        max_features=10000,
        use_idf=1,
        smooth_idf=1, sublinear_tf=1,
        stop_words="english",tokenizer=lu.tokenize
    )
    vec.fit(df.loc[seen_index,'content'].astype(str))
    ab_X = vec.transform(df['content'].astype(str))
    with open ('vec.pickle','wb') as f: #changed dir from data/vec.pickle to just vec.pickle
        pickle.dump(vec, f)
    import scipy.sparse
    scipy.sparse.save_npz('X.npz', ab_X)
else:
    with open ('vec.pickle','rb') as f:
        vec = pickle.load(f)
    ab_X = scipy.sparse.load_npz('X.npz')#changed dir from data/vec.pickle to just vec.pickle
    
print(ab_X.shape)<jupyter_output>/home/galm/software/django/tmv/venv/lib/python3.6/site-packages/sklearn/feature_extraction/text.py:484: UserWarning: The parameter 'token_pattern' will not be used since 'tokenizer' is not None'
  warnings.warn("The parameter 'token_pattern' will not be used"
/home/galm/software/django/tmv/venv/lib/python3.6/site-packages/sklearn/feature_extraction/text.py:386: UserWarning: Your stop_words may be inconsistent with your preprocessing. Tokenizing the stop words generated tokens ['make'] not in stop_words.
  'stop_words.' % sorted(inconsistent))<jupyter_text>We build 3 types of classifiers, and test these with k-fold cross validation<jupyter_code>def cross_validate_models(X,y,clf_models, seen_index, n_splits=10, classes=None,
    upsample=False,roundup=False, df=None,
    stratified_k=False, test_index=None, p_threshold=None):    
    
    if stratified_k:
        label_encoder = LabelEncoder()
        kf = StratifiedKFold(n_splits=n_splits)
        kfs = kf.split(X[seen_index],label_encoder.fit_transform(y[seen_index]))
    else:
        kf = KFold(n_splits=n_splits)
        kfs = kf.split(X[seen_index],y[seen_index])
    i=0    
    
    def tpr(y_true, y_pred):
        return roc_curve(y_true, y_pred)[1]
    def fpr(y_true, y_pred):
        return roc_curve(y_true, y_pred)[0]
    def prec(y_true, y_pred):
        return precision_recall_curve(y_true,y_pred)[0]
    def rec(y_true, y_pred):
        return precision_recall_curve(y_true,y_pred)[1]    
    
    scores = [
        # name, function, on y when multiclas, on each y when multiclass, # proba
        ('p',precision_score, True, True, False),
        ('r',recall_score, True, True, False),
        ('f1', f1_score, True, True, False),
        ('e', accuracy_score, True, True, False),
        ('i', None, False, False, False),
        ('auc', roc_auc_score, True, True, True),
        ('tpr', tpr, False, True, True),
        ('fpr', fpr, False, True, True),
        ('prec', prec, False, True, True),
        ('rec', rec, False, True, True)
    ]    
    
    if classes:
        scores += [
            ('cov_err', coverage_error, True, False, False),
            ('LRAP', label_ranking_average_precision_score, True, False, False),
            ('LRL', label_ranking_loss, True, False, False)
        ]    
        
    for model in clf_models:
        for m in scores:
            model[m[0]] = []
        metrics = ['e']
        if classes:
            for j, y_class in enumerate(classes):
                for m in scores:
                    if m[1]:
                        model[f'{m[0]}\n{y_class}'] = []                
                        
                metrics += [f'p\n{y_class}', f'r\n{y_class}']    
                
    if test_index is not None:
        test_preds = []    
    for k_train, k_test in kfs:
        k_train = seen_index[k_train]
        k_test = seen_index[k_test]
        if test_index is not None:
            k_test = test_index
        if upsample:
            ros = RandomOverSampler(random_state=42)
            if classes:
                lp = LabelPowerset()
                yt = lp.transform(y)
                X_train, y_resampled = ros.fit_resample(X[k_train],yt[k_train])
                y_train = lp.inverse_transform(y_resampled).todense()
            else:
                X_train, y_train = ros.fit_resample(X[k_train],y[k_train].todense())
        else:
            X_train = X[k_train]
            y_train = y[k_train]
        i+=1
        print(i)
        for model in clf_models:
            clf = model['model']
            model['i'].append(i)
            #clf = SVC(kernel='rbf',probability=True)
            clf.fit(X_train,y_train)
            predictions = clf.predict(X[k_test])
            try:
                predictions_proba = clf.predict_proba(X[k_test])
            except:
                predictions_proba = predictions
                print("WARNING! Can't predict probabilities with this model, just using binary predictions")
                
            if hasattr(predictions_proba,"todense"):
                predictions_proba = predictions_proba.todense()
            if hasattr(predictions,"todense"):
                predictions = predictions.todense()
            if test_index is not None:
                test_preds.append(predictions_proba)
            if p_threshold is not None:
                predictions = np.where(predictions_proba>=p_threshold,1,0)[:,1]           
            
            if classes:
                if roundup:
                    for j, c in enumerate(predictions_proba.argmax(axis=1)):
                        predictions[j,c] = 1
                for m in scores:
                    if m[4]:
                        y_pred = predictions_proba
                    else:
                        y_pred = predictions                    
                    if not m[1] or not m[2]:
                        continue
                    try:
                        model[m[0]].append(m[1](y[k_test],y_pred,average="weighted"))
                    except TypeError:
                        model[m[0]].append(m[1](y[k_test],y_pred))
                    except ValueError:
                        pass
                    
                for j, y_class in enumerate(classes):
                    # if y[k_train,i].sum() == 0:
                    #     print("no labels for {y_class}")
                    for m in scores:
                        if not m[1]:
                            continue
                        if m[3]: # if do this metric on each class
                            if m[4]: # if use probabilities
                                y_pred = predictions_proba
                            else:
                                y_pred = predictions
                            try:
                                model[f'{m[0]}\n{y_class}'].append(m[1](y[k_test,j],y_pred[:,j]))
                            except:
                                model[f'{m[0]}\n{y_class}'].append(None)
                    if df is not None:
                        df.loc[k_test,f"{y_class} - k_prediction"] = predictions_proba[:,j]
                        df.loc[k_test,f"{y_class} - k_prediction_binary"] = predictions[:,j]
            else:
                for m in scores:
                    if not m[1]:
                        continue
                    model[m[0]].append(m[1](y[k_test],predictions))
                if df is not None:
                    df.loc[k_test, "y_k_prediction"] = predictions_proba[:,1]    
                    
    if classes:
        if df is not None:
            return clf_models, metrics, df
        return clf_models, metrics
    else:
        if df is not None:
            return clf_models, df
        elif test_index is not None:
            return clf_models, np.array(test_preds)
        return clf_models
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import KFold
from importlib import reload  
import scoping.utils.learning_utils as lu
reload(scoping.utils.learning_utils)



#NB: keep the preferred model last if you want to use the y_k prediction as it overwrites with every new model
models = [
    {
        'title':"Neural Net",
        'model': MLPClassifier(
            solver='lbfgs', 
            alpha=0.1,
            hidden_layer_sizes=(10), 
            random_state=2020
        )
    },
    {
        'title': 'Bayes', 
        'model': MultinomialNB(alpha=0.1)     
    },
    {
        'title': 'RandForest',
        'model': RandomForestClassifier(n_estimators=1000, max_depth=None,
                                        max_features="sqrt", min_samples_split=2, random_state= 2020)
    },
    {
        'title': 'SVM-rbf',
        'model': SVC(kernel='rbf',class_weight='balanced',probability=True, random_state= 2020)
    },
]

models, df = lu.cross_validate_models(ab_X, y, models, seen_index,df=df)
fig, axs = plt.subplots(1, len(models),dpi=125,figsize=(8,5),sharey=True)              
lu.plot_model_output(models,  ['p','r','e'], fig, axs)
plt.xlabel(['precision', 'recall', 'accuracy'])
plt.show()
model_df = pd.DataFrame(models)
model_df.to_csv(f"Ar&Rev_models_relevance_{now}.csv",index=False) #removed doctables/ from dir<jupyter_output><empty_output><jupyter_text>We can choose the best performing model (for our purposes) and fit it on our unseen data, in order to predict the probability that our unseen documents are relevant<jupyter_code>clf = SVC(kernel='rbf',class_weight='balanced',probability=True)

y = df['relevant']

clf.fit(ab_X[seen_index],y[seen_index])

y_pred = clf.predict_proba(ab_X[unseen_index])
plt.close()

threshold = 0.33

df['0 - relevance - prediction'] = None
df.loc[unseen_index,'0 - relevance - prediction'] = y_pred[:,1]

predicted_index = df[df['0 - relevance - prediction']>threshold].index

fig, ax = plt.subplots()

pdf = df.sort_values('0 - relevance - prediction').reset_index(drop=True)
ax.plot(pdf.index,pdf['0 - relevance - prediction'])
ax.axhline(threshold,linestyle="--",color="grey")

rel = df[pdf['0 - relevance - prediction']>threshold].shape[0]
tot = pdf[pdf['0 - relevance - prediction'].notna()].shape[0]


print(f"{rel} out of {tot} documents are predicted to be relevant")

df[['id','0 - relevance - prediction']].to_csv(f'predictions_{now}.csv',index=False)#removed doctables/ from dir

plt.show()<jupyter_output>/home/galm/software/django/tmv/venv/lib/python3.6/site-packages/pandas/core/indexes/range.py:720: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version.  Convert to a numpy array before indexing instead.
  return super().__getitem__(key)
/home/galm/software/django/tmv/venv/lib/python3.6/site-packages/matplotlib/cbook/__init__.py:1402: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version.  Convert to a numpy array before indexing instead.
  ndim = x[:, None].ndim
/home/galm/software/django/tmv/venv/lib/python3.6/site-packages/matplotlib/axes/_base.py:278: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version.  Convert to a numpy array before indexing instead.
  y = y[:, np.newaxis]<jupyter_text>Saving the outcome<jupyter_code>#Added to write relevant docs to csv
#threshold = 0.33 # set if you didn't set it above
RelevantDocs  = df[(df['0 - relevance - prediction'] >= threshold) | (df['relevant'] == 1)]
RelevantDocs[['id', 'PY', 'title', 'content']].to_csv(f'{now}_Ar&REv_PredictedRelevantSVM.csv', index = False)

RelevantDocs.shape
#writing a smaller DF with just the IDs, human-made labels, and machine prediction
outdf = df[['id', 'relevant', 'seen', '0 - relevance - prediction']]
outdf.loc[outdf['seen'] ==0, 'relevant'] = None
outdf = outdf[['id', 'relevant', '0 - relevance - prediction']]
outdf.columns = ['id', 'label', 'prediction']

outdf.to_csv(f'IdLabelPrediction_{now}.csv', index=False)
outdf.head()
#create a new query that can be used for topic modelling for example
newq, created = Query.objects.get_or_create(
    project=p,
    creator=User.objects.get(username="<EMAIL>"),
    title="The final run v2"
)

T = Doc.query.through
dqs = [T(doc_id=did, query=newq) for did in list(RelevantDocs['id'])]
T.objects.bulk_create(dqs)
newq.r_count = newq.doc_set.count()
newq.save()
# Create a tag of positvely predicted so we can screen based on those if we want
tag, created = Tag.objects.get_or_create(
    title="positively predicted".format(now),
    query=q
)
tag.save()

for i, row in RelevantDocs.iterrows(): #using the index of the newly created dataframe on the full df
    d = Doc.objects.get(pk=row['id'])
    d.tag.add(tag)
tag.update_tag()<jupyter_output><empty_output><jupyter_text>Evaluation of model<jupyter_code>#SIMPLE EXAMPLE TO GET AN IDEA FOR HOW MUCH ADDED TRAINING DATA ADDS
#Requires re-training so will take a while
from sklearn.model_selection import cross_val_score
from sklearn import metrics

#clf = SVC(kernel='rbf',class_weight='balanced',probability=True)

#y = df['relevant']

#WITH A RANDOM SAMPLE
scores = []
steps = np.linspace(0.05,1,10,endpoint=True)

#basically: take a random sample from the relevance column, use the index to create a matching sample from ab_x for different sizes
for n in steps:
    ysample = y[seen_index].sample(n=math.floor(y[seen_index].shape[0]*n))
    score = cross_val_score(clf,
                            ab_X[ysample.index],
                            ysample, 
                            cv=5,  scoring='f1')
    scores.append([len(ysample), score.mean(), score.std()*2])

scores
plt.close()
fig, ax = plt.subplots(figsize=(10,6))

xplot = [score[0] for score in scores]
yplot = [score[1] for score in scores]
yerr = [score[2] for score in scores]

ax.plot(xplot, yplot, 'k-o')

ax.fill_between(xplot,np.subtract(yplot,yerr), np.add(yplot, yerr),
                 color='gray', alpha=0.2)

ax.set_ylabel('f1 score')
ax.set_xlabel('Number of documents in sample')
ax.set_title('f1 scores for increasing sizes of test set', fontsize=14)
ax.grid()

plt.show()
#Using the original k-fold crossvalidation: Area under curve plots
from scipy import interp
fig = plt.figure(dpi=125, figsize=(8,5))
for j, model in enumerate(models):
    ax = fig.add_subplot(2,3,j+1)
    mean_fpr = np.linspace(0, 1, 100)
    mean_auc = np.mean(model[f'auc'])
    std_auc = np.std(model[f'auc'])
    tprs = []
    for i in range(len(model['i'])):
        viz = ax.plot(
            model[f'fpr'][i],
            model[f'tpr'][i],
            lw=0.5,
            alpha=0.5,
            color=f"C1"
        )
        interp_tpr = interp(mean_fpr, model[f'fpr'][i], model[f'tpr'][i])
        interp_tpr[0] = 0.0
        tprs.append(interp_tpr)
        
    mean_tpr = np.mean(tprs, axis=0)
    mean_tpr[-1] = 1.0
    ax.plot(mean_fpr, mean_tpr, color='b',
    label=r'Mean ROC (AUC = %0.2f $\pm$ %0.2f)' % (mean_auc, std_auc),
    lw=1, alpha=.8)
    ax.plot([0, 1], [0, 1], linestyle='--', lw=1, color='r',
        label='Chance', alpha=.8)

    ax.legend(fontsize=5)
    ax.set_aspect("equal")
    ax.set_title(model['title'],fontsize=8)
    ax.set_ylabel('True positive rate')
    ax.set_xlabel('False positive rate')
    
    std_tpr = np.std(tprs, axis=0)
    tprs_upper = np.minimum(mean_tpr + std_tpr, 1)
    tprs_lower = np.maximum(mean_tpr - std_tpr, 0)
    ax.fill_between(mean_fpr, tprs_lower, tprs_upper, color='grey', alpha=.2,
                    label=r'$\pm$ 1 std. dev.')
    
        
fig.tight_layout()
#Precission-recall plots 
from scipy import interp
fig = plt.figure(dpi=125, figsize=(8,5))
for j, model in enumerate(models):
    ax = fig.add_subplot(2,3,j+1)
    mean_fpr = np.linspace(0, 1, 100)
    mean_auc = np.mean(model[f'auc'])
    std_auc = np.std(model[f'auc'])
    tprs = []
    for i in range(len(model['i'])):
        viz = ax.plot(
            model[f'rec'][i],
            model[f'prec'][i],
            lw=0.5,
            alpha=0.5,
            color=f"C1"
        )
    ax.set_ylabel('precision')
    ax.set_xlabel('recall')
    ax.set_title(model['title'],fontsize=8)
        
fig.tight_layout()
from sklearn.metrics import f1_score, auc

#predict probabilities
#lr_probs = clf.predict_proba(testX)
#r_probs = lr_probs[:, 1] #keep positive outcomes

fig, ax = plt.subplots(dpi=150)

lr_probs = df.loc[seen_index,'y_k_prediction']
lr_true = df.loc[seen_index,'relevant']

# predict class values
#yhat = clf.predict(testX)

# summarize scores
lr_precision, lr_recall, _ = lr_precision, lr_recall, thresholds = precision_recall_curve(lr_true, lr_probs)

#lr_f1, lr_auc = f1_score(lr_true, lr_probs), auc(lr_true, lr_probs)
#print('For sample: f1=%.3f auc=%.3f' % (lr_f1, lr_auc))

# plot the precision-recall curves
no_skill = len(lr_true[lr_true==1]) / len(lr_true)
plt.plot([0, 1], [no_skill, no_skill], linestyle='--', label='No Skill')
plt.plot(lr_recall, lr_precision, marker='.', label='SVM')
# axis labels
plt.xlabel('Recall')
plt.ylabel('Precision')
# show the legend
plt.legend()
# show the plot
plt.show() 
plt.plot(thresholds,lr_precision[:-1])
plt.plot(thresholds,lr_recall[:-1])
#If we want to play with the thresholds, we can use this function to see what the added/removed data looks like
import seaborn as sns

def precision_recall_bins(lr_true, lr_probs, bt):    
    fig, ax = plt.subplots(dpi=150)    
    categories = []
    
    for i in range(len(bt)-1):
        categories.append((bt[i+1],bt[i]))    
        
        l = 0
    for i,c in enumerate(categories):
        y_thresh_pred = np.where(lr_probs>c[0],1,0)
        y_bin_pred = np.where(np.logical_and(lr_probs<=c[1],lr_probs>c[0]),1,0)        
        
        bin_p = precision_score(lr_true,y_bin_pred)
        thresh_p = precision_score(lr_true,y_thresh_pred)   
        thresh_r = recall_score(lr_true,y_thresh_pred)           
        width = y_bin_pred.sum()/lr_true.shape[0]        
        
        if i==0:
            rl = "Relevant"
            il = "Irrelevant"
        else:
            rl=None
            il=None
            
        ax.bar(l, 1-bin_p,bottom=bin_p, width=width, align="edge",color="#fc8d59",ec="grey",alpha=0.5,label=il)
        ax.bar(l, bin_p, width=width, align="edge",color="#91bfdb",ec="grey",alpha=0.5,label=rl)        
        
        ly = 1.15+0.2*i        
        ax.annotate(
            f"$p>{c[0]}$ \nrecall: {thresh_r:.0%}\nprecision: {thresh_p:.0%}",
            (l+width,ly),(0,ly),
            arrowprops=dict(facecolor='black', shrink=0.04,width=0.2,headwidth=4),
            va="bottom", fontsize=6
        )
        
        l += width    
        
        sns.despine()    
        
        ax.spines['left'].set_bounds(0, 1)    
        
        ax.set_ylim(ymax=ly)
    ax.set_yticks([t for t in ax.get_yticks() if t <=1])
    ax.set_ylabel("Bin composition")
    ax.set_xlabel("Cumulative proportion of documents at threshold")    
    ax.legend(loc="right")    
    fig.tight_layout()
    
precision_recall_bins(lr_true, lr_probs, [1,0.68, 0.5, 0.35, 0]) #Start with 1, end with 0 so we see whole dataset
 precision_recall_bins(lr_true, lr_probs, [1,0.35, 0])
plt.plot(thresholds,lr_recall[:-1])
from sklearn.metrics import f1_score, auc

#predict probabilities
#lr_probs = clf.predict_proba(testX)
#r_probs = lr_probs[:, 1] #keep positive outcomes

lr_probs = df.loc[seen_index,'y_k_prediction']
lr_true = df.loc[seen_index,'relevant']

# predict class values
#yhat = clf.predict(testX)

# summarize scores
lr_precision, lr_recall, _ = lr_precision, lr_recall, thresholds = precision_recall_curve(lr_true, lr_probs)

#lr_f1, lr_auc = f1_score(lr_true, lr_probs), auc(lr_true, lr_probs)
#print('For sample: f1=%.3f auc=%.3f' % (lr_f1, lr_auc))

# plot the precision-recall curves
no_skill = len(lr_true[lr_true==1]) / len(lr_true)
plt.plot([0, 1], [no_skill, no_skill], linestyle='--', label='No Skill')
plt.plot(lr_recall, lr_precision, marker='.', label='SVM')
# axis labels
plt.xlabel('Recall')
plt.ylabel('Precision')
# show the legend
plt.legend()
# show the plot
plt.show() 
plt.plot(thresholds,lr_precision[:-1])<jupyter_output><empty_output><jupyter_text>Getting false positives/-negatives for a given threshold<jupyter_code>#highest predicted
df.sort_values("y_k_prediction",ascending=False).head()
#playing with the threshold (kind of irrelevant given additional analysis above, but easy reference)
threshold = 0.36
likely = np.where(df['y_k_prediction']>threshold,1,0)
p = precision_score(df.loc[seen_index,"relevant"],likely[seen_index])
r = recall_score(df.loc[seen_index,"relevant"],likely[seen_index])

print(f"Threshold: {threshold}, precision: {p:.2g}, recall: {r:.2g}")
#print false positives just above the border
#threshold = 0.33
fp = df[(df['relevant']==0) & (df[f'y_k_prediction']>threshold)]
print(fp.shape)
django.db.connection.close()
for i, row in fp.sort_values('y_k_prediction',ascending=True).head(3).iterrows():
    d = Doc.objects.get(pk=row['id'])
    print()
    print(row['title'])
    print(row['content'])
    print(row['id'])
    print(row['y_k_prediction'])
    print(d.docownership_set.filter(query=qid).values('user__username','relevant','finish','tag__title'))
#False postives
fp = df[(df['relevant']==0) & (df[f'y_k_prediction']>0.5)]
print(fp.shape)
django.db.connection.close()

#print highest scoring ones
for i, row in fp.sort_values('y_k_prediction',ascending=False).head(3).iterrows():
    d = Doc.objects.get(pk=row['id'])
    print()
    print(row['title'])
    print(row['content'])
    print(row['id'])
    print(row['y_k_prediction'])
    print(d.docownership_set.filter(query=qid).values('user__username','relevant','finish','tag__title'))
#False negatives
fn = df[(df['relevant']==1) & (df[f'y_k_prediction']<0.5)]
print(fn.shape)
django.db.connection.close()
for i, row in fn.sort_values('y_k_prediction',ascending=True).head(3).iterrows():
    d = Doc.objects.get(pk=row['id'])
    print()
    print(row['title'])
    print(row['content'])
    print(row['id'])
    print(row['y_k_prediction'])
    print(d.docownership_set.filter(query=qid).values('user__username','relevant','finish','tag__title'))
#create a tag so I can check them out on the platform
tempq = merged_query #8510

tag, created = Tag.objects.get_or_create(
    title=f"false negatives{now}".format(now),
    query=tempq
)
tag.save()

for i, row in fn.sort_values('y_k_prediction',ascending=True).head(20).iterrows(): #using the index of the newly created dataframe on the full df
    d = Doc.objects.get(pk=row['id'])
    d.tag.add(tag)
tag.update_tag()
tag, created = Tag.objects.get_or_create(
    title=f"false positives{now}".format(now),
    query=tempq
)
tag.save()

for i, row in fp.sort_values('y_k_prediction',ascending=False).head(20).iterrows():
    d = Doc.objects.get(pk=row['id'])
    d.tag.add(tag)
tag.update_tag()<jupyter_output><empty_output><jupyter_text>Now we want to look at impact, mitigation and adaptation<jupyter_code>df['categories'] = df['1 - Mitigation'] + df['1 - Adaptation'] + df['1 - Impact']
print("seen documents:",df.loc[df['seen']==1,].shape[0])
print("seen mitigation documents:", df.loc[df['seen']==1,'1 - Mitigation'].sum())
print("seen adaptation documents", df.loc[df['seen']==1,'1 - Adaptation'].sum())
print("seen impact documents", df.loc[df['seen']==1,'1 - Impact'].sum())
print("seen in multiple", df.loc[df['categories']>=2,].shape[0])
print("seen non-mit/ad documents",df.loc[(df['seen']==1) & (df['categories']==0),].shape[0])
fig, ax = plt.subplots(dpi=150)

mitigation = df.loc[df['seen']==1,'1 - Mitigation'].sum()
adaptation = df.loc[df['seen']==1,'1 - Adaptation'].sum()
impact = df.loc[df['seen']==1,'1 - Impact'].sum()
overlap = df.loc[df['categories']>=2,].shape[0]
nonmitad = df.loc[(df['seen']==1) & (df['categories']==0),].shape[0]

lw = 1
alpha=0.5
ec = "grey"

ax.bar(
    0,mitigation,
    alpha=alpha, label="Mitigation",
    lw=lw, edgecolor=ec, tick_label = str(mitigation)
)
ax.bar(
    0,adaptation,bottom=mitigation-overlap/3,
    alpha=alpha, label="Adaptation",
    lw=lw, edgecolor=ec
)
ax.bar(
    0,impact,bottom=mitigation-2*overlap/3+adaptation ,
    alpha=alpha, label="impact",
    lw=lw, edgecolor=ec
)
ax.bar(
    0,nonmitad,bottom=mitigation-overlap+adaptation+impact,
    alpha=alpha, label="Not relevant",
    lw=lw, edgecolor=ec,
)

ax.legend()<jupyter_output><empty_output><jupyter_text>We want to do multiclass prediction nowso we name the classes and create a matrix of class labels<jupyter_code>classes = ["1 - Mitigation","1 - Adaptation", "1 - Impact"]
y = np.matrix(df[classes])
print(y.shape)

rel_index = df.loc[df["relevant"]==1,].index<jupyter_output>(286144, 3)<jupyter_text>and cross validate a model to see how well we can predict<jupyter_code>from importlib import reload  
import scoping.utils.learning_utils as lu
reload(scoping.utils.learning_utils)

from sklearn.model_selection import KFold
multimodels = [
    {
        'title': 'One vs rest SVC balanced', 
        'model': OneVsRestClassifier(SVC(kernel='linear', class_weight="balanced", probability=True))
    },
]

multimodels, metrics, df = lu.cross_validate_models(ab_X, y, multimodels, rel_index, classes=classes, df=df)

fig, axs = plt.subplots(1, len(multimodels),dpi=125,figsize=(8,5),sharey=True)              
lu.plot_model_output(multimodels, metrics, fig, [axs])
plt.show()
descriptives = [ #'e', 'p', 'r', 
                'e\n1 - Mitigation', 'p\n1 - Mitigation', 'r\n1 - Mitigation',
                'e\n1 - Adaptation', 'p\n1 - Adaptation', 'r\n1 - Adaptation',
                'r\n1 - Impact', 'p\n1 - Impact', 'r\n1 - Impact'
               ]
for des in descriptives:
    print(f'{des}  - {statistics.mean(multimodels[0][des])}')
    print(multimodels[0][des])
    print()
#With better formatting and the overall relevance included in one graph

boxplotdf = pd.DataFrame()
boxplotdf['Category'] = descriptives
gendf = pd.DataFrame({'Category' : ['e\n1 - Relevance', 'p\n1 - Relevance', 'r\n1 - Relevance']})

#get scores as list
boxplotdf['Scores'] = [multimodels[0][des] for des in boxplotdf['Category']]
gendf['Scores'] = [models[len(models)-1][des] for des in ['e', 'p', 'r']]

#merge the general and the multi-model
boxplotdf = pd.concat([gendf, boxplotdf], axis=0, ignore_index=True).reset_index()

#make a new row for each in list (basically: unfold the list to new variable, add new variable as column & drop old column)
s = boxplotdf.apply(lambda x: pd.Series(x['Scores']),axis=1).stack().reset_index(level=1, drop=True)
s.name = 'Scores'
boxplotdf = boxplotdf.drop('Scores', axis=1).join(s)

#Get the category and the metric separately
new = boxplotdf["Category"].str.split("\n1 - ", n = 1, expand = True)
boxplotdf["Metric"]= new[0]
boxplotdf.drop(columns =["Category"], inplace = True) 
boxplotdf["Category"]= new[1] 

#rename
di = {'e': "Accuracy",
      'p': "Precision",
      'r': "Recall"}
boxplotdf.replace({"Metric": di}, inplace=True)

#plot 
sns.set(style="whitegrid")
fig, ax = plt.subplots(figsize=(10,6),dpi=150)
ax = sns.boxplot(x="Category", y="Scores", hue="Metric",
                 data=boxplotdf, palette="Set2", linewidth = 2)

ax = sns.swarmplot(x="Category", y="Scores", hue="Metric", dodge=True, #add points
                   data=boxplotdf, color="grey", alpha = 0.75, size=4
                  ) 


ax.set_ylabel('Score', fontsize = 10, weight = 'bold')
ax.set_xlabel('Classifier', fontsize = 10, weight = 'bold')

title = 'Performance of Category Classifier '
subtitle = 'Accuracy, precision and recall per classifier for 10 k-fold cross-validation'
ax.text(x=0.5, y=1.1, s=title, fontsize=12, weight='bold', ha='center', va='bottom', transform=ax.transAxes)
ax.text(x=0.5, y=1.05, s=subtitle, fontsize=10, alpha=0.75, ha='center', va='bottom', transform=ax.transAxes)

handles, labels = ax.get_legend_handles_labels()
#l = plt.legend()
ax.legend(handles[0:3], labels[0:3],
          fontsize=10, title_fontsize='10', loc = 'lower left',
          facecolor = 'white', edgecolor='grey')


fig.tight_layout(rect=[0, 0, 1, 0.98])
fig.savefig(f'{title}_{subtitle}.png')

for des in descriptives: 
    print(f'{des} - - {statistics.mean(multimodels[0][des])}')
    print()
    
print ('- - - - - ')
for des in ['e', 'p', 'r']:
    print(f'{des} - - {statistics.mean(multimodels[0][des])}')
    print()
    
statistics.mean(multimodels[0]['r\n1 - Impact'])
boxplotdf = pd.DataFrame()
boxplotdf['Category'] = descriptives
gendf = pd.DataFrame({'Category' : ['e\n1 - Relevance', 'p\n1 - Relevance', 'r\n1 - Relevance']})

#get scores as list
boxplotdf['Scores'] = [multimodels[0][des] for des in boxplotdf['Category']]
gendf['Scores'] = [models[0][des] for des in ['e', 'p', 'r']]

#merge the general and the multi-model
boxplotdf = pd.concat([gendf, boxplotdf], axis=0, ignore_index = True)

#make a new row for each in list (basically: unfold the list to new variable, add new variable as column & drop old column)
s = boxplotdf.apply(lambda x: pd.Series(x['Scores']),axis=1).stack().reset_index(level=1, drop=True)
s.name = 'Scores'
boxplotdf = boxplotdf.drop('Scores', axis=1).join(s)

boxplotdf
from sklearn.model_selection import cross_val_score
cross_val_score(clf, ab_X[seen_index],y[seen_index], cv=5, scoring='recall_samples')<jupyter_output>/home/galm/software/django/tmv/venv/lib/python3.6/site-packages/sklearn/model_selection/_validation.py:552: FitFailedWarning: Estimator fit failed. The score on this train-test partition for these parameters will be set to nan. Details: 
Traceback (most recent call last):
  File "/home/galm/software/django/tmv/venv/lib/python3.6/site-packages/sklearn/model_selection/_validation.py", line 531, in _fit_and_score
    estimator.fit(X_train, y_train, **fit_params)
  File "/home/galm/software/django/tmv/venv/lib/python3.6/site-packages/sklearn/svm/_base.py", line 162, in fit
    accept_large_sparse=False)
  File "/home/galm/software/django/tmv/venv/lib/python3.6/site-packages/sklearn/base.py", line 432, in _validate_data
    X, y = check_X_y(X, y, **check_params)
  File "/home/galm/software/django/tmv/venv/lib/python3.6/site-packages/sklearn/utils/validation.py", line 73, in inner_f
    return f(**kwargs)
  File "/home/galm/software/django/tmv/venv/lib/python3.6/site-packages/sklearn/utils/v[...]<jupyter_text>Multi-matrix<jupyter_code>import seaborn as sns
sns.reset_defaults()
classes = ["1 - Mitigation","1 - Adaptation", "1 - Impact"]
cols = [f"{c} - k_prediction_binary" for c in classes]
print(cols)
y_pred = np.array(list(df.loc[seen_index,cols].values))
   
y_true = y[seen_index]

confusion_matrix = np.zeros(shape=(len(classes)+1,len(classes)+1),dtype=float)

norm = True

for j, c in enumerate(classes): #rows
    for k, c in enumerate(classes): #cols
        n = 0
        for i in range(len(seen_index)):
            if j==k:
                if y_true[i,j]+y_pred[i,k]==2:
                    n+=1
            else:
                if y_pred[i,j]+y_true[i,k]==2 and y_pred[i,k]+y_true[i,k]!=2:
                    n+=1/y_pred[i,].sum()
        confusion_matrix[j,k]=n
        
    sums_pred = y_pred.sum(axis=1).ravel()
    confusion_matrix[k+1,j]=y_true[np.argwhere(sums_pred==0).ravel(),j].sum()
    sums_true = y_true.sum(axis=1).ravel()
    confusion_matrix[j,k+1]=y_pred[np.argwhere(sums_true==0).ravel(),j].sum()
    confusion_matrix[j+1,k+1] = y_pred[np.argwhere(sums_true+sums_pred==0)].shape[0]


if norm is not False:
    confusion_matrix = confusion_matrix/confusion_matrix.sum(axis=0,keepdims=True)


fig, ax = plt.subplots(dpi=125, figsize=(8,8))


if norm is not False:
    mat = ax.imshow(confusion_matrix, cmap="Blues", vmin=0, vmax=1, origin="lower")
else:
    mat = ax.imshow(confusion_matrix, cmap="Blues",origin="lower")
    
classes = classes + ["None"]
for j, c in enumerate(classes):
    for k, c in enumerate(classes):
        n  = confusion_matrix[j,k]
        if n > np.max(confusion_matrix)/2:
            color="#d9d9d9"
        else:
            color="#949680"
        if norm is not False:
            ax.text(k,j, f"{confusion_matrix[j,k]:.0%}",ha="center",color=color)
        else:
            ax.text(k,j, f"{confusion_matrix[j,k]:.1f}",ha="center",color=color)

    if j < len(classes)-1:    
        ax.text(j,len(classes)-0.25,f"{y_true[:,j].sum():,.0f}",ha="center",va="center")
        ax.text(len(classes)-0.25, j, f"{y_pred[:,j].sum():,.0f}",ha="center",va="center")
        
ax.set_xticks(list(range(len(classes))))
ax.set_xticklabels([x.replace("12 - ","") for x in classes],rotation=45,ha="right")

ax.set_yticks(list(range(len(classes))))
ax.set_yticklabels([x.replace("12 - ","") for x in classes],ha="right")
ax.set_xlabel("True value", fontsize = 10, weight = 'bold')
ax.set_ylabel("Predicted value", fontsize = 10, weight = 'bold')

#ax.set_xlim(xmax=len(classes)+1)

cbar = fig.colorbar(mat, ax=ax, pad=0.15, shrink = 0.5)
if norm is not False:
    cbar.set_label("% of True values")
    
title = "Confusion matrix for predicted categories"
ax.text(x=0.5, y=1.1, s=title, fontsize=12, weight='bold', ha='center', va='bottom', transform=ax.transAxes)

fig.tight_layout(rect=[0, 0, 1, 0.98])
fig.savefig(f'confusion_matrix.png')
fig.show
model_df = pd.DataFrame(multimodels)
#model_df.to_csv(f"multimodels_adaptation_mitigation_{now}.csv",index=False)<jupyter_output><empty_output><jupyter_text>Fitting multi-class<jupyter_code>predicted_index = df[df['0 - relevance - prediction']>0.33].index
LabeledRelevant_index = df[df['relevant']==1].index

clf = OneVsRestClassifier(SVC(kernel='linear', class_weight="balanced", probability=True))
classes = ["1 - Mitigation","1 - Adaptation", "1 - Impact"]
y = np.matrix(df[classes])
y.shape

clf.fit(ab_X[seen_index],y[seen_index]) #fit on the labeled as relevant data
predictions = clf.predict_proba(ab_X[predicted_index]) #get probability scores for those predicted to be relevant
binary = clf.predict(ab_X[predicted_index]) #get binary scores too
#Add predictions to the dataframe or nan if no prediction
for i,c in enumerate(classes):
    df[f"{c} - prediction"] = np.NaN
    
    df.loc[predicted_index, f"{c} - prediction"] = predictions[:, i]
#Check
print(df[df['1 - Adaptation - prediction'] > .1].shape)
df[df['1 - Adaptation - prediction'] > .1].head(20)<jupyter_output>(3281, 23)<jupyter_text>False positives & NegativesFirst, let's look into the high number of adaptation docs being classified as impact<jupyter_code>#Just because I got tired of getting the titles confused, import markdown display
from IPython.display import Markdown, display
def printmd(string): 
    display(Markdown(string))
    
confusion = df[
   (df['1 - Adaptation'] == 1) 
   & (df['1 - Adaptation - k_prediction_binary'] == 0) 
   & (df['1 - Impact - k_prediction_binary'] == 1)
    ]

for i, row in confusion.sort_values('1 - Impact - k_prediction',ascending=False).head(5).iterrows():
    d = Doc.objects.get(pk=row['id'])
    print()
    printmd(f"**{row['title']}**")
    print(row['content'])
    #print(row['id'])
    printmd(f"*Impact prediction: {row['1 - Impact - k_prediction']} -- Adaptation prediction: {row['1 - Adaptation - k_prediction']}*" )
    #print(d.docownership_set.filter(query=qid).values('user__username','relevant','finish','tag__title'))
#Listing all false positives and false negatives
#NB: this includes docs belonging to multiple classes where only one was classified wrongly
classes_pred = [(c, c+" - k_prediction_binary") for c in classes[0:-1]]

c_fp = pd.DataFrame(columns = df.columns)
c_fn = pd.DataFrame(columns = df.columns)

for c, c_p in classes_pred:
    c_fn = c_fn.append(df[(df[c]==1) & (df[c_p]==0)])
    c_fp = c_fp.append(df[(df[c]==0) & (df[c_p]==1)])
    
print(f"False negatives: {c_fn.shape[0]} \nFalse positives: {c_fp.shape[0]}")
df['categories'] = df['1 - Mitigation'] + df['1 - Adaptation'] + df['1 - Impact']
print(f"Total categorised: {df.loc[df['categories']>=1].shape[0]}")
#create a tag so I can check them out on the platform
tempq = merged_query

tag, created = Tag.objects.get_or_create(
    title=f"Categories f neg {now}".format(now),
    query=tempq
)
tag.save()

for i, row in c_fn.iterrows(): #using the index of the newly created dataframe on the full df
    d = Doc.objects.get(pk=row['id'])
    d.tag.add(tag)
tag.update_tag()
#same for false positives
tag, created = Tag.objects.get_or_create(
    title=f"Categories f pos {now}".format(now),
    query=tempq
)
tag.save()

for i, row in c_fp.iterrows(): #using the index of the newly created dataframe on the full df
    d = Doc.objects.get(pk=row['id'])
    d.tag.add(tag)
tag.update_tag()<jupyter_output><empty_output><jupyter_text>To CSV<jupyter_code>columns = ['1 - Adaptation', '1 - Mitigation', '1 - Impact' ]
ycols = [col + ' - k_prediction' for col in columns]
pcols = [col + ' - prediction' for col in columns]
outdf = df[['id', 'relevant', 'seen']]
outdf[ycols] = df[ycols]
outdf[pcols] = df[pcols]
outdf.head()
#writing a smaller DF with just the IDs, human-made labels, and machine prediction
#outdf = df[['id', 'relevant', 'seen', '0 - relevance - prediction']]
outdf.loc[outdf['seen'] ==0, 'relevant'] = None
for n, ycol in enumerate(ycols):
    outdf[ycol] = np.where(outdf['seen'] == 1, outdf[ycol], outdf[pcols[n]])
       
outdf = outdf.rename(columns={
    'relevant': 'label',
    '1 - Impact - k_prediction': 'impact',
    '1 - Adaptation - k_prediction': 'adaptation',
    '1 - Mitigation - k_prediction': 'mitigation'
})

outdf = outdf[['id', 'label', 'adaptation', 'mitigation', 'impact']].dropna(thresh=5)

outdf.to_csv(f'IdCategoryPrediction.csv', index=False)
outdf.head()
outdf.shape
outdf.to_csv(f'IdCategoryPrediction.csv', index=False)
outdf.shape
outdf[outdf['mitigation'] == 1].shape<jupyter_output><empty_output>