# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Overview
#
# This notebook contains all code and discussions for the __intermediate iteration__ of the research question involving __lesion coloration and malignancy__. The specific question is whether, to a statistically significant degree, coloration (*i.e.*, darkness or intensity) can be used to predict malignancy (*i.e.*, occurrence of malignant diagnosis). At this level, the question involves ANOVA and logistic regression techniques. It is an advancement of [the novice module](../novice/Q1.ipynb) on this topic.
#
# # Table of Content
#
# I. [Setup](#s)
#
# II. [Data Retrieval](#dr)
#
# 1. [File Access](#fa)
# 2. [Loading & Processing](#lp)
#
# III. [Analysis](#anl)
#
# 1. [ANOVA](#anova)
# 2. [Logistic Regression](#lr)
#
#
# IV. [Discussion](#discuss)

# # Setup <a id='s'></a>
#
# Refer to [this module](../novice/Data_Download.ipynb) for replicating the data-download process using the ISIC Python tool. The command is repeated below, although file/directory locations may need adjustment. 
#
# `
# python download_archive.py \
# --num-images=50 \ 
# --images-dir ../../sample_imgs \ 
# --descs-dir ../../sample_dscs -s \ 
# --seg-dir ../../sample_segs --seg-skill expert
# `
#
# The following are necessary imports for this module.

# +
# file access
import json
import glob 

# image handling
from PIL import Image

# data manipulation and visualization
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

# analysis
from scipy.stats import f_oneway
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import train_test_split
from sklearn.model_selection import GridSearchCV
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import f1_score
from sklearn.metrics import plot_roc_curve
# -

# # Data Retrieval <a id='dr'></a>
#
# ## File Access <a id='fa'></a>
#
# As seen in novice modules, all three pieces of data are needed in this question: images, segmentations, and descriptions. Load as before.

# +
img_filepaths = glob.glob('../../sample_imgs/*.jp*')
seg_filepaths = glob.glob('../../sample_segs/*.png')
dsc_filepaths = glob.glob('../../sample_dscs/*')

print(len(img_filepaths), 'images')
print(len(seg_filepaths), 'segmentation files')
print(len(dsc_filepaths), 'descriptions')
# -

# Loop through the image filepaths to make sure the identifying numbers in every $k^{th}$ image file matches the number in the corresponding $k^{th}$ segmentation and description files. That is, check that the corresponding elements of the three lists refer to the same data point.

for i in range(len(img_filepaths)):
    img_num = int(img_filepaths[i].split('_')[-1].split('.')[0])
    seg_num = int(seg_filepaths[i].split('_')[-2])
    dsc_num = int(dsc_filepaths[i].split('_')[-1])   
    wrong_dsc = (dsc_num != img_num) or (dsc_num != seg_num)
    wrong_seg = seg_num != img_num  
    if wrong_dsc or wrong_seg:
        print("Error in file order")
        break
print('Validated.')

# ## Loading & Processing <a id='lp'></a>
#
# There's only one required part of the description files, which is the benign/malignant indicator. Descriptions (JSON format) have this information nested in the `clinical`-`meta`-`benign_malignant` key.

json.load(open(dsc_filepaths[0], 'r'))

# Create a new list where each $k^{th}$ description is represented by the $k^{th}$ malignancy value.

# +
dscs = [json.load(open(x, 'r'))['meta']['clinical'] for x in dsc_filepaths]
flags = [x['benign_malignant'] for x in dscs]

# sample output
for i in range(3):
    print('Data point #%d: %s' % (i, flags[i]))
# -

# The image files need to be opened and made into matrices. The matrices have the following structure:
#
# - *m* rows
#     - *n* columns per row
#         - *3* pixel values (RGB order or red, green, blue pixel values) per column 
#         
# All images can be opened according to this pattern.

images = [Image.open(x).resize((400,300)) for x in img_filepaths]
images = [np.array(x.convert('RGB')) for x in images]

# Segmentation files, as seen in previous modules, can be opened without RGB notation. Their structure is:
#  
# - *m* rows
#     - *n* columns
#         - *1* value (0 for benign, 1 for malignant) per column (not an array of 3 like in the images)

segs = [Image.open(x).resize((400,300)) for x in seg_filepaths]
segs = [np.array(x) for x in segs]

# The next step is to use the segmentations to extract only the lesion part of the images. Notice how the segmentation creates a boundary; coordinates in the dark area are not relevant.

fig, axs = plt.subplots(1,2)
axs[0].imshow(images[0])
axs[1].imshow(segs[0])
plt.show()

# The above images give an idea of how to extract the meaningful portion. Although the segmentation appears colored, it is actually black and white. At any location *m,n*, a 0 indicates a black (non-relevant pixel), and non-0 value indicates lesion area. Test by locating points that are visibly inside the lesion area.

# +
# naked-eye test for location in the lesion area at row 100, column 300
x1, s1 = images[0][100,300], segs[0][100,300]
print('At row 100, column 300 the segmentation value is non-black/non-0:', s1)

# naked-eye test for location outside at row 10, column 50
x1, s1 = images[0][10,50], segs[0][10,50]
print('At row 10, column 50 the segmentation value is black/0:', s1)
# -

# If a black pixel occurs within the lesion, it is part of the lesion's coloring, and should be counted (to distinguish between very dark and very very dark cells, for example). With this in mind, use indexing to extract the lesion-only portion for a sample image.

test_img = images[0]
test_seg = segs[0]

# Find all non-0 (non-black) occurrences in the segmentation. WHAT DOES THE RESULT MEAN

idxs = np.where(test_seg)
idxs

# Use these locations to extract only the lesion portions of the image.

test_rdx = test_img[idxs]
for i in range(3):
    rgb = test_rdx[i]
    print('RGB at location #%d where segmentation is non-0: %s' % (i, rgb))

# Apply this filtering to all images. Each element of the resulting list represents one image, and consists of a non-constant amount of triplets representing the relevant RGB pixels for that image. 

image_rdxs = []
for i in range(len(segs)):
    idxs = np.where(segs[i]) # clarify that default is checking 0 value
    image_rdx = images[i][idxs] # test code that shows how np.where successfully tells you the right locations
    image_rdxs.append(image_rdx)

# # Analysis <a id='anl'></a>
#
# ## ANOVA <a id='anova'></a>
#
# If the above list represents image "subsets", ANOVA can be used to identify significant differences in color among these subsets, given some grouping thereof. First, the subsets have to be boiled down from lists of RGB values to one value. See the first data point:

image_rdxs[0]

# In each RGB triplet, the first value is the red pixel value, the second the green value, and lastly the blue value. Access the mean red, green, and blue values across all these triplets. 

R, G, B = np.mean(image_rdxs[0], axis=0)
print('Mean R: %0.1f, mean G: %0.1f, mean B: %0.1f' % (R,G,B))

# This can be turned into a single, hexadecimal number (how colors are commonly represented), yielding one numeric summary of the averaged RGB-combination in this image.

hex_value = '%02x%02x%02x' % (int(R), int(G), int(B))
print('The hex representation is', hex_value)

# Apply to all images.

# +
rgb_means = [np.mean(x, axis=0) for x in image_rdxs]
rgb_means = [(int(r), int(g), int(b)) for (r,g,b) in rgb_means]
hex_values = ['%02x%02x%02x' % (r,g,b) for (r,g,b) in rgb_means]

# sample output
for i in range(3):
    print('Image #%d as hex: %s' % (i, hex_values[i]))
# -

# Now apply a grouping to these colors. Turn the hex values into their integer counterparts, and sort the images from light to dark (*e.g.*, lowest int values to darkest).

hex_to_int = [int(x,16) for x in hex_values]
sorted_ints = sorted(hex_to_int)

# An example grouping is given: the lightest four images are the 'light' group, the next four are the 'medium' group, the remaining darkest colors are the 'dark' group.

light, med, dark = sorted_ints[:4], sorted_ints[4:8], sorted_ints[8:]

# Performing one-way ANOVA (no response) gives a test statistic and p-value. If p < 0.05 (the common standard), reject the null hypothesis (no significant difference among groups) and argue in favor of the alternative hypothesis.
#
# __*Note:*__ the alternative hypothesis in this case argues that there is significant difference in __color__ among these images. It says that there are non-trivial differences in color among the groups, and so the data is not an essentially identical set of images.

f_oneway(light, med, dark)

# ## Logistic Regression <a id='lr'></a>
#
# To predict malignancy based on color, notice the binary response (benign/malignant) and multiple predictors (mean R, mean G, mean B). Logistic Regression is perfect for such a binary classification problem. The hexadecimal operations aren't even needed, since each color channel can be entered as a separate attribute. 
#
# Begin with the response variables. They are still in the order in which they were loaded from files, and can be encoded into binary numeric outcomes. Benign is 0, malignant is 1. 

le = LabelEncoder()
y = le.fit_transform(flags)
y.shape

# Mean R,G,B values have already been obtained for each image. Refer to the step prior to sorting color values.

for i in range(3):
    print('RGB for Image %d is %s' % (i, rgb_means[i]))

# Transform into a numeric array, and then create splits of the data for training, validating, and testing the model.

# +
X = np.array(rgb_means)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1000)

# verify results
print(X_train.shape, X_test.shape)
print(y_train.shape, y_test.shape)
# -

# Set up Logistic Regression using the F1 metric for the binary classification problem (this is a weighted mean of precision and recall, and relevant to this task). Create a 'grid' of hyper-parameters for the internal Logistic Regression training process, so that the best combinations thereof can be found. Finally, train on training data and show progress.

# +
# l1 refers to lasso regression, l2 to ridge regression
hyperparameters = {'C': np.logspace(-3,2,6), 'penalty':['l1', 'l2']} 

# default 5-fold CV is used
clf = GridSearchCV(LogisticRegression(solver='liblinear'), hyperparameters, scoring='f1') # is scoring by f1 wrong?
clf.fit(X_train, y_train)
print('Optimized parameters', clf.best_params_)

lr = LogisticRegression(solver='liblinear', penalty='l1', C=0.001, max_iter=200)
lr.fit(X_train, y_train)
# -

# Output F1-based performance from the model, on test data to determine generalizability. This value should be close to 1 (the maximum) and far from 0. 

preds_test = lr.predict(X_test)
f1 = f1_score(y_test, preds_test) # try again
print('F1:', f1)

print(y_test)
print(preds_test)
print(np.mean(y_test == preds_test))
# calc prec/recall yourself and see if different

# For a more visual result, consider the ROC plot, which should maximize the Area Under Curve (AUC).

plot_roc_curve(clf, X_test, y_test)

# # Discussion <a id='discuss'></a>
#
# *__ANOVA__*
#
# *Discuss the findings revealed by ANOVA on the manual groupings of images by coloration.*
#
# *__Logistic Regression__*
#
# *Discuss the impact and accuracy of Logistic Regression on malignancy using colors as predictors.*

StarcoderdataJupyterScriptsDedupFiltered