Fluorescence spectroscopy and machine learning methods for detection of Alzheimer's disease from circulating white blood cells.

Code ↔ Paper

4 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 4 matches

[1] § Methods › Machine learning methods and statistics for spectral image analysis ↔ stys_WTv5xFAD.py, lines 46–86 · score 0.82 · dimensional space, wavelet transform, stronger, weight, Euclidean, selection
[2] § Methods › Machine learning methods and statistics for spectral image analysis ↔ stys_WTv5xFAD.py, lines 88–128 · score 0.75 · cross validated, decision boundary, Mac, rbf, minimizing, gamma
[3] § Methods › Machine learning methods and statistics for spectral image analysis ↔ stys_WTv5xFAD.py, lines 1–44 · score 0.69 · surface scan, variance, UMAP1, MANOVA, SVM, subsets
[4] § Results › Machine learning-based analysis of human PBMC immunoprecipitates ↔ stys_WTv5xFAD.py, lines 1027–1091 · score 0.59 · support vector regressor, cross validation, hyperparameters, splits, fitting, score

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 4,178 lines · 195 KB · CC-BY-4.0 · 4 matches

"""
Python script used for wavelet decomposition>UMAP dim reduction>SVM classification_report
Title: Fluorescence spectroscopy and machine learning methods for detection of Alzheimer's disease from circulating white blood cells
Authors: Shigeki Tsutsui1, Anastasiia A. Stepanchuk, Julian P. Stys, Stefanie A.G. Black, George W. Templeton, Russell Greiner, Peter K. Stys
Correspondence:
Peter K. Stys, MD
Hotchkiss Brain Institute, Department of Clinical Neurosciences, Cumming School of Medicine, University of Calgary,
3330 Hospital Drive, NW
Calgary, AB Canada T2N 4N1
Email: [email hidden]
rev.2025-12-23
"""
nShuffles = 10 # set to >0 to perform a final RAND sanity check (typically 5-10)
dimReductionMethod = 'UMAP'
nComponents = 2 # best total # of components for the dim reduction by SFS (must be >1 and <=maxComponents)
maxComponents = 3 # number of components to be calculated in the reduction, from which the top nComponents are selected. Increasing maxComponents may yield better results but interpret with caution if the PCA % variance is low. If you want to default to the old behavior of selecting the FIRST nComponents set maxComponents=nComponents
CWT_type = ['fbsp1-0.40-0.048'] # you can enter an explicit list of wavelet names here, otherwise 'deep' tells the script to scan a standard set for the best one
CWT_N_scales = 32 # no. of scales = pixel height of scalogram (CAUTION: increasing this value excessively will place greater demand on memory, possibly crashing the system)
scalogramCompressionPowerFactors_start = -0.2 # setting both scalogramCompressionPowerFactors_start & _end to negative values generates a series of only negative compression levels (sign of scalogramCompressionPowerFactors_nsteps is ignored)
scalogramCompressionPowerFactors_end = -0.2
### SETTINGS ###
# *** IMPORTANT: set CWT_N_scales=32 below if you change ExportKernelwiseReductions=True else you will almost certainly crash your system because of out-of-memory
ExportKernelwiseReductions = False # set to False if you only want dim reductions on subject-wise average spectra, otherwise True will compute reductions on each kernel for surface scan import but will take a long time and is very memory and compute intensive. Note: not all wavelets will generate discernible clusters on the imported data surfaces, consider starting with sourceData=0 below.
ShuffleMode = 0 # 0: only the classIDs are shuffled: some names will end up in the other class; 1: shuffle lambda rows: names retain their class IDs but inherit someone else's spectrum
ClassPair = [0,1] # the 2 classIDs to compare (CSVs may have more than class 0 & 1; must have at least 2 elements)
kNormMode = 1 # 0: don't normalize; 1: normalize each spectrum to its own peak
StandardScaleReductions = True # if true, stdscales dim-reduced arrays after pca, umap, etc
LambdaInputRange = [-1, 800] # restrict wavelength range [start nm,end nm] (spectral data outside this range will be deleted after load, but BEFORE normalization). Set start = -1 to include all wavelengths (NOTE: scalogramLambdaStartEndRange determines the actual range used for analyses; LambdaInputRange is mainly for display purposes)
FirstLambdaCol = -1 # 0-based index of 1st wavelength column; typically 3; use -1 for auto-detect from column headers
nLambdas = 32 # no. of wavelength columns in input csv; if FirstLambdaCol = -1 nLambdas will be auto-detected
SubSampleSize = int(0) # 0 to use all rows in csvIN, otherwise set # of rows to use as a subset. Approx. 31GB of memory required per 1M input kernels
ReplaceNameUnderscoresWithDashes = True # set to True so that names with underscores are not misinterpreted as replicates of same subject (_ char typically represents replicates of same subject)
LimitMemoryUse = True # applies only when ExportKernelwiseReductions = True: set to False to bypass this script's attempt at limiting SubSampleSize to avoid OOM crashes (in which case watch memory pressure very carefully!)
T = [1] # if csv has a 'T' column, retain only these T points (1-based). Set to -1 to retain all T-points. Can do eg T = [1,10] to retain an inclusive range of T points. Ignored if CSV has no T column.
D = 0 # 0: use raw 0-order spectra; >0: higher order spectral derivatives
# Savitzky–Golay filter params (set w_size = 0 for no filtering, applied to average spectra only, mainly applicable to Mk6 and Duetta spectra). You can examine what the filters will do using the CSV summary script first.
SG_w_size = 0 # window size (must be odd); 0: no filtering; otherwise a value of 17 is a good choice for Duetta spectra; 31-41 for Mk5
SG_polyorder = 3 # polynomial order
componentList = [] # explicit 1-based list of components to use, overrides nComponents option. IMPORTANT: you may need to pass an explicit maxComponents (> greatest component in componentList) because UMAP (not sure about PCA) will generate different results for eg UMAP1 & 3 if maxComponents=3 vs 4 for instance. Settings maxComponents=0 will use the largest value in componentList as maxComponents.
tie_breaker = 'dbi' # manova, silhouette (not as good?: permutation, eta_squared, mahalanobis, linsep, auc, dbi [Davies-Bouldin Index])
# CWT OPTIONS (these only apply if sourceData = 1 below):
CWT_scales_start = 0.3 # higher frequencies
CWT_scales_end = 150 # lower frequencies
CWT_scales_power = 2 # 1: linear, >1 more and more convex upward (see def custom_growth); controls the spread of scales/frequencies along the y-axis of the scalogram, higher powers give more grain at higher frequencies (smaller scales)
CWT_generateStdWaveletList_n_linspace_points = 6 # only with CWT_type = 'deep' & 'fbsp-xxx'; see CustomWavelet.generateStdWaveletList (results can be very sensitive to small adjustments in wavelet type)
# log sequence of power factors: smaller scalogramCompressionPowerFactors (more compression) will put more weight on lower amplitude components of your spectrum, which may or may not aid with group separation
scalogramCompressionPowerFactors_nsteps = -7 # setting nsteps to a negative value will additionally append an equivalent set of -ve compression levels e.g. -4 will generate a series of 8 compression levels, 4 negative, 4 positive
scalogramLambdaStartEndRange = [400, 800] # enter a restricted lambda range (in nm) for analysis/scalogram calculation. Not equivalent to LambdaInputRange if normalization is used.
# scalogram_wsd_threshold =
# None: omit WSD
# 0.01-0.99: absolute normalized scalogram difference threshold: Lower threshold: More inclusive, keeps more regions of the scalogram; higher threshold (e.g., 0.6): more selective, focuses only on regions with the strongest differences
# >= 1 interpreted as the percentage strongest total difference scalogram pixels to include
# -1: auto-determine threshold
scalogram_wsd_threshold = -40 # None: omit WSD feature engineering; -1 auto-determine threshold; 0.01-0.99 used as an absolute threshold value on the normalized scalogram difference; >= 1 interpreted as the percentage strongest total difference scalogram pixels to include; <-1 abs number of strongest pixels in scalogram diff to use (let's you set the abs # of extracted "features" from the scalograms) [see def CV_nested_cached]
scalogram_wsd_window_percent = 20 # width/height of sliding window when computing mean scalogram diffs, as a % of nLambdas & CWT_N_scales
scalogram_wsd_mask_power = 0 # scalogram mask power (0 for binary mask; > 0 graded mask values raised to this power e.g. 1 for untransformed graded mask; see class WaveletTransformer_cached_WSD)
# UMAP options
umap_repeats = 5 # UMAP appears to return highly variable solutions, many more repeats may be necessary to find even better solutions (NOTE: UMAP library has a memory leak so too many repeats will run out of mem)
umap_n_neighbors = 30 # Higher values (e.g., 30–100) emphasize inter-class separation by capturing broader patterns in the data. Limited to nSamples.
umap_min_dist = 1.0 # Higher values (0.5–1.0) spread clusters apart, enhancing inter-class separation. Use values > 1 with caution (check silhouette and DBI cluster metrics)
umap_metric = 'cosine' # Defines the distance metric for high-dimensional space (e.g., euclidean, cosine). Non-Euclidean metrics (e.g., cosine for text) can better separate classes with non-linear relationships (https://umap-learn.readthedocs.io/en/latest/api.html)
umap_init = 'spectral' # 'spectral','random','pca','tswspectral' (https://umap-learn.readthedocs.io/en/latest/api.html)
# CNN_AE options (only when dimReductionMethod = 'CNN_AE')
cnn_hidden_dims = [32] # determines no. of Conv2/3D layers (= no. of elements in the list). Add additional elements to this list to add additional Conv2/3D layers
cnn_epochs=100
cnn_batch_size=32
cnn_use_attention = False # if True sets scalogram_wsd_threshold = -1 to omit WSD mechanism, can't use both (develop attention mechanism later, seems to run but needs optimization)
# CNN_AE_PCA options (only when dimReductionMethod = 'CNN_AE_PCA')
cnn_latent_dims = 64 # for CNN_AE_PCA only: latent dim of AE portion of the pipeline, which are then further reduced to nComponents,maxComponents by a final PCA
# CNN1D_AE_PCA options (only when dimReductionMethod = 'CNN1D_AE_PCA' and sourceData = 0)
cnn1d_pca_hidden_dims = [16]
cnn1d_pca_latent_dims = 8 # latent dim of AE stage of the pipeline, which is then further reduced to nComponents,maxComponents by a final PCA. may need to increase this to 32 or 64 for Mk6/Duetta spectra
cnn1d_pca_epochs = 500
cnn1d_pca_batch_size=32
cnn1d_pca_kernel_size=3
cnn1d_pca_stride=2
# LSTM_AE options
lstm_bidirectional = True
lstm_hidden_layers = 64
lstm_epochs = 200
lstm_batch_size = 16 # also for CNN_LSTM_AE
lstm_learning_rate=1e-3 # or None
lstm_patience=10
lstm_min_delta=2e-4 # early stopping delta
# CNN_LSTM_AE & CNN_LSTM_AE_PCA options
cnn_lstm_conv_sizes = [32] # determines no. of Conv2 layers (= no. of elements in the list). Add additional elements to this list to add additional Conv2 layers
cnn_lstm_hidden_size = 64 # number of units/neurons in the LSTM layer:
cnn_lstm_latent_dims = 64 # for CNN_LSTM_AE_PCA only: latent dims of AE portion of the pipeline, which are then further reduced to nComponents,maxComponents by a final PCA reduction
cnn_lstm_epochs = 200
cnn_lstm_batch_size = 16
# SVM/CV options
import numpy as np
svm_param_grid = {'C': np.logspace(-2, 2, 6), 'gamma': ['scale'], 'kernel': ['linear'], 'degree': [3]} # expanding C np.logspace, with larger range and more grain may find better solutions; gamma='scale' lets the SVM decide. Only a single entry for 'kernel' is allowed eg: linear, rbf (tends to heavily overfit for small sample sizes), sigmoid, poly. degree applies to 'poly' kernels only
svm_MarginScorer_margin_weight = 0.01 # margin weight for SVM scoring function, see MarginScorer class. Higher values will sacrifice acc too much in favor of wider SVM margins. Small values will ensures that acc is first. Set to 0 to use 'accuracy' scoring only.
svm_distance_mode = 1 # 0: distances are computed using svc.decision_function (this mode is forced for rbf kernels); 1: distances are computed manually using decision boundary (only applies to linear kernels)
svm_iters = 5000 # max iters of base SVM in __compute_separation_sfs_xxx methods
svm_random_state = 42
k_q = 15 # threshold for LOOCV (N<k_q) vs RepeatedStratifiedKFold cross-validation (see CV_nested function)
k_val_acc_tiebreaker_tolerance = 0.9 # see CV_nested method. Set to 1.0 to always use max val accuracy for tiebreaker pipeline set
# auto-clustering options:
nClusters = 2 # number of expected clusters in dataset (only applies to Clustering acc result)
# plotting/output options:
LabelDatapoints = False # label datapoints on the graphs with subject IDs
subjectEmphasisList = [''] # enter a list of subject names to plot their spectra and markers in bold for identification e.g. subjectEmphasisList = ['HC-DN318','AD-BG0911']
saveSVGgraphs = True # set to True to also export svg version of all graphs for import into graphics apps and editing for publication figs
saveModel = True # set to True to also export svg version of all graphs for import into graphics apps and editing for publication figs
# DEBUG
kUseAUGforAllMethods = False # experimental: uses augmented datasets for all reducers, not just NN-based that require train-val splits. Currently only PCA is supported (in addition to NN-based reducers)
kVerbose = False
kMP = True # multiprocessing switch
kGPU = True # use GPU if available (on Silicon Macs choose Stats Buddy>Hide Others and minimize all Stats Buddy windows to maximize GPU availability)
#################
# full path to input CSV (leave empty for SB to replace from GUI @ script launch; path cannot contain single quotes):
csvIN='/Volumes/LabStuff/ITKs & images/test CSVs/AS_WT v 5XFAD_hi gain_PLAQUES OMITTED_K5.csv'
# full path to input directory enclosing 1 or more CSVs (takes precedence over csvIN if both are defined). NOTE: Finder aliases will NOT be resolved:
dirIN=''
# full path to output directory (leave empty for SB to replace from GUI @ script launch; path cannot contain single quotes):
dirOUT='/Users/pstys/Documents/dirOUT/'
# full path to output file, if defined all strings sent to printSB will also be written to this file (useful for ARC jobs to write interim results to a file that don't seem to be printed in stdOut until the end of the job):
# printSB_fOut=''
import sys
import time
import math
from scipy import stats
import pandas as pd
import os
import sklearn
import matplotlib.pyplot as plt
from pathlib import Path
import seaborn as sns
import pywt
from datetime import datetime
import random
import operator
import multiprocessing
from multiprocessing import Manager
from scipy.signal import correlate
import psutil
from itertools import combinations
from collections import defaultdict
from sklearn import metrics, decomposition
from sklearn.metrics import mean_squared_error, r2_score, make_scorer
from sklearn.model_selection import permutation_test_score
# from sklearn.base import BaseEstimator, RegressorMixin, TransformerMixin
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import MiniBatchKMeans, KMeans, SpectralClustering
from sklearn.metrics import classification_report, confusion_matrix, accuracy_score
from sklearn.decomposition import NMF
import umap
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers, models, Model
########################### STANDARD_FUNCTION_BLOCK ###########################
import math
import warnings
import os
import numpy as np
import pandas as pd
from datetime import datetime
from statsmodels.multivariate.manova import MANOVA
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import GridSearchCV, train_test_split, LeaveOneOut, cross_val_score, StratifiedKFold, RepeatedStratifiedKFold, cross_val_predict
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_squared_error, accuracy_score, silhouette_score, roc_auc_score, roc_curve, auc, davies_bouldin_score
from sklearn.ensemble import RandomForestClassifier
import random
import re
import tensorflow as tf
from tensorflow.keras import layers, models
from scipy.optimize import minimize
import seaborn as sns
warnings.filterwarnings("ignore") # suppress all warnings for this script run
def formatPstring(P, includePrefix=True):
prefix = ''
if includePrefix: prefix = 'P='
if (P is None) or (P > 1.0) or (P < 0) or math.isnan(P):
Pstr = prefix + "???"
elif P > 0.1:
Pstr = prefix + "{:.2f}".format(P)
elif P == 0.1:
Pstr = prefix + "0.1"
elif P > 0.01:
Pstr = prefix + "{:.3f}".format(P)
elif P == 0.01:
Pstr = prefix + "0.01"
elif P > 0.001:
Pstr = prefix + "{:.4f}".format(P)
elif P == 0.001:
Pstr = prefix + "0.001"
elif P > 0.0001:
#Pstr = prefix + "{:.5f}".format(P)
Pstr = prefix + "{:.1e}".format(P)
elif P == 0.0001:
#Pstr = prefix + "0.0001"
Pstr = prefix + "1e-4"
elif P > 1e-50:
Pstr = prefix + "{:.1e}".format(P)
else:
Pstr = "P ≈ 0"
# strip leading 0s in exponent
Pstr = Pstr.replace("e+0", "e+")
Pstr = Pstr.replace("e-0", "e-")
return Pstr
### t-test
# pass means_df that already contains the averaged components from all pixels/kernels, for each subject
# see stats_tTest_raw for passing raw kernel-wise data
def stats_tTest(means_df,defaultLabel='PC1',twoTail=True, classList=[0,1]):
import math
from scipy import stats
from sklearn import metrics
try:
means2class_df = means_df[means_df['classID'].isin(classList)] # drop all rows that are not in classList
means2class_df.reset_index(inplace=True, drop=True) # don't forget this if you want to re-index rows from 0
a = means2class_df.loc[means2class_df['classID'] == classList[0]][defaultLabel] # extract all lda_means where classID=0
b = means2class_df.loc[means2class_df['classID'] == classList[1]][defaultLabel] # extract all lda_means where classID=1
MD = abs(round(np.average(b) - np.average(a),3))
if twoTail:
tStat, pValue = stats.ttest_ind(a, b, equal_var = True, alternative='two-sided') #run independent 2 sample T-Test
else:
tStat, pValue = stats.ttest_ind(a, b, equal_var = True, alternative='less') #a < b
P = pValue
Pstr = formatPstring(P)
# AUC https://scikit-learn.org/stable/modules/generated/sklearn.metrics.auc.html
labels_0 = np.zeros(len(a)).astype(int)
labels_1 = np.ones(len(b)).astype(int)
labels = np.concatenate((labels_0,labels_1))
aANDb = np.concatenate((a,b))
y = labels
pred = aANDb
fpr, tpr, thresholds = metrics.roc_curve(y, pred, pos_label=1)
AUC = metrics.auc(fpr, tpr)
if AUC == 1.0:
AUCstr = 'AUC=1.0'
else:
AUCstr = 'AUC={:.2f}'.format(AUC)
CM = max(0,AUC-0.5)*(-math.log(max(P,1e-32))) # composite metric: AUC and P only (limit to 1e-32 to avoid math error)
#CM = round(AUC*MD*(-math.log(P))) # composite metric
#CM = round(AUC*MD*(CA*100-50)*(-math.log(P))) # composite metric, this should match Orange
if P < 0:
CM = -1
CMstr = 'CM:ERR (P<0)'
elif P == 0:
CM = 1e32
CMstr = 'CM≈ ∞ (P=0)'
else:
CM = max(0,AUC-0.5) * (-math.log(P))
CMstr = 'CM={:.1f}'.format(CM)
except: # catch malformed inputs eg. a single class, too few Ns, etc
P = -1
Pstr = '*** ERROR ***'
CM = -1
CMstr = '*** ERROR ***'
MD = -1
AUC = -1
AUCstr = '*** ERROR ***'
return P, Pstr, CM, CMstr, MD, AUC, AUCstr
### t-test on raw kernel-wise (not averaged) input data
# pass raw_df that containes kernel-wise (not averaged)data, for each subject
# nComponents: # PCs to calculate, must be <= # features, >= than the default colLabel requested for the t-test
def stats_tTest_raw(raw_df, defaultLabel='PC1', twoTail=True, classList=[0,1]):
means2class_df = raw_df[raw_df['classID'].isin(classList)] # drop all rows that are not in classList
means2class_df.reset_index(inplace=True, drop=True) # don't forget this if you want to re-index rows from 0
names_unique = means2class_df['name'].unique()
N_names_unique = names_unique.shape[0] # how many subjects?
classIDarr = np.zeros(N_names_unique).astype(int) # classIds will match names_unique
# prepare a df to accept the subject-wise PCA means
PCAmeansDF = pd.DataFrame(0.0,index=range(N_names_unique),columns=[defaultLabel]) # 0.0 forces floats
for nameCtr in range(N_names_unique):
n = names_unique[nameCtr]
subarr = means2class_df[means2class_df['name'] == n] # extract all rows from raw_df for this name
subarr.reset_index(inplace = True, drop = True)
classID = subarr['classID'].iloc[0]
classIDarr[nameCtr] = classID
mean_pca_n = subarr[defaultLabel].mean()
PCAmeansDF[defaultLabel].at[nameCtr] = mean_pca_n
# assemble results into a df
names_df = pd.DataFrame(data=names_unique,columns=['name'])
classID_df = pd.DataFrame(data=classIDarr,columns=['classID'])
PCAout_df = pd.concat([names_df,classID_df,PCAmeansDF], axis = 1)
# now do the t-test on the single requested PC
return stats_tTest(PCAout_df,defaultLabel=defaultLabel,twoTail=twoTail, classList=classList)
def stats_tTest_np(group1,group2, twoTail=True):
# group1&2 are 1D np vectors
# twoTail=False: group2 > group1
import math
from scipy import stats
try:
if twoTail:
tStat, pValue = stats.ttest_ind(group1.reshape(-1), group2.reshape(-1), equal_var = True, alternative='two-sided') # run independent 2 sample T-Test; .reshape(-1) because we may pass a (N,1) 2D array
else:
tStat, pValue = stats.ttest_ind(group1.reshape(-1), group2.reshape(-1), equal_var = True, alternative='less') # group1 < group2
P = pValue
Pstr = formatPstring(P)
except:
P = -1
Pstr = '*** ERROR ***'
finally:
return P, Pstr
### MANOVAs
def stats_MANOVA(nComponents,pca_df,defaultLabel='PC',PClist=[], epsilon=1e-7): # epsilon=1e-7: anything smaller and the MANOVA P again returns 1.0
# compute SUBJECT-WISE PCA component means
# pass nComponents=0 and a non-empty PClist if you want to select specific PCs for the MANOVA
# defaultLabel (col names) cannot begin with a number
# see stats_MANOVA_np for adding jitter
try:
pca2class_df = pca_df.drop(pca_df[pca_df['classID'] > 1].index) # drop all rows that are not classID 0 or 1
pca2class_df.reset_index(inplace=True, drop=True) # don't forget this if you want to re-index rows from 0
# names_unique = pca_df['name'].unique()
names_unique = pca2class_df['name'].unique()
N_names_unique = names_unique.shape[0] # how many subjects?
classIDarr = np.zeros(N_names_unique).astype(int) # classIds will match names_unique
# prepare a df to accept the subject-wise PCA means
if nComponents == 0: # specific PCs from a list
PCAmeansDF = pd.DataFrame(0.0,index=range(N_names_unique),columns=[defaultLabel+'%i' % PClist[i] for i in range(len(PClist))]) # 0.0 forces floats
else:
PCAmeansDF = pd.DataFrame(0.0,index=range(N_names_unique),columns=[defaultLabel+'%i' % i for i in range(1,nComponents+1)]) # 0.0 forces floats
for nameCtr in range(N_names_unique):
n = names_unique[nameCtr]
subarr = pca_df[pca_df['name'] == n] # extract all rows from pca_df for this name
subarr.reset_index(inplace = True, drop = True)
classID = subarr['classID'].iloc[0]
classIDarr[nameCtr] = classID
if nComponents == 0: # specific PCs from a list
for componentCtr in PClist:
pca_label = defaultLabel+str(componentCtr)
# printSB(pca_label)
mean_pca_n = subarr[pca_label].mean()
if epsilon>0: mean_pca_n += np.random.uniform(-epsilon, epsilon)*mean_pca_n
PCAmeansDF[pca_label].at[nameCtr] = mean_pca_n
else: # 1st nComponents PCs
for componentCtr in range(nComponents):
i=componentCtr+1
pca_label = defaultLabel+'%i' % i
mean_pca_n = subarr[pca_label].mean()
if epsilon>0: mean_pca_n += np.random.uniform(-epsilon, epsilon)*mean_pca_n
PCAmeansDF[pca_label].at[nameCtr] = mean_pca_n
# assemble results into a df
names_df = pd.DataFrame(data=names_unique,columns=['name'])
classID_df = pd.DataFrame(data=classIDarr,columns=['classID'])
PCAout_df = pd.concat([names_df,classID_df,PCAmeansDF], axis = 1)
# compute MANOVA on PCA means
if (np.count_nonzero(classIDarr==0) < 3) or (np.count_nonzero(classIDarr==1) < 3): # must have a minimum in each class?
raise Exception
independent_variable = 'classID'
formulaStr = ''
if nComponents == 0: # specific PCs from a list
for componentCtr in PClist:
pca_label = defaultLabel+str(componentCtr)
formulaStr = formulaStr + pca_label
if componentCtr != PClist[-1]:
formulaStr = formulaStr + ' + '
else:
for componentCtr in range(nComponents):
i=componentCtr+1
# pca_label = 'PC%i' % i
pca_label = defaultLabel+'%i' % i
formulaStr = formulaStr + pca_label
if componentCtr < nComponents-1:
formulaStr = formulaStr + ' + '
formulaStr = formulaStr + ' ~ ' + independent_variable
# printSB(formulaStr)
fit = MANOVA.from_formula(formulaStr, data=PCAout_df)
# extract the P value
test_id = 1 # index of testID e.g. Pillai's trace = 1 from above table
P = fit.mv_test().results[independent_variable]['stat'].values[test_id, 4] # 4 is index of the P value column
Pstr = formatPstring(P)
except: # catch malformed inputs eg. a single class, or N<3, etc
P = -1
Pstr = '*** ERROR ***'
PCAout_df = None
return P, Pstr, PCAout_df
def stats_MANOVA_np(group1, group2, epsilon=1e-7): # epsilon=1e-7: anything smaller and the MANOVA P again returns 1.0
# group1&2 are 2D np arrays (Nsamples,NdependentVariables)
try:
if (len(group1) < 3) or (len(group2) < 3): return -1,'*** ERROR ***'
if epsilon > 0: # because MANOVA will incorrectly return P=1 when all values in a column are equal, apply a tiny, random jitter
random_values = np.random.uniform(-epsilon, epsilon, size=group1.shape) * group1 # * group1 ensures that the jitter is always 1e-7x smaller than the magnitude of each element, ensuring that we don't appreciably alter the statistics
# Add the random values to group_0_data
group1_jitter = group1.copy() + random_values
random_values = np.random.uniform(-epsilon, epsilon, size=group2.shape) * group2
group2_jitter = group2.copy() + random_values
else:
group1_jitter = group1
group2_jitter = group2
# Combine the arrays and create a grouping factor
data = np.concatenate((group1_jitter, group2_jitter), axis=0)
groups = np.repeat(['Group 1', 'Group 2'], repeats=[len(group1_jitter), len(group2_jitter)])
# Create a DataFrame for MANOVA (dynamically generate column names)
num_columns = data.shape[1] # Get the number of columns
column_names = [f'Var{i+1}' for i in range(num_columns)]
df = pd.DataFrame(data, columns=column_names)
df['Group'] = groups
# Construct the formula string dynamically
formula_str = ' + '.join(column_names) + ' ~ Group'
# Fit the MANOVA model
manova_model = MANOVA.from_formula(formula_str, data=df)
# Print the MANOVA results
#printSB(manova_model.mv_test())
# Extract and print the P-values
results_summary = manova_model.mv_test()
pillai_trace_pvalue = results_summary.results['Group']['stat'].loc['Pillai\'s trace']['Pr > F']
wilks_lambda_pvalue = results_summary.results['Group']['stat'].loc['Wilks\' lambda']['Pr > F']
hotelling_trace_pvalue = results_summary.results['Group']['stat'].loc['Hotelling-Lawley trace']['Pr > F']
roy_largest_root_pvalue = results_summary.results['Group']['stat'].loc['Roy\'s greatest root']['Pr > F']
P = pillai_trace_pvalue
Pstr = formatPstring(P)
except:
P = -1
Pstr = '*** ERROR ***'
finally:
return P, Pstr
def stats_MANOVA_Xy_np(X, y):
# pass both classes in a single 2D array X, y is used to separate by class
unique_classes = np.unique(y)
if len(unique_classes)>2: raise ValueError(f'stats_MANOVA_Xy_np: y must have only 2 unique classes (unique_classes={unique_classes})')
X0 = X[y == unique_classes[0]]
X1 = X[y == unique_classes[1]]
return stats_MANOVA_np(X0, X1)
### MANOVA on univariate Duetta spectr
def stats_MANOVA_Duetta(nComponents, pca_df, classList=[0,1], defaultLabel='PC'): # compute MANOVA on PCA means
try:
pca2class_df = pca_df[pca_df['classID'].isin(classList)] # drop all rows that are not in classList
pca2class_df.reset_index(inplace=True, drop=True) # don't forget this if you want to re-index rows from 0
names_unique = pca2class_df['name'].unique()
N_names_unique = names_unique.shape[0] # how many subjects?
classIDarr = pca2class_df['classID']
if (np.count_nonzero(classIDarr==classList[0]) < 3) or (np.count_nonzero(classIDarr==classList[1]) < 3): # must have a minimum in each class?
raise Exception
independent_variable = 'classID'
formulaStr = ''
for componentCtr in range(nComponents):
i=componentCtr+1
pca_label = defaultLabel+'%i' % i
formulaStr = formulaStr + pca_label
if componentCtr < nComponents-1:
formulaStr = formulaStr + ' + '
formulaStr = formulaStr + ' ~ ' + independent_variable
# printSB(formulaStr)
fit = MANOVA.from_formula(formulaStr, data= pca2class_df)
# extract the P value
test_id = 1 # index of testID e.g. Pillai's trace = 1 from above table
P = fit.mv_test().results[independent_variable]['stat'].values[test_id, 4] # 4 is index of the P value column
Pstr = formatPstring(P)
except: # catch malformed inputs eg. a single class, or N<3, etc
P = -1
Pstr = '*** ERROR ***'
return P, Pstr
# calculates a composite metric. Higher values of p_value_logbase emphasize AUC over P value. p_value_logbase = 0 to use AUC only
def calcCM(p_value, AUC_value, accuracy=-1, p_value_logbase=10, multiplier=1):
# pass AUC_value=-1 or accuracy=-1 to omit from CM calc
# new multiplier can be passed to multiply the computed CM eg. inter-group delta
if accuracy < 0:
acc = 1.0 # omit from calc
else:
acc = accuracy
if AUC_value < 0:
AUC = 1.0 # omit from calc
else:
AUC = AUC_value
if math.isnan(p_value):
CM = -1
CMstr = 'CM: ERR (P<0)'
elif p_value_logbase <= 0: # AUC-only
CM = AUC
CM *= multiplier
CMstr = 'CM={:.2f}'.format(CM)
elif p_value < 0:
CM = -1
CMstr = 'CM: ERR (P<0)'
elif p_value == 0:
CM = 100 # arbitrary ceiling
CMstr = 'CM≈ ∞ (P=0)'
else:
CM = AUC * (-math.log(p_value,p_value_logbase)) * acc # 2024-03-24: now log base 10 to further de-emphasize effect of p_value
CM *= multiplier # typically absMeanGrpDiff/2 to capture spread between the group means after regression
CMstr = 'CM={:.1f}'.format(CM)
return CM, CMstr
### conditions raw Duetta/Mk5 spectra
def conditionSpectra(spectralData, # input df, 1st 3 cols typically name, classID, T
normMode, # 0: don't normalize; 1: normalize each spectrum to its own peak; 2: normalize each T series to T0 spectrum
firstSpectralColumnIx, # 0-based index of 1st wavelength column
baselineStart_nm, # baseline extends from baselineStart_nm baselineEnd_nm; enter 0 to begin at start of spectrum
baselineEnd_nm, # enter 0 for no baseline subtraction
spectrumStart_nm, # nominal start of spectrum with leading baseline omitted; enter 0 to include entire spectrum including baseline
spectrumEnd_nm, # nominal end of spectrum; enter 0 to include entire spectrum to end
Tmin=1, # a subset of time series is extracted: this is the first t-point (1-based index)
Tmax=-1, # a subset of time series is extracted: this is the last t-point (1-based index); -1 for all
SG_w_size = 0, # Sav-Gol filter window, 0 for no filtering
SG_polyorder = 3, # Sav-Gol polyorder
decimate_q = 1, # decimate factor if > 1 (https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.decimate.html)
replaceNameUnderscoresWithDashes = False, # _ represent replicates of same subject, set to true to treat these as different subjects
variablesColID='T', # label of true time dimension ie complete spectra acquired at nT time points during the PC
graphOutPDF = '' # file path to save conditioned spectral overlay
):
from scipy import signal
# printSB('Conditioning input spectra...')
# extract the wavelength vector
cols = spectralData.columns
nm = cols[firstSpectralColumnIx:].astype('float') # wavelengths
if Tmax == -1: Tmax = spectralData[variablesColID].max()+1 # 1-based
spectralDataSubset = spectralData.loc[spectralData[variablesColID] < Tmax] # extract subset of T-points
spectralDataSubset = spectralDataSubset.loc[spectralDataSubset[variablesColID] >= Tmin-1]
spectralDataSubset.reset_index(inplace=True,drop=True)
### extract the row metas
IDarr = spectralDataSubset.loc[:,('name','classID',variablesColID)] # variablesColID = 'T'
if replaceNameUnderscoresWithDashes:
IDarr['name'] = IDarr['name'].str.replace('_','-')
spectralDataSubset = spectralDataSubset.iloc[:,firstSpectralColumnIx:] # we will operate on just the spectral data for simplicity
### Sav-Gol filter for normalizing (to avoid spurious peaks)
polyorder = SG_polyorder # polynomial order [3]
w_size = SG_w_size
if SG_w_size == 0: # if not filtering of raw spectra use this default filter for norm peak detection only
w_size = int(nm.shape[0] / 20) # window size [33]
w_size = max(SG_polyorder +1,w_size) # polyorder must be less than window_length
### baseline subtract
spectralDataSubset_bsub = spectralDataSubset.copy()
if nm[0] < baselineEnd_nm:
baselineIndexesGT = np.where(nm>=baselineStart_nm)[0]
baselineIndexesLT = np.where(nm<=baselineEnd_nm)[0]
baselineIndexes = np.intersect1d(baselineIndexesGT, baselineIndexesLT)
firstBaselineIx = baselineIndexes[0] # 0-based ix of first bin of baseline
lastBaselineIx = baselineIndexes[-1] # 0-based ix of last bin of baseline
# subtract baseline
for j in range(spectralDataSubset.shape[0]):
# baseline = spectralDataSubset.iloc[j,firstBaselineIx:lastBaselineIx].mean(axis=0)
baseline = spectralDataSubset.iloc[j,firstBaselineIx:lastBaselineIx].mean(axis=0)
spectralDataSubset_bsub.iloc[j,:] -= baseline
### normalize
spectralDataSubset_bsub_norm = spectralDataSubset_bsub
# printSB(spectralDataSubset_bsub_norm.max().max())
if normMode == 0: # don't normalize
pass # do nothing more
elif normMode == 1: # normalize each spectrum to its own peak
spectralDataSubset_bsub_filt = spectralDataSubset_bsub.copy()
for rowCtr in range(spectralDataSubset_bsub.shape[0]):
y = spectralDataSubset_bsub.iloc[rowCtr,:].to_numpy()
y_filt = signal.savgol_filter(y, w_size, polyorder, mode="nearest")
spectralDataSubset_bsub_filt.iloc[rowCtr,:] = y_filt
# normalize each spectrum to its own (filtered) peak value
for j in range(spectralDataSubset_bsub.shape[0]):
maxRowVal = spectralDataSubset_bsub_filt.iloc[j,:].max(axis=0)
spectralDataSubset_bsub_norm.iloc[j,:] /= maxRowVal
else: # normalize each spectral series to Tmin
# first extract all T0 spectra
# spectralDataSubset_bsub_T0 = spectralDataSubset_bsub.loc[spectralDataSubset_bsub[variablesColID]==0] # all T0 spectra
spectralDataSubset_bsub_Tmin = spectralDataSubset_bsub[IDarr[variablesColID] == Tmin-1] # extract all rows from spectralDataSubset_bsub where corresponding rows in IDarr['T'] == Tmin-1
spectralDataSubset_bsub_Tmin.reset_index(inplace=True,drop=True)
IDarr_Tmin = IDarr.loc[IDarr[variablesColID] == Tmin-1] # all Tmin meta
IDarr_Tmin.reset_index(inplace=True,drop=True)
# S-G filter T0 spectra to get true noise-free maxval for norm
spectralDataSubset_bsub_Tmin_filt = spectralDataSubset_bsub_Tmin.copy()
for rowCtr in range(spectralDataSubset_bsub_Tmin.shape[0]):
y = spectralDataSubset_bsub_Tmin.iloc[rowCtr,:].to_numpy()
y_filt = signal.savgol_filter(y, w_size, polyorder, mode="nearest")
spectralDataSubset_bsub_Tmin_filt.iloc[rowCtr,:] = y_filt
for nameCtr in range(IDarr_Tmin.shape[0]):
name = IDarr_Tmin.iloc[nameCtr,IDarr_Tmin.columns.get_loc('name')]
maxVal_Tmin = spectralDataSubset_bsub_Tmin_filt.iloc[nameCtr,:].max(axis=0)
#printSB(name + ' max: ' + str(maxVal_T0))
for rowCtr in range(spectralDataSubset_bsub_norm.shape[0]):
if IDarr.iloc[rowCtr,IDarr.columns.get_loc('name')] == name: # does the row name match name?
spectralDataSubset_bsub_norm.iloc[rowCtr,:] /= maxVal_Tmin
#printSB(spectralDataSubset_bsub_norm.max().max())
### S-G filter raw spectral data
#printSB(firstSpectralColumnIx)
if SG_w_size > 0:
for rowCtr in range(spectralDataSubset_bsub_norm.shape[0]):
y = spectralDataSubset_bsub_norm.iloc[rowCtr,:].to_numpy()
y_filt = signal.savgol_filter(y, SG_w_size, polyorder, mode="nearest")
spectralDataSubset_bsub_norm.iloc[rowCtr,:] = y_filt
### downsample
if decimate_q > 1:
y = spectralDataSubset_bsub_norm.iloc[0,:].to_numpy()
y_filt = signal.decimate(y, q=decimate_q) # dummy (ftype = "fir" doesn't work!)
#printSB(y_filt.shape)
# Initialize a 2D Pandas DataFrame with all zeroes
spectralDataSubset_bsub_norm_dec = pd.DataFrame(0, index=np.arange(spectralDataSubset_bsub_norm.shape[0]), columns=np.arange(y_filt.shape[0]))
for rowCtr in range(spectralDataSubset_bsub_norm.shape[0]):
y = spectralDataSubset_bsub_norm.iloc[rowCtr,:].to_numpy()
y_filt = signal.decimate(y, q=decimate_q, n=SG_w_size) # ftype = "fir" doesn't work!
spectralDataSubset_bsub_norm_dec.iloc[rowCtr,:] = y_filt
# nm = signal.decimate(nm, q=decimate_q, n=0) # must decimate the nm vector too! This imparts a significant phase shift int he nm values
nm_dec = nm[::decimate_q]
spectralDataSubset_bsub_norm = spectralDataSubset_bsub_norm_dec
else:
nm_dec = nm
spectralDataSubset_bsub_norm.columns = spectralDataSubset_bsub_norm.columns.astype(str) # make sure col names are type str
# rename columns with nm
# nm_str = nm_dec.astype(str)
# spectralDataSubset_bsub_norm.columns.values[firstSpectralColumnIx:] = nm_str ### THIS SHOULD WORK BUT CRASHES THE KERNEL!!!
for j in range(nm_dec.shape[0]):
spectralDataSubset_bsub_norm.columns.values[j] = "{:.1f}".format(nm_dec[j])
### extract required wavelength subrange
if spectrumStart_nm == 0: spectrumStart_nm = nm_dec[0] # from start
if spectrumEnd_nm == 0: spectrumEnd_nm = nm_dec[-1] # to end
specIndexesGT = np.where(nm_dec>=spectrumStart_nm)[0]
specIndexesLT = np.where(nm_dec<=spectrumEnd_nm)[0]
specIndexes = np.intersect1d(specIndexesGT, specIndexesLT)
first_nmIx = specIndexes[0]
last_nmIx = specIndexes[-1]
spectraExtracted = spectralDataSubset_bsub_norm.iloc[:,first_nmIx:last_nmIx+1]
spectraExtracted.reset_index(inplace=True,drop=True)
#printSB(spectraExtracted.max().max())
### reassemble into a final output df
conditioned_df = pd.DataFrame(data=IDarr)
conditioned_df = pd.concat([conditioned_df,spectraExtracted], axis = 1)
#printSB(conditioned_df.iloc[:,firstSpectralColumnIx:].max().max())
# nm_conditioned = spectraExtracted.columns.astype('float') # wavelengths
nm_conditioned = nm_dec[first_nmIx:last_nmIx+1] # wavelength subrange
if graphOutPDF != '': # plot conditioned UNfiltered spectra
plt.style.use('bmh')
plt.figure(figsize=(6, 4))
spectraExtracted_T0 = conditioned_df.loc[conditioned_df[variablesColID] == 0] # extract T0 spectra
spectraExtracted_T0.reset_index(inplace=True,drop=True)
# there must be a vectorized way to do the below:
labels = spectraExtracted_T0['classID']
colors = ['green','red']
colorArr = []
for i in range(labels.shape[0]):
colorArr.append(colors[labels[i]])
x = nm_conditioned
# we plot the T0 spectra only
for i in range(spectraExtracted_T0.shape[0]):
y = spectraExtracted_T0.iloc[i,firstSpectralColumnIx:].transpose()
plt.plot(x, y, color=colorArr[i],linewidth=1.0)
#plt.show()
plt.savefig(graphOutPDF)
printSB('Conditioned T0 spectral graph saved to: ' + graphOutPDF)
plt.close()
conditioned_df.reset_index(inplace=True, drop=True)
return conditioned_df, nm_conditioned
### extract root name e.g. abc_def_1 returns abc_def; abc-def-1 returns abc-def-1,
def extractRootname(name, separator = '_'):
rpart = name.rpartition(separator)
if rpart[0] == '':
# return rpart[2],'' # separator does not exist: complete name is in [2]
return rpart[2] # separator does not exist: complete name is in [2]
else:
# return rpart[0],rpart[2] # [0] contains everything ahead of last separator, [2] contains suffix after last seprator
return rpart[0] # [0] contains everything ahead of last separator
### extract root names from a column named 'name' in df
### use this to get subject rootnames from an augmented dataset (must call BEFORE underscore replacement ie: if kReplaceNameUnderscoresWithDashes: imageData['name'] = imageData['name'].str.replace('_','-') )
def extractRootnames(df, separator = '_AUG_'): # df must have a 'name' and 'classID' column
rootNames = []
rootClassIDs = []
uniqueNames = df['name'].unique()
for n in uniqueNames:
rootNames.append(extractRootname(n, separator))
classID_value = df.loc[df['name'] == n, 'classID'].values[0]
rootClassIDs.append(classID_value)
rootNames_df = pd.DataFrame(list(zip(rootNames, rootClassIDs)), columns=['name','classID'])
rootNames_df.drop_duplicates(subset='name', inplace=True) # extract unique names in df
rootNames_df.sort_values(by=['classID', 'name'], inplace=True)
rootNames_df.reset_index(drop=True, inplace=True)
rootNames = list(set(rootNames)) # extract unique names in list
return rootNames, rootNames_df # return as a list and a df
### summarizes the extracted root subjects and replicate instances for a set of Duetta/Mk5 input spectra
def summarizeSubjects(input_df): # pass a pandas df with a column 'name'
instance_names_all = input_df['name']
instance_names_all
instance_classIDs_all = input_df['classID']
instance_classIDs_all
instance_names_unique = instance_names_all.unique()
instance_names_unique
instance_classIDs_unique = []
for unique_name in instance_names_unique:
for k in range(instance_names_all.shape[0]):
if instance_names_all[k] == unique_name:
instance_classIDs_unique.append(instance_classIDs_all[k])
break # continue with next unique_name
nClass0 = instance_classIDs_unique.count(0)
nClass1 = instance_classIDs_unique.count(1)
root_names = []
root_classIDs = []
for j in range(instance_names_unique.shape[0]):
root_names.append(extractRootname(instance_names_unique[j]))
root_classIDs.append(instance_classIDs_unique[j])
root_names_unique = np.unique(root_names)
root_classIDs_unique = []
for unique_rootname in root_names_unique:
for k in range(len(root_names)):
if root_names[k] == unique_rootname:
root_classIDs_unique.append(root_classIDs[k])
break # continue with next unique_rootname
nClass0_root = root_classIDs_unique.count(0)
nClass1_root = root_classIDs_unique.count(1)
maxNameLen = len(max(instance_names_unique, key=len)) # used for padding below
result = 'All instances'.ljust(max(3+maxNameLen+2,len('All instances')+1)) + 'classID (' + str(nClass0) + '+' + str(nClass1) + '):\n'
for j in range(instance_names_unique.shape[0]):
result = result + ' ' + instance_names_unique[j].ljust(max(maxNameLen+2,len('All instances')+1-3)) + str(instance_classIDs_unique[j]) + '\n'
maxRootnameLen = len(max(root_names_unique, key=len)) # used for padding below
result = result + 'Root names'.ljust(max(3+maxRootnameLen+2,len('Root names')+1)) + 'classID (' + str(nClass0_root) + '+' + str(nClass1_root) + '):\n'
for j in range(root_names_unique.shape[0]):
result = result + ' ' + root_names_unique[j].ljust(max(maxRootnameLen+2, len('Root names')+1-3)) + str(root_classIDs_unique[j]) + '\n'
result = result + '\n'
return result
### strips the Orange tokens from col headers
def fixOrangeHeaders_OLD(df):
cols = df.columns
colsAdj = []
for headerItem in cols:
# printSB(headerItem)
if '#' in headerItem:
colsAdj.append(headerItem.rsplit('#')[1])
else:
colsAdj.append(headerItem)
df.columns = colsAdj
return df
# strip Orange tokens in-place, returns the adjusted col names, 0-based ix of 1st lambda col, number of lambda cols
def fixOrangeHeaders(df):
"""
Regular Expression Pattern: The r'^.*#' regular expression does the following:
^: Matches the start of the string (column name).
.*: Matches any character (.) zero or more times (*).
#: Matches the literal '#' character.
Together, this pattern matches everything from the beginning of the column name up to and including the first '#' character it finds.
https://g.co/gemini/share/bc7dc4af6baf
"""
def starts_with_L_and_numeric(s):
# Regular expression pattern to match 'L' followed by a valid numeric value
pattern = r'^L\d+(\.\d+)?$'
# Use re.match to check if the string matches the pattern
return bool(re.match(pattern, s))
columns = df.columns.str.replace(r'^.*#', '', regex=True).tolist()
df.columns = columns
# now auto-detect start and end of lambda cols
first_numeric = -1
last_numeric = -1
for i, col in enumerate(columns):
try:
# Try converting the string to a float
float(col)
if first_numeric == -1: first_numeric = i
last_numeric = i
except ValueError:
pass # Ignore non-numeric strings
# lambda cols might be 'L400', 'L410', etc
if first_numeric == -1: # we never found a numeric column with orange tokens above
for i, col in enumerate(columns):
if starts_with_L_and_numeric(col):
if first_numeric == -1: first_numeric = i
last_numeric = i
columns[i] = columns[i].lstrip('L')
if first_numeric > 0: df.columns = columns # update the column names with the stripped Ls if we found any
firstLambdaCol = first_numeric
nLambdas = last_numeric - first_numeric + 1
return columns, firstLambdaCol, nLambdas # columns is a list of string
def get_child_files(directory, fIn='', extension='csv'): # optional fIn returned if directory=''; typically pass dirIN, csvIN
# returns a list of fullpath(s)
# new version: fIn can contain a tab-delim list of full paths to accommodate the new csvInList() array in SB
file_paths = []
# non-nil enclosing directory takes precedence
if directory!='': # extract all csv files from directory
for root, dirs, files in os.walk(directory):
for file in files:
if file.endswith('.' + extension): file_paths.append(os.path.join(root, file))
elif "\t" in fIn: # if fIn contains a tab it's a list of files
file_paths = fIn.split("\t")
else: # append only a single fIn to list
file_paths.append(fIn)
return file_paths
### load a flow cytometry 3.0 fcs file
# call with path to fcs3.0 file, returns a dataframe with raw data as floats
# Time column is removed
# min/maxLimits needed to strip occasional rogue values
def load_fcs(fcsPath, minRogueLimit=1, maxRogueLimit=1e6):
f = open(fcsPath, "rb")
# The format name is encoded in 6 letters
# An ASCII letter is coded with one octet
#file_format = "".join([f.read(1) for __ in range(6)])
file_format = f.read(6).decode('ascii')
#sys.stdout.write("Format: %s\n" % file_format)
#printSB("Format: %s\n" % file_format)
# The format descriptions reserves 4 octets that we skip
skip = f.read(4)
# 8 octet chunks encode the start and end positions
# of different parts of the data
text_start = int(f.read(8).decode('ascii').strip())
text_end = int(f.read(8).decode('ascii').strip())
data_start = int(f.read(8).decode('ascii').strip())
data_end = int(f.read(8).decode('ascii').strip())
analysis_start = int(f.read(8).decode('ascii').strip())
analysis_end = int(f.read(8).decode('ascii').strip())
####################################################
# Here starts the parsing of the "TEXT" portion #
# which describes how the data proper is organized #
####################################################
f.seek(text_start)
# The first character in the primary TEXT segment is the ASCII delimiter character.
sep = f.read(1).decode('ascii')
text_segment = f.read(text_end - text_start).decode('ascii')
fields = text_segment.split(sep)
info = {} # dictionary
i = 0
while i < len(fields) - 1:
key = fields[i]
i += 1
val = fields[i]
i += 1
# Keywords are case insensitive, they may be written in a file in lower case, upper case, or a
# mixture of the two. However, an FCS file reader must ignore keyword case. A keyword value may
# be in lower case, upper case or a mixture of the two. Keyword values are case sensitive.
info[key.upper()] = val
# extract parameter names
parameters = []
# indices of the parameters
p_indices = range(1, int(info["$PAR"]) + 1)
for i in p_indices:
p_name = info["$P%dN" % i]
parameters.append(p_name)
if info["$BYTEORD"] == "4,3,2,1":
endianness = ">"
else:
endianness = "<"
assert info["$BYTEORD"] == "1,2,3,4"
# Type of data:
if info["$DATATYPE"] == "F":
nParams = int(info["$PAR"])
nEvents = int(info["$TOT"]) # no of floats should be nParams * nEvents
f.seek(data_start)
flatFloats = np.fromfile(f, dtype=np.float32, count=nParams*nEvents).byteswap().newbyteorder(sys.byteorder)
# Reshape the array into a 2D array with nEvents rows and 8 columns
reshaped_flatFloats = flatFloats.reshape(nEvents, nParams)
# Convert the reshaped array into a DataFrame
df = pd.DataFrame(data=reshaped_flatFloats,columns=parameters)
if 'Time' in df.columns: df.drop('Time', axis=1, inplace=True) # Remove the 'Time' column, messes up limits below
df = df[~(df > maxRogueLimit).any(axis=1)]
df = df[~(df < minRogueLimit).any(axis=1)]
df.reset_index(inplace=True, drop=True)
else:
df = pd.DataFrame() # empty placeholder
printSB('*** Only float data are supported ***')
f.close()
return df
def cleanAUGfname(fBaseNoExt): # removes the unique '_NNNN' suffix appended to AUGnX fnames
if (fBaseNoExt[-6:-4] == 'X_') and (fBaseNoExt[-10:-7] == 'AUG'):
return fBaseNoExt[0:-5]
else:
return fBaseNoExt
### computes support vector classification or regression on X & Y
### needs global vars: maxIter_svm, param_grid for grid search
def SVCR(kSVregression, X, Y, kGridSearch=True, SVhyperparameters = [10, 1, 'linear', 0.0001], n_splits=10, n_repeats=5): # SVhyperparameters: C, gamma, kernel (other kernel options: 'rbf'), ignored when kGridSearch=True
from sklearn.svm import SVC,SVR
from sklearn.model_selection import GridSearchCV, RepeatedKFold
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
C_svm = SVhyperparameters[0]
gamma_svm = SVhyperparameters[1]
kernel_svm = SVhyperparameters[2]
tol_svm = SVhyperparameters[3]
rkf = RepeatedKFold(n_splits=min(n_splits,X.shape[0]), n_repeats=n_repeats)
printSB()
if kSVregression: # SVR
if kGridSearch == False:
svm_model = SVR(kernel=kernel_svm,verbose=3,C=C_svm, gamma=gamma_svm, max_iter=maxIter_svm, tol=tol_svm) # https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVR.html
svm_model.fit(X, Y)
#predictions
y_pred = svm_model.predict(X) # make predictions on the input set to check accurcy of linear separability
else: # grid search
svm_model = SVR() # Create the SVR regressor
grid_search = GridSearchCV(estimator=svm_model, param_grid=param_grid, cv=rkf)
grid_search.fit(X, Y)
printSB("Grid-search best SVR parameters (NOTE: rbf kernel tends to overfit, check with RAND):", grid_search.best_params_)
y_pred = grid_search.predict(X) # make predictions on the input set to check accurcy of linear separability
decision_function_df = pd.DataFrame(data=y_pred,columns=['score'])
accPercent = -1 # undefined for SVR
else: # SVC
if kGridSearch == False:
# svm_model = SVC(kernel= kernel_svm,verbose=3,C=C_svm, gamma=gamma_svm, max_iter=maxIter_svm, tol=tol_svm, random_state=42) # https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.html?highlight=svc#sklearn.svm.SVC
svm_model = SVC(kernel= kernel_svm,verbose=3,C=C_svm, gamma=gamma_svm, max_iter=maxIter_svm) # https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.html?highlight=svc#sklearn.svm.SVC
svm_model.fit(X, Y)
#predictions
y_pred = svm_model.predict(X) # make predictions on the input set to check accurcy of linear separability
decision_function = svm_model.decision_function(X)
if kStandardScaleResult:
decision_function = scaler.fit_transform(decision_function.reshape(-1, 1)) # often decision_function contains tiny values so rescaling standardizes when averaging with different hyperparameters
else: # grid search
svm_model = SVC() # Create the SVM classifier
grid_search = GridSearchCV(estimator=svm_model, param_grid=param_grid, cv=rkf)
grid_search.fit(X, Y)
printSB("Grid-search best SVC parameters (NOTE: rbf kernel tends to overfit, check with RAND):", grid_search.best_params_)
y_pred = grid_search.predict(X) # make predictions on the input set to check accurcy of linear separability
decision_function = grid_search.decision_function(X)
if kStandardScaleResult:
decision_function = scaler.fit_transform(decision_function.reshape(-1, 1)) # often decision_function contains tiny values so rescaling standardizes when averaging with different hyperparameters
decision_function_df = pd.DataFrame(data=decision_function,columns=['score'])
accPercent = 100*accuracy_score(Y, y_pred)
printSB()
return decision_function_df, accPercent
### computes support vector regression with cross-validation on X & Y
### needs global vars: param_grid for grid search
def SVR_CV(X, Y, n_splits=10, n_repeats=5):
from sklearn.svm import SVR
from sklearn.model_selection import train_test_split, GridSearchCV, cross_val_score, RepeatedKFold
### Sandardization of data ###
from sklearn.preprocessing import StandardScaler
PredictorScaler=StandardScaler()
TargetVarScaler=StandardScaler()
# Storing the fit object for later reference
PredictorScalerFit=PredictorScaler.fit(X)
TargetVarScalerFit=TargetVarScaler.fit(Y.values.reshape(-1,1))
# Generating the standardized values of X and y
Xstd=PredictorScalerFit.transform(X)
Ystd=TargetVarScalerFit.transform(Y.values.reshape(-1,1))
# Ystd=TargetVarScalerFit.transform(Y.values)
# Ystd = np.ravel(Ystd)
rkf = RepeatedKFold(n_splits=n_splits, n_repeats=n_repeats)
svm_model = SVR() # Create the SVR regressor
grid_search = GridSearchCV(estimator=svm_model, param_grid=param_grid, cv=rkf, n_jobs=-1)
grid_search.fit(Xstd, Ystd)
printSB("Grid-search best SVR parameters (NOTE: rbf kernel tends to overfit, check with RAND):", grid_search.best_params_)
y_pred = grid_search.predict(Xstd) # make predictions on the input set
# Scaling the predicted Price data back to original price scale
y_pred=TargetVarScalerFit.inverse_transform(y_pred.reshape(-1,1))
decision_function_df = pd.DataFrame(data=y_pred,columns=['score'])
return decision_function_df
### computes support vector classification on X & Y, selects the best hyperparams based on CM (composite metric)
### pass subject-wise mean PCs in PCmeans_df: 1st col is name, 2nd col is classID, rest are PCs
### needs global vars: maxIter_svm, param_grid for grid search, kStandardScaleResult
# see also SVC_CM2
# 2024-06-15: pass an explicit list of 0-based featureColumnIndexes to use for SVC instead of trying to determine these from defaultLabel
def SVC_CM1(PCmeans_df, param_grid, n_splits=10, n_repeats=5, suppressWarnings=False, defaultLabel='PC', featureColumnIndexes=[], p_value_logbase=10):
from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV, RepeatedKFold
from sklearn import metrics
from sklearn.preprocessing import StandardScaler
C_arr = param_grid.get('C')
gamma_arr = param_grid.get('gamma')
kernel = param_grid.get('kernel')
tol = param_grid.get('tol')
rkf = RepeatedKFold(n_splits=min(n_splits,PCmeans_df.shape[0]), n_repeats=n_repeats)
scaler = StandardScaler()
CMbest = -1e32
dict_best = {'CM':-1} # in case no solution/exception
#X = PCmeans_df.iloc[:,2:]
if len(featureColumnIndexes) == 0:
first_PC_index = next((i for i, col in enumerate(PCmeans_df.columns) if col.startswith(defaultLabel)), -1) # find index of 1st PCx or UMAPx etc column
if first_PC_index == -1:
printSB('SVC_CM1: PCx columns could not be found (fatal error)')
return dict_best, pd.DataFrame()
X = PCmeans_df.iloc[:,first_PC_index:]
else: # user passed explicit 0-based feature column indexes
try:
X = PCmeans_df.iloc[:,featureColumnIndexes]
#printSB(X)
except:
printSB('Invalid featureColumnIndexes')
return dict_best, pd.DataFrame()
#Y = PCmeans_df.iloc[:,1]
Y = PCmeans_df['classID']
svm_model = SVC() # Create the SVM classifier
for C in C_arr:
for gamma in gamma_arr:
d = {'C': [C], 'gamma': [gamma], 'kernel': kernel, 'tol': tol} # create a 4-element dictionary for this C,gamma pair
#printSB(d)
grid_search = GridSearchCV(estimator=svm_model, param_grid=d, cv=rkf)
grid_search.fit(X, Y)
y_pred = grid_search.predict(X) # make predictions on the input set to check accurcy of linear separability
decision_function = grid_search.decision_function(X)
if kStandardScaleResult:
decision_function = scaler.fit_transform(decision_function.reshape(-1, 1)) # often decision_function contains tiny values so rescaling standardizes when averaging with different hyperparameters
decision_function_df = pd.DataFrame(data=decision_function,columns=['score'])
# now compute t-test for dec fn, AUC -> CM
a = decision_function_df.loc[PCmeans_df['classID'] == 0]['score']
b = decision_function_df.loc[PCmeans_df['classID'] == 1]['score']
_, pValue = stats.ttest_ind(a, b, equal_var = True, alternative='less') #run independent 2 sample T-Test
y_classification = np.where(decision_function_df['score'] > 0, 1, 0) # so convert it to a predicted 0 or 1 class, depending on <0 vs >0
# AUC
predicted_probabilities = np.where(decision_function_df['score'] > 0, 1, 0) # so convert it to a predicted 0 or 1 class, depending on <0 vs >0
actual_labels = PCmeans_df['classID']
AUC = metrics.roc_auc_score(actual_labels, predicted_probabilities)
if AUC == 1.0:
AUCstr = 'AUC = 1.0'
else:
AUCstr = 'AUC = {:.2f}'.format(AUC)
# accuracy
accuracy = metrics.accuracy_score(actual_labels, y_classification)
if accuracy == 1.0:
accuracyStr = 'accuracy=1.0'
else:
accuracyStr = 'accuracy={:.2f}'.format(accuracy)
# f1 score
f1 = metrics.f1_score(actual_labels, y_classification)
# composite metric
CM, CMstr = calcCM(pValue, AUC, accuracy, p_value_logbase=p_value_logbase) # Higher values of p_value_logbase emphasize AUC over P value; p_value_logbase=0 uses AUC only
#printSB('nuSVC_CM: ',CMstr)
if CM > CMbest:
CMbest = CM
if p_value_logbase <= 0: CM, CMstr = calcCM(pValue, AUC, accuracy, p_value_logbase=10) # for scoring/optimization purposes we use the AUC-only CM, but for reporting we recalc it the regular way so it's comparable to other runs
#dict_best = {'nu':nu, 'pValue':pValue, 'Pstr':formatPstring(pValue), 'AUC':AUC, 'AUCstr':AUCstr, 'accuracy':accuracy, 'accuracyStr':accuracyStr, 'f1':f1, 'CM':CM, 'CMstr':CMstr, 'decision_function_df':decision_function_df, 'y_pred':y_classification}
dict_best = {'C':C, 'gamma':gamma, 'kernel':kernel, 'tol':tol, 'pValue':pValue, 'Pstr':formatPstring(pValue), 'AUC':AUC, 'AUCstr':AUCstr, 'accuracy':accuracy, 'accuracyStr':accuracyStr, 'f1':f1, 'CM':CM, 'CMstr':CMstr, 'decision_function_df':decision_function_df, 'y_pred':y_classification}
#printSB('SVC_CM1:dict_best',dict_best)
return dict_best # return the dict with best results (see above for keys). get elements of dict like this: dict_best['pValue']
### computes nu-support vector classification on X & Y, selects the best hyperparams based on CM (composite metric)
### pass subject-wise mean PCs in PCmeans_df: 1st col is name, 2nd col is classID, rest are PCs
### needs global vars: maxIter_svm, param_grid for grid search?, kStandardScaleResult
# 2024-04-23: CM now includes accuracy as well as log(P) and AUC
# 2024-06-15: pass an explicit list of 0-based featureColumnIndexes to use for SVC instead of trying to determine these from defaultLabel
def SVC_CM2(PCmeans_df, param_grid, n_splits=10, n_repeats=5, suppressWarnings=False, defaultLabel='PC', featureColumnIndexes=[], p_value_logbase=10):
# PCmeans_df can be UMAP means, etc
from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV, RepeatedKFold
from sklearn import metrics
from sklearn.preprocessing import StandardScaler
rkf = RepeatedKFold(n_splits=min(n_splits,PCmeans_df.shape[0]), n_repeats=n_repeats)
scaler = StandardScaler()
dict_best = {'CM':-1} # in case no solution/exception
if len(featureColumnIndexes) == 0:
first_PC_index = next((i for i, col in enumerate(PCmeans_df.columns) if col.startswith(defaultLabel)), -1) # find index of 1st PCx or UMAPx etc column
if first_PC_index == -1:
printSB('SVC_CM2: PCx columns could not be found (fatal error)')
return dict_best, pd.DataFrame()
X = PCmeans_df.iloc[:,first_PC_index:]
else: # user passed explicit 0-based feature column indexes
try:
X = PCmeans_df.iloc[:,featureColumnIndexes]
#printSB(X)
except:
printSB('Invalid featureColumnIndexes')
return dict_best, pd.DataFrame()
Y = PCmeans_df['classID']
svm_model = SVC()
try: # some values of nu are illegal and will throw an exception: catch it and return what we have
grid_search = GridSearchCV(estimator=svm_model, param_grid=param_grid, cv=rkf)
grid_search.fit(X, Y)
y_pred = grid_search.predict(X) # make predictions on the input set to check accurcy of linear separability
decision_function = grid_search.decision_function(X)
if kStandardScaleResult:
decision_function = scaler.fit_transform(decision_function.reshape(-1, 1)) # often decision_function contains tiny values so rescaling standardizes when averaging with different hyperparameters
decision_function_df = pd.DataFrame(data=decision_function,columns=['score'])
#printSB('nuSVC_CM: about to stats.ttest_ind')
# now compute t-test for dec fn, AUC -> CM
a = decision_function_df.loc[PCmeans_df['classID'] == 0]['score']
b = decision_function_df.loc[PCmeans_df['classID'] == 1]['score']
_, pValue = stats.ttest_ind(a, b, equal_var = True, alternative='less') #run independent 2 sample T-Test
#printSB('nuSVC_CM: about to metrics.roc_auc_score')
# AUC
#predicted_probabilities = decision_function_df['score'] THIS IS WRONG: because dec_func varies from <0 to >0 (either side of the hyperplane) it is not a 0-1 probability which is what metrics.roc_auc_score expects
y_classification = np.where(decision_function_df['score'] > 0, 1, 0) # so convert it to a predicted 0 or 1 class, depending on <0 vs >0
actual_labels = PCmeans_df['classID']
AUC = metrics.roc_auc_score(actual_labels, y_classification)
if AUC == 1.0:
AUCstr = 'AUC = 1.0'
else:
AUCstr = 'AUC = {:.2f}'.format(AUC)
#printSB('nuSVC_CM: about to metrics.accuracy_score')
# accuracy
accuracy = metrics.accuracy_score(actual_labels, y_classification)
if accuracy == 1.0:
accuracyStr = 'accuracy=1.0'
else:
accuracyStr = 'accuracy={:.2f}'.format(accuracy)
# f1 score
f1 = metrics.f1_score(actual_labels, y_classification)
# composite metric
CM, CMstr = calcCM(pValue, AUC, accuracy, p_value_logbase=p_value_logbase) # Higher values of p_value_logbase emphasize AUC over P value
dict_best = {'C':grid_search.best_params_['C'], 'gamma':grid_search.best_params_['gamma'], 'kernel':grid_search.best_params_['kernel'], 'tol':grid_search.best_params_['tol'], 'pValue':pValue, 'Pstr':formatPstring(pValue), 'AUC':AUC, 'AUCstr':AUCstr, 'accuracy':accuracy, 'accuracyStr':accuracyStr, 'f1':f1, 'CM':CM, 'CMstr':CMstr, 'decision_function_df':decision_function_df, 'y_pred':y_classification}
except:
if suppressWarnings == False: printSB('SVC_CM2 exception')
# sort by CM
#sorted_dList = sorted(dict_out_list, key=lambda x: x['CM'], reverse=True) # sort by CM: sorted_dList[0] is the best dict
return dict_best # return the dict with best results (see above for keys). get elements of dict like this: d_best.get('decision_function_df')
### computes nu-support vector classification on X & Y, selects the best hyperparams based on CM (composite metric)
### pass subject-wise mean PCs in PCmeans_df: 1st col is name, 2nd col is classID, rest are PCs
### needs global vars: maxIter_svm, param_grid for grid search, kStandardScaleResult
# 2024-04-23: CM now includes accuracy as well as log(P) and AUC
# 2024-05-11: now must pass explict param_grid argument, else risk of altering the global and having obscure bad results
# 2024-06-15: pass an explicit list of 0-based featureColumnIndexes to use for SVC instead of trying to determine these from defaultLabel
def nuSVC_CM(PCmeans_df, param_grid, n_splits=10, n_repeats=5, suppressWarnings=False, defaultLabel='PC', featureColumnIndexes=[], p_value_logbase=10):
# PCmeans_df can be UMAP means, etc
from sklearn.svm import NuSVC
from sklearn.model_selection import GridSearchCV, RepeatedKFold
from sklearn import metrics
from sklearn.preprocessing import StandardScaler
nu_arr = param_grid.get('nu')
kernel = param_grid.get('kernel')
tol = param_grid.get('tol')
rkf = RepeatedKFold(n_splits=min(n_splits,PCmeans_df.shape[0]), n_repeats=n_repeats)
scaler = StandardScaler()
CMbest = -1e32
dict_best = {'CM':-1} # in case no solution/exception
if len(featureColumnIndexes) == 0:
first_PC_index = next((i for i, col in enumerate(PCmeans_df.columns) if col.startswith(defaultLabel)), -1) # find index of 1st PCx or UMAPx etc column
if first_PC_index == -1:
printSB('nuSVC_CM: PCx columns could not be found (fatal error)')
return dict_best, pd.DataFrame()
X = PCmeans_df.iloc[:,first_PC_index:]
else: # user passed explicit 0-based feature column indexes
try:
X = PCmeans_df.iloc[:,featureColumnIndexes]
#printSB(X)
except:
printSB('Invalid featureColumnIndexes')
return dict_best, pd.DataFrame()
Y = PCmeans_df['classID']
# svm_model = SVC() # Create the SVM classifier
# svm_model = NuSVC(tol=0.0001, random_state=42)
svm_model = NuSVC(tol=tol)
try: # some values of nu are illegal and will throw an exception: catch it and return what we have
for nu in nu_arr:
d = {'nu': [nu], 'kernel': kernel, 'tol': [tol]} # create a 2-element dictionary for this C,gamma pair
#printSB('nuSVC_CM: ',d)
grid_search = GridSearchCV(estimator=svm_model, param_grid=d, cv=rkf)
grid_search.fit(X, Y)
y_pred = grid_search.predict(X) # make predictions on the input set to check accurcy of linear separability
decision_function = grid_search.decision_function(X)
if kStandardScaleResult:
decision_function = scaler.fit_transform(decision_function.reshape(-1, 1)) # often decision_function contains tiny values so rescaling standardizes when averaging with different hyperparameters
decision_function_df = pd.DataFrame(data=decision_function,columns=['score'])
#accPercent = 100*accuracy_score(Y, y_pred)
a = decision_function_df.loc[PCmeans_df['classID'] == 0]['score']
b = decision_function_df.loc[PCmeans_df['classID'] == 1]['score']
_, pValue = stats.ttest_ind(a, b, equal_var = True, alternative='less') #run independent 2 sample T-Test
y_classification = np.where(decision_function_df['score'] > 0, 1, 0) # so convert it to a predicted 0 or 1 class, depending on <0 vs >0
actual_labels = PCmeans_df['classID']
AUC = metrics.roc_auc_score(actual_labels, y_classification)
if AUC == 1.0:
AUCstr = 'AUC = 1.0'
else:
AUCstr = 'AUC = {:.2f}'.format(AUC)
#printSB('nuSVC_CM: about to metrics.accuracy_score')
# accuracy
accuracy = metrics.accuracy_score(actual_labels, y_classification)
if accuracy == 1.0:
accuracyStr = 'accuracy=1.0'
else:
accuracyStr = 'accuracy={:.2f}'.format(accuracy)
# f1 score
f1 = metrics.f1_score(actual_labels, y_classification)
# composite metric
if p_value_logbase <= 0: # by AUC first (but not only, so we also find best P value)
CM = AUC + (-math.log(pValue,1e300)) / 1000 # de-emphasize log(P) more by /1000 because e.g. log1e300(1e-20) = -0.06666 may still have an effect on "AUC-only" CM (AUC+logP, not AUC*logP)
CM /= 100 # reduce further because above log may still return values large enough to affect the AUC-only CM for very small Ps (but we still want a little effect of P on CM so we can find the best AUC first, then the best P)
CMstr = str(CM) # placeholder, not useful
else:
CM, CMstr = calcCM(pValue, AUC, accuracy, p_value_logbase=p_value_logbase) # Higher values of p_value_logbase emphasize AUC over P value; p_value_logbase=0 uses AUC only
#printSB('nuSVC_CM: ',CMstr)
if CM > CMbest:
CMbest = CM
CM_logbase10, CM_logbase10str = calcCM(pValue, AUC, accuracy, p_value_logbase=10) # for scoring/optimization purposes we may the AUC-only CM, but for reporting we always return the CM_logbase10 so it's comparable to other runs
dict_best = {'nu':nu, 'pValue':pValue, 'Pstr':formatPstring(pValue), 'AUC':AUC, 'AUCstr':AUCstr, 'accuracy':accuracy, 'accuracyStr':accuracyStr, 'f1':f1, 'CM':CM, 'CMstr':CMstr, 'CM_logbase10':CM_logbase10, 'CM_logbase10str':CM_logbase10str, 'decision_function_df':decision_function_df, 'y_pred':y_classification}
except:
if suppressWarnings == False: printSB('nuSVC_CM exception at nu = {:.2f}'.format(nu) + ': method terminating early (this may be OK)')
# sort by CM
#sorted_dList = sorted(dict_out_list, key=lambda x: x['CM'], reverse=True) # sort by CM: sorted_dList[0] is the best dict
return dict_best # return the dict with best results (see above for keys). get elements of dict like this: d_best.get('decision_function_df')
### computes nu-support vector classification on X & Y, selects the best hyperparams based on CM (composite metric)
### pass subject-wise mean PCs in PCmeans_df: 1st col is name, 2nd col is classID, rest are PCs
### needs global vars: maxIter_svm, param_grid for grid search, kStandardScaleResult
# 2024-04-23: CM now includes accuracy as well as log(P) and AUC
def SVC_CM(PCmeans_df, param_grid, n_splits=10, n_repeats=5, suppressWarnings=False, defaultLabel='PC', p_value_logbase=10):
# PCmeans_df can be UMAP means, etc
from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV, RepeatedKFold
from sklearn import metrics
from sklearn.preprocessing import StandardScaler
rkf = RepeatedKFold(n_splits=min(n_splits,PCmeans_df.shape[0]), n_repeats=n_repeats)
scaler = StandardScaler()
dict_best = {'CM':-1} # in case no solution/exception
first_PC_index = next((i for i, col in enumerate(PCmeans_df.columns) if col.startswith(defaultLabel)), -1) # find index of 1st PCx or UMAPx etc column
if first_PC_index == -1:
printSB('SVC_CM: PCx columns could not be found (fatal error)')
return dict_best, pd.DataFrame()
X = PCmeans_df.iloc[:,first_PC_index:]
Y = PCmeans_df['classID']
svm_model = SVC()
try: # some values of nu are illegal and will throw an exception: catch it and return what we have
grid_search = GridSearchCV(estimator=svm_model, param_grid=param_grid, cv=rkf)
grid_search.fit(X, Y)
y_pred = grid_search.predict(X) # make predictions on the input set to check accurcy of linear separability
decision_function = grid_search.decision_function(X)
if kStandardScaleResult:
decision_function = scaler.fit_transform(decision_function.reshape(-1, 1)) # often decision_function contains tiny values so rescaling standardizes when averaging with different hyperparameters
decision_function_df = pd.DataFrame(data=decision_function,columns=['score'])
#printSB('nuSVC_CM: about to stats.ttest_ind')
# now compute t-test for dec fn, AUC -> CM
a = decision_function_df.loc[PCmeans_df['classID'] == 0]['score']
b = decision_function_df.loc[PCmeans_df['classID'] == 1]['score']
_, pValue = stats.ttest_ind(a, b, equal_var = True, alternative='less') #run independent 2 sample T-Test
#printSB('nuSVC_CM: about to metrics.roc_auc_score')
# AUC
#predicted_probabilities = decision_function_df['score'] THIS IS WRONG: because dec_func varies from <0 to >0 (either side of the hyperplane) it is not a 0-1 probability which is what metrics.roc_auc_score expects
y_classification = np.where(decision_function_df['score'] > 0, 1, 0) # so convert it to a predicted 0 or 1 class, depending on <0 vs >0
actual_labels = PCmeans_df['classID']
AUC = metrics.roc_auc_score(actual_labels, y_classification)
if AUC == 1.0:
AUCstr = 'AUC = 1.0'
else:
AUCstr = 'AUC = {:.2f}'.format(AUC)
#printSB('nuSVC_CM: about to metrics.accuracy_score')
# accuracy
accuracy = metrics.accuracy_score(actual_labels, y_classification)
if accuracy == 1.0:
accuracyStr = 'accuracy=1.0'
else:
accuracyStr = 'accuracy={:.2f}'.format(accuracy)
# f1 score
f1 = metrics.f1_score(actual_labels, y_classification)
# composite metric
CM, CMstr = calcCM(pValue, AUC, accuracy, p_value_logbase=p_value_logbase) # Higher values of p_value_logbase emphasize AUC over P value
dict_best = {'C':grid_search.best_params_['C'], 'gamma':grid_search.best_params_['gamma'], 'kernel':grid_search.best_params_['kernel'], 'pValue':pValue, 'Pstr':formatPstring(pValue), 'AUC':AUC, 'AUCstr':AUCstr, 'accuracy':accuracy, 'accuracyStr':accuracyStr, 'f1':f1, 'CM':CM, 'CMstr':CMstr, 'decision_function_df':decision_function_df, 'y_pred':y_classification}
except:
if suppressWarnings == False: printSB('SVC_CM exception')
# sort by CM
#sorted_dList = sorted(dict_out_list, key=lambda x: x['CM'], reverse=True) # sort by CM: sorted_dList[0] is the best dict
return dict_best # return the dict with best results (see above for keys). get elements of dict like this: d_best.get('decision_function_df')
# scikit-learn.org/stable/modules/generated/sklearn.discriminant_analysis.QuadraticDiscriminantAnalysis.html
def QDA(PCmeans_df, p_value_logbase=10):
from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis
from sklearn import metrics
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X = PCmeans_df.iloc[:,2:]
Y = PCmeans_df.iloc[:,1]
clf = QuadraticDiscriminantAnalysis()
clf.fit(X, Y)
y_pred = clf.predict(X)
class1probs = clf.predict_proba(X)[:,1] # probabilities of class1 membership https://scikit-learn.org/stable/modules/generated/sklearn.discriminant_analysis.QuadraticDiscriminantAnalysis.html#sklearn.discriminant_analysis.QuadraticDiscriminantAnalysis.predict_proba
AUC = metrics.roc_auc_score(Y, y_pred)
if AUC == 1.0:
AUCstr = 'AUC = 1.0'
else:
AUCstr = 'AUC = {:.2f}'.format(AUC)
#accuracy = metrics.accuracy_score(y, y_pred)
accuracy = clf.score(X, Y)
if accuracy == 1.0:
accuracyStr = 'accuracy=1.0'
else:
accuracyStr = 'accuracy={:.2f}'.format(accuracy)
y_classification = np.where(class1probs > 0.5, 1, 0) # class0 if probability <0.5, else class1
probability_accuracy = metrics.accuracy_score(Y, y_classification)
if probability_accuracy == 1.0:
probability_accuracyStr = 'accuracy=1.0'
else:
probability_accuracyStr = 'accuracy={:.2f}'.format(probability_accuracy)
decision_function = clf.decision_function(X)
decision_function = scaler.fit_transform(decision_function.reshape(-1, 1)) # often decision_function contains tiny values so rescaling standardizes when averaging with different hyperparameters
decision_function_df = pd.DataFrame(data=decision_function,columns=['score'])
# now compute t-test for dec fn, AUC -> CM
a = decision_function_df.loc[PCmeans_df['classID'] == 0]['score']
b = decision_function_df.loc[PCmeans_df['classID'] == 1]['score']
_, pValue = stats.ttest_ind(a, b, equal_var = True, alternative='less') #run independent 2 sample T-Test
# composite metric
CM, CMstr = calcCM(pValue, AUC, accuracy, p_value_logbase=p_value_logbase) # Higher values of p_value_logbase emphasize AUC over P value
# f1 score
f1 = metrics.f1_score(Y, y_pred)
dict_out = {'pValue':pValue, 'Pstr':formatPstring(pValue), 'AUC':AUC, 'AUCstr':AUCstr, 'accuracy':accuracy, 'accuracyStr':accuracyStr, 'f1':f1, 'CM':CM, 'CMstr':CMstr, 'decision_function_df':decision_function_df, 'y_pred':y_pred, 'class1probs':class1probs, 'probability_accuracy':probability_accuracy, 'probability_accuracyStr':probability_accuracyStr}
return dict_out
# pass a multivariate ie spectra over several T points as in a Duetta/Mk5 PC series
# will average the spectra at the chosen T points and subract this average from all spectra, returning a new bsub'd dataframe (df_multivariate_ALLcols_bsub)
# also returns the mean baseline spectrum used for the bsub for each name (df_baselineSpectra)
# df_baselineSourceSpectra are spectral T series from which averages will be calculated, typically df_multivariate but with ALL T points included (ie Tend=-1 equivalent)
def baselineSubtract(df_multivariate, df_baselineSourceSpectra, Tavg_startIx, Tavg_endIx): # Tpoints Tavg_startIx, Tavg_endIx (0-based, inclusive) used to compute mean spectrum for baseine subtraction
namesUnique = df_multivariate['name'].unique() # unique names from 'name' col
colNames = df_multivariate.columns # original colnames, ALL cols
for j in range(namesUnique.shape[0]):
#for j in range(1):
name = namesUnique[j]
# printSB(j,name)
df_baselineSourceSpectra_byname = df_baselineSourceSpectra[df_baselineSourceSpectra['name']==name] # extract all rows at all Ts by name from the ALL-in df
df_baselineSourceSpectra_byname_spectra = df_baselineSourceSpectra_byname.iloc[:,spectrum1stColIndex:] # just the spectral cols
baselineRows = df_baselineSourceSpectra_byname.iloc[Tavg_startIx:Tavg_endIx+1,] # the rows that will be averaged and used for subtraction
baselineRows.reset_index(inplace=True,drop=True)
baselineRows_spectra = baselineRows.iloc[:,spectrum1stColIndex:] # just the spectral cols
baselineMean = baselineRows_spectra.mean(axis=0)
baselineMean_df = pd.DataFrame(baselineMean).T # make a df and transpose from 730 rows,1 col to 1 row,730 cols
df_multivariate_byname = df_multivariate[df_multivariate['name']==name] # extract all rows at all Ts by name
df_multivariate_byname_spectra = df_multivariate_byname.iloc[:,spectrum1stColIndex:] # just the spectral cols
df_multivariate_byname_spectra_bsub = df_multivariate_byname_spectra.copy() # this will be the bsub'd table
df_multivariate_byname_spectra_bsub -= baselineMean # bsub'd spectral rows for name at all Ts
metaCols_byname = df_multivariate[df_multivariate['name']==name].iloc[:,0:spectrum1stColIndex] # the 1st 3 meta cols (name, classID, T) at all Ts for name
df_multivariate_byname_bsub = pd.concat([metaCols_byname,df_multivariate_byname_spectra_bsub], axis=1, ignore_index=True) # concat (by columns) the metaColumns and the bsub'd spectral columns
df_multivariate_byname_bsub.columns = colNames # essential else concats fail below
metaCols_byname_1stRow = metaCols_byname.iloc[0:1] # just 1st row, for building df_baselineSpectra
metaCols_byname_1stRow.reset_index(inplace=True,drop=True)
baselineSpectrum_byname = pd.concat([metaCols_byname_1stRow,baselineMean_df], axis=1)
baselineSpectrum_byname.columns = colNames
if j == 0:
df_multivariate_ALLcols_bsub = pd.DataFrame(data=df_multivariate_byname_bsub)
df_multivariate_ALLcols_bsub.columns = colNames # essential else concats fail below
df_baselineSpectra = pd.DataFrame(data=baselineSpectrum_byname)
df_baselineSpectra.columns = colNames # essential else concats fail below
else:
df_multivariate_ALLcols_bsub = pd.concat([df_multivariate_ALLcols_bsub, df_multivariate_byname_bsub], axis=0, ignore_index=True) # concat (by rows) into the master bsub'd df
df_baselineSpectra = pd.concat([df_baselineSpectra, baselineSpectrum_byname], axis=0, ignore_index=True)
df_multivariate_ALLcols_bsub.reset_index(inplace=True,drop=True)
df_baselineSpectra.reset_index(inplace=True,drop=True)
return df_multivariate_ALLcols_bsub, df_baselineSpectra
def generate_balanced_random_0sand1s(N):
import random
if N % 2 != 0: N+= 1
half_N = N // 2
zeros = ones = half_N
result = [0] * zeros + [1] * ones
random.shuffle(result)
return result
def printNamesAndClassIDs(df, prependLF=False, appendLF=False):
if prependLF: printSB()
#printSB(df.drop_duplicates(subset=['name', 'classID']).copy().iloc[:,0:2].to_string(index=False))
# Group by 'name' and 'classID', count the occurrences, and reset the index
result_df = df.groupby(['name', 'classID']).size().reset_index(name='count')
# Rename the columns for clarity
result_df = result_df.rename(columns={'name': 'Name', 'classID': 'ClassID', 'count': 'Count'})
# Print the result in a nicely formatted table
printSB(result_df.to_string(index=False))
if appendLF: printSB()
return
def randomizeImageData(df, unique_names_classIDs, randflag=1, printNamesAndIDs=True, verbose=1):
# randomizes df classIDs by name (ie multiple rows with same name will be randomized together to same new classID), returns a copy of df and the appropriate randSuffix
# randflag selects 1 of 3 modes (see code)
# requires: unique_names_classIDs, typically: imageData.drop_duplicates(subset=['name', 'classID']).copy().iloc[:,0:2] # unique names and matching classID cols
# set verbose=0 to suppress all printing
# all columns in df_rand will be the same as df, except the 'classID' col will be shuffled
randSuffix = ''
df_rand = df.copy() # if randflang=0 just returns a copy of df
if randflag == 1: # rand using generate_balanced_random_0sand1s()
if verbose > 0: printSB('Randomizing classIDs using generate_balanced_random_0sand1s...')
randSuffix = ' RAND'
randClassList = generate_balanced_random_0sand1s(unique_names_classIDs.shape[0])
for j, name in enumerate(unique_names_classIDs['name']):
randClassID = randClassList[j] # random 0 or 1 sequence (balanced)
df_rand.loc[df['name'] == name,['classID']] = randClassID # assign classID = randClassID to all rows whose name = name
elif randflag == 2: # mod 2 classID
if verbose > 0: printSB('Randomizing classIDs using mod 2...')
randSuffix = ' RAND'
unique_names_classIDs_RAND = unique_names_classIDs.sample(frac=1).reset_index(drop=True) # shuffle the rows of unique_names_classIDs into a copy
unique_names_classIDs_RAND.sort_values(by=['classID'],inplace=True,ignore_index=True) # sort by classID only: names will remain shuffled WITHIN a class, eliminating the same randomization for each pass mod 2
for j, name in enumerate(unique_names_classIDs_RAND['name']):
randClassID = j % 2 # mod 2 classID for name
df_rand.loc[df['name'] == name,['classID']] = randClassID # assign classID = randClassID to all rows whose name = name
elif randflag == 3: # mod 2 classID but names sorted and NOT randomized before classID sort
if verbose > 0: printSB('Randomizing classIDs using mod 2...')
randSuffix = ' RAND'
unique_names_classIDs_RAND = unique_names_classIDs.copy().reset_index(drop=True) # shuffle the rows of unique_names_classIDs into a copy
unique_names_classIDs_RAND.sort_values(by=['classID','name'],inplace=True,ignore_index=True) # sort by classID first (so mod 2 re-assignment will produce a balanced RAND), then by name so same randomization will occur in contrast to randflag = 2
for j, name in enumerate(unique_names_classIDs_RAND['name']):
randClassID = j % 2 # mod 2 classID for name
df_rand.loc[df['name'] == name,['classID']] = randClassID # assign classID = randClassID to all rows whose name = name
if printNamesAndIDs and (verbose > 0): printNamesAndClassIDs(df_rand, True, True)
return df_rand, randSuffix # df was randomized as a copy
def bsub(raw_df):
# detects spectra with classID = -1, averages all these as a unvariate set
# subtracts this mean dark spectrum from all other rows
# build NON-norm univariate vectors
# build a 2D univariate features array of shape (n_instances, nL*nT)
names_unique = raw_df['name'].unique()
n_subjects = names_unique.shape[0]
nT_original = raw_df['T'].unique().shape[0]
n_lambdas_raw = raw_df.columns[spectrum1stColIndex:].shape[0] # lambdas from raw input data before range pruning etc
spectralData_NONnorm_univariate = np.zeros((n_subjects, n_lambdas_raw*nT_original)) # 2D np.ndarray,
#names_unique
names = []
rootnames = []
classIDs = []
for instanceCtr, name in enumerate(names_unique): # counts subjects
rowsByName = spectralData_raw.loc[spectralData_raw['name']==name]
rowsByName.reset_index(inplace=True, drop=True) # all rows with name = name
classID = rowsByName['classID'].iloc[0]
names.append(name)
rootnames.append(extractRootname(name)) # Split the string at the last occurrence of sep, and return a 3-tuple containing the part before the separator, the separator itself, and the part after the separator (https://docs.python.org/2/library/stdtypes.html#str.rpartition)
classIDs.append(classID)
spectralValues_NONnorm_allPerSubjectList = [] # we will flatten all spectra from all T points into a single row for each subject
for Tctr in range(nT_original):
rowByT = rowsByName.iloc[Tctr]
spectralValuesAtT = rowByT[spectrum1stColIndex:spectrum1stColIndex+n_lambdas_raw].to_numpy() # extracted spectrum at time Tctr
spectralValues_NONnorm_allPerSubjectList.append(spectralValuesAtT)
spectralValues_NONnorm_allPerSubject_np = np.asarray(spectralValues_NONnorm_allPerSubjectList).flatten()
spectralData_NONnorm_univariate[instanceCtr,:] = spectralValues_NONnorm_allPerSubject_np
### spectralData_cond_univariate now contains the univariate 2D array: dim 1=subjects (incl replicates); dim 2=dummy; dim 3=lambda bins
# assemble into a df
spectrum1stColIndex_univariate = 2
names_df = pd.DataFrame(data=names,columns=['name'])
classIDs_df = pd.DataFrame(data=classIDs,columns=['classID'])
spectralData_NONnorm_univariate_df = pd.DataFrame(data=spectralData_NONnorm_univariate)
spectralData_NONnorm_univariate_df = pd.concat([names_df,classIDs_df,spectralData_NONnorm_univariate_df],axis=1)
# extract all the univariate dark spectra ie classID=-1
dark_df = spectralData_NONnorm_univariate_df.loc[spectralData_NONnorm_univariate_df['classID'] == -1]
# extract spectral columns
dark_spectra_df = dark_df.iloc[:,2:]
# average all the dark rows
dark_mean_df = dark_spectra_df.mean()
# extract all the univariate NONdark spectra ie classID>-1 as copies
NONdark_df = spectralData_NONnorm_univariate_df.loc[spectralData_NONnorm_univariate_df['classID'] > -1].copy()
NONdark_df.reset_index(inplace=True, drop=True)
NONdark_spectra_df = NONdark_df.iloc[:,2:] # spectral cols only
NONdark_bsub_spectra_df = NONdark_spectra_df - dark_mean_df # dark subtract
NONdark_bsub_spectra_df[NONdark_bsub_spectra_df < 0] = 0 # fix neg values
NONdark_df = pd.concat([NONdark_df.iloc[:,0:2],NONdark_bsub_spectra_df],axis=1) # reassemble with meta columns
#NONdark_df.columns = cols_orig # restore col names ie wavelengths
# reshape into a multivariate df by T
raw_multivariate_df = raw_df.iloc[0:1,:].copy() # copy 1st row of spectralData_raw to get the corrrect cols and their names: we'll replace the values from bsub'd df and append new ones
raw_multivariate_1stRow = raw_multivariate_df.copy()
nT_original = raw_df['T'].unique().shape[0]
n_lambdas_raw = raw_df.columns[spectrum1stColIndex:].shape[0] # no. lambdas from raw input data before range pruning etc
for rowCtr in range(NONdark_df.shape[0]):
row_univariate = NONdark_df.iloc[rowCtr,:] # next entire univariate row for next subject
name = row_univariate.iloc[0]
classID = row_univariate.iloc[1]
# printSB(name)
for tCtr in range(nT_original):
if (rowCtr==0) and (tCtr==0): # 1st multivariate row is special
raw_multivariate_df.iloc[0,0:1] = name
raw_multivariate_df.iloc[0,1:2] = classID
raw_multivariate_df.iloc[0,2:3] = tCtr
raw_multivariate_df.iloc[0,3:] = row_univariate[2+tCtr*n_lambdas_raw:2+(tCtr+1)*n_lambdas_raw]
else:
raw_multivariate_newRow = raw_multivariate_1stRow.copy()
raw_multivariate_newRow.iloc[0,0:1] = name
raw_multivariate_newRow.iloc[0,1:2] = classID
raw_multivariate_newRow.iloc[0,2:3] = tCtr
raw_multivariate_newRow.iloc[0,3:] = row_univariate[2+tCtr*n_lambdas_raw:2+(tCtr+1)*n_lambdas_raw]
raw_multivariate_df = pd.concat([raw_multivariate_df,raw_multivariate_newRow], axis=0, ignore_index=True) # append the new bsub'd multivariate row
return raw_multivariate_df
# load multiple CSVs into a single df. Typically get csvList from get_child_files()
def loadCSVs(csvList):
for n, csv in enumerate(csvList):
# printSB(n,csv)
if n == 0:
# read 1st line: do we have a comment?
with open(csv,'r') as f:
csvLine1 = f.readline() # csvLine1 might be a comment
if csvLine1[:1] != "'": # must match app.kPandasCSVcommentChar in ITRK and SB
csvLine1 = '' # if not a comment we don't want to dbl-write line1 below
printSB('Loading: ' + os.path.basename(csv) + '...')
if n == 0:
df = pd.read_csv(csv, comment="'") # skip lines with '
columns0 = df.columns
else: # append
dfN = pd.read_csv(csv, comment="'")
if (columns0 != dfN.columns).any(): # if any of the columns differ
printSB('*** WARNING: columns for ' + os.path.basename(csv) + ' differ from 1st csv: resetting to columns of 1st csv. Make sure this is correct!')
dfN.columns = columns0
df = pd.concat([df,dfN],axis=0,ignore_index=True)
return df, csvLine1 # return the combined df and the csv line1 comment, if any
def loadCSVs_v2(csvList, startChar='L', stripLeadingLs=True, printLoading=True):
# v2 autodetects 'L' lambda columns
for n, csv in enumerate(csvList):
# printSB(n,csv)
if n == 0:
# read 1st line: do we have a comment?
with open(csv,'r') as f:
csvLine1 = f.readline() # csvLine1 might be a comment
if csvLine1[:1] != "'": # must match app.kPandasCSVcommentChar in ITRK and SB
csvLine1 = '' # if not a comment we don't want to dbl-write line1 below
if printLoading: printSB('Loading: ' + os.path.basename(csv) + '...')
if n == 0:
df = pd.read_csv(csv, comment="'") # skip lines with '
columns0 = df.columns
else: # append
dfN = pd.read_csv(csv, comment="'")
if (columns0 != dfN.columns).any(): # if any of the columns differ
printSB('*** WARNING: columns for ' + os.path.basename(csv) + ' differ from 1st csv: resetting to columns of 1st csv. Make sure this is correct!')
dfN.columns = columns0
df = pd.concat([df,dfN],axis=0,ignore_index=True)
firstLambdaCol, nLambdas, lambdaList = get_column_indices(df)
if nLambdas == 0: # 'L' cols were not detected
firstLambdaCol, nLambdas = -1, -1
elif stripLeadingLs: # 'L' cols were detected, strip the leading 'L's from their names
df.columns = df.columns.str.replace('^' + startChar, '', regex=True)
return df, csvLine1, firstLambdaCol, nLambdas, lambdaList # return the combined df and the csv line1 comment, if any
def loadCSVs_v3(csvList, startChar='L', stripLeadingLs=True, printLoading=True):
# v2 autodetects 'L' lambda columns
# v3: returns raw col names for reference
for n, csv in enumerate(csvList):
# printSB(n,csv)
if n == 0:
# read 1st line: do we have a comment?
with open(csv,'r') as f:
csvLine1 = f.readline() # csvLine1 might be a comment
if csvLine1[:1] != "'": # must match app.kPandasCSVcommentChar in ITRK and SB
csvLine1 = '' # if not a comment we don't want to dbl-write line1 below
if printLoading: printSB('Loading: ' + os.path.basename(csv) + '...')
if n == 0:
df = pd.read_csv(csv, comment="'") # skip lines with '
columns0 = df.columns
else: # append
dfN = pd.read_csv(csv, comment="'")
if (columns0 != dfN.columns).any(): # if any of the columns differ
printSB('*** WARNING: columns for ' + os.path.basename(csv) + ' differ from 1st csv: resetting to columns of 1st csv. Make sure this is correct!')
dfN.columns = columns0
df = pd.concat([df,dfN],axis=0,ignore_index=True)
firstLambdaCol, nLambdas, lambdaList = get_column_indices(df)
if nLambdas == 0: # 'L' cols were not detected
firstLambdaCol, nLambdas = -1, -1
elif stripLeadingLs: # 'L' cols were detected, strip the leading 'L's from their names
df.columns = df.columns.str.replace('^' + startChar, '', regex=True)
return df, csvLine1, firstLambdaCol, nLambdas, lambdaList, columns0 # return the combined df and the csv line1 comment, if any. Also return the raw col names
# main csv loader function
def readCSVs(csvIN, dirIN, dirOUT, nowStr, classList, intensityLO, intensityHI, lambdaSTART, lambdaEND, replaceNameUnderscoresWithDashes, subSampleSize, augX=1, T=[0]):
# delete all rows whose classID is not in classList (pass None or empty classList to accept all classes)
# if csv has a 'T' column, include only rows where T=T; pass -1 to include all T columns;can pass T=[startIx,endIx] with inclusive 0-based indexes of T columns to retain
# if augX < 0 uses augPow2df augmentation
if dirOUT=='':
parentFolder = os.path.dirname(csvIN)
else:
parentFolder = dirOUT
if dirIN == '':
fBaseExt = os.path.basename(csvIN) # includes .csv extension
fBaseNoExt = Path(csvIN).stem # w/o extension
else:
fBaseExt = 'MULTIPLE_FROM:' + Path(dirIN).stem
fBaseNoExt = fBaseExt
if nowStr == '':
datetimeSuffix = ''
else:
datetimeSuffix = '_' + nowStr
#printSB("Loading " + csvIN + "...")
csvList = get_child_files(dirIN, csvIN, 'csv')
imageData_df, csvLine1, firstLambdaCol, nLambdas, _ = loadCSVs_v2(csvList) # v2 will autodetect new 'L' lambda columns, if not firstLambdaCol, nLambdas returned as -1, so regular code below will do the old auto-detection
if len(csvList) == 1:
printSB("1 CSV LOADED. Total {:,} rows".format(imageData_df.shape[0]))
else:
printSB(str(len(csvList)) + " CSVs LOADED. Total {:,} rows".format(imageData_df.shape[0]))
headerRow, firstLambdaCol, nLambdas = fixOrangeHeaders(imageData_df)
for h in headerRow:
if ('unique4' in h):
imageData_df.rename(columns = {h:'unique4'}, inplace = True)
break
allClasses_original = imageData_df.classID.unique() # list of all classes in the original csv's before filtering
#imageData_df.loc[imageData_df['classID'] == 2, 'classID'] = 1 # you can reset one class to another eg all class2 rows reset to class1 (say you want to combine class1&2 in the analysis)
ignoreClassList = (classList is None) or (len(classList)==0)
if ignoreClassList == False:
imageData_df = imageData_df[imageData_df['classID'].isin(classList)] # filter only rows whose classIDs are in the classPair list
nClasses = imageData_df.classID.unique().shape[0]
hasTcol = ('T' in headerRow)
if hasTcol:
if (len(T)==1) and (T[0] > -1): # is T a single integer?
imageData_df = imageData_df[imageData_df['T'] == T[0]] # filter rows by a single T
elif len(T)==2: # list of start/end indexes?
# Filter the DataFrame to keep rows where T is within the inclusive range
imageData_df = imageData_df[(imageData_df['T'] >= T[0]) & (imageData_df['T'] <= T[1])]
else:
raise ValueError("Illegal T param in readCSVs method")
adjustColumnTypes(imageData_df, firstLambdaCol, nLambdas, cleanLambdaColHeaders=True, promoteLambdaColsToFloat64=True) # reduce memory footprint
# restrict by intensity, etc
#imageData_df = imageData_df.loc[(imageData_df['intensity'] > 70)]
if 'intensity' in imageData_df.columns: # mk5/6 csv's may not have intensity col
imageData_df = imageData_df.loc[(imageData_df['intensity'] >= intensityLO) & (imageData_df['intensity'] < intensityHI)]
imageData_df.reset_index(inplace=True, drop=True)
nLambdas, imageData_df = restrictLambdaRange(imageData_df, lambdaSTART, lambdaEND, firstLambdaCol, nLambdas)
printSB('firstLambdaCol: ' + str(firstLambdaCol))
printSB('nLambdas: ' + str(nLambdas))
# Extract lambda headers
lambda_headers = imageData_df.columns[firstLambdaCol:firstLambdaCol + nLambdas]
# Convert headers to floats
lambdas_np = np.array(lambda_headers, dtype=float)
rootNames, rootNames_df = extractRootnames(imageData_df, '_AUG_') # in case need pre-augmented subject rootnames
if replaceNameUnderscoresWithDashes:
imageData_df['name'] = imageData_df['name'].str.replace('_','-')
# subsample for debugging
if subSampleSize > 0:
subSampleSize = min(subSampleSize,imageData_df.shape[0])
imageData_df = imageData_df.sample(n=subSampleSize, ignore_index=True, random_state=123) # ignore_index=True is essential
imageData_df.sort_values(by=['classID','name','unique4'],inplace=True,ignore_index=True) # for human user, not necessary for embedding. ignore_index=True is critical else you end up with a randomized frame somehow
if augX < 0:
augPow2 = -augX
printSB(f'Augmenting data by pow2X ({2**augPow2} total uniques)...')
imageData_df, _ = augPow2(dfIN, augPow2) # pow2 data augmentation
elif augX > 1:
printSB('Augmenting data ' + str(augX) + 'X...')
imageData_df, _ = augXdf(imageData_df, augX) # data augmentation
elif augX == 0: # special case for Mark6
printSB('Augmenting each kernel to a unique instance...')
imageData_df = aug0df(imageData_df) # per-kernel data augmentation
names_unique = imageData_df['name'].unique()
N_names_unique = names_unique.shape[0]
unique_names_classIDs = imageData_df.drop_duplicates(subset=['name', 'classID']).copy().iloc[:,0:2] # unique names and matching classID cols
unique_names_classIDs.reset_index(inplace=True, drop=True)
d = {}
d['imageData'] = imageData_df
d['fBaseExt'] = fBaseExt
d['fBaseNoExt'] = fBaseNoExt
d['parentFolder'] = parentFolder
d['datetimeSuffix'] = datetimeSuffix
d['csvLine1'] = csvLine1
d['firstLambdaCol'] = firstLambdaCol
d['nLambdas'] = nLambdas
d['lambdas_np'] = lambdas_np # lambdas an a np array
d['nClasses'] = nClasses # after filtering
d['allClasses_original'] = allClasses_original # before filtering
d['rootNames'] = rootNames
d['names_unique'] = names_unique
d['N_names_unique'] = N_names_unique
d['unique_names_classIDs'] = unique_names_classIDs
return d # return results as a dict
# merge multiple CSVs in a folder and write to a new csv in same folder
def mergeCSVs(directory_path, dirOut=''):
"""
Combines all CSV files within a directory into a single DataFrame, then saves it.
Args:
directory_path (str): The path to the directory containing the CSV files.
"""
all_files = os.listdir(directory_path)
csv_files = [file for file in all_files if file.endswith('.csv')]
dataframes = []
for n, file in enumerate(csv_files):
filepath = os.path.join(directory_path, file)
printSB('Loading: ' + os.path.basename(file) + '...')
df = pd.read_csv(filepath)
if n == 0:
columns0 = df.columns
else:
if (columns0 != df.columns).any(): printSB('*** WARNING: columns for ' + os.path.basename(file) + ' differ from 1st csv: resetting to columns of 1st csv. Make sure this is correct!') # if any of the columns differ
dataframes.append(df)
combined_df = pd.concat(dataframes, ignore_index=True)
output_filename = f"{len(csv_files)} merged.csv"
if dirOut == '':
output_filepath = os.path.join(directory_path, output_filename)
else:
output_filepath = os.path.join(dirOut, output_filename)
printSB('Saving merged CSV to: ' + output_filename + '...')
combined_df.to_csv(output_filepath, index=False)
return output_filepath
# merge multiple CSVs in a list and write to a new csv in same folder
def mergeCSVsFromList(csvList, dirOut):
"""
Combines all CSV files in csvList into a single DataFrame, then saves it.
Args:
listOfCSVs (str list):
"""
imageData, csvLine1, firstLambdaCol, nLambdas, _ = loadCSVs_v2(csvList)
output_filename = f"{len(csvList)} merged.csv"
output_filepath = os.path.join(dirOut, output_filename)
printSB('\nSaving merged CSV to: ' + output_filename + ' (' + str(imageData.shape[0]) + ' rows)...')
imageData.to_csv(output_filepath, index=False)
return output_filepath
def Raman_blank(spectra_df, spectrum1stColIndex, raman_nm_start, raman_nm_end):
# blanks the Raman water spectrum by linearly interpolating all spectra in spectra_df between raman_nm_start & raman_nm_end nm
# spectrum1stColIndex: index of 1st column of spectral data after all meta columns
if (raman_nm_start<=0) or (raman_nm_end<=0): return spectra_df # set raman_nm_start or raman_nm_end to 0 to omit blanking
nmArr = spectra_df.columns[spectrum1stColIndex:].astype('float') # wavelengths
metaCols = spectra_df.iloc[:,:spectrum1stColIndex] # non-spectral meta cols only
df = spectra_df.iloc[:,spectrum1stColIndex:] # spectral cols only
df.reset_index(inplace=True, drop=True) # don't forget this if you want to re-index rows from 0
# find the index of the column @ raman_nm_start
for j in range(nmArr.shape[0]):
if nmArr[j] >= raman_nm_start:
raman_nm_start_ix = j
break
# find the index of the column @ raman_nm_end
for j in range(nmArr.shape[0]-1, -1, -1):
if nmArr[j] <= raman_nm_end:
raman_nm_end_ix = j
break
new_df = df.copy() # make a copy, we'll replace only the columns within the Raman band
nmRaman = nmArr[raman_nm_start_ix:raman_nm_end_ix] # nm values within the Raman band
raman_start_y = new_df.iloc[:,raman_nm_start_ix] # y-values from each row @ raman_nm_start
raman_end_y = new_df.iloc[:,raman_nm_end_ix] # y-values from each row @ raman_nm_end
raman_start_y.reset_index(inplace=True, drop=True) # don't forget this if you want to re-index rows from 0
raman_end_y.reset_index(inplace=True, drop=True) # don't forget this if you want to re-index rows from 0
x = np.array([nmRaman[0], nmRaman[-1]]) # 2 x-values (nm) defining the start & end of the Raman band
for rowCtr in range(new_df.shape[0]):
y = np.array([raman_start_y[rowCtr], raman_end_y[rowCtr]]) # x,y are 2 points that define the line thru the Raman band
ynew = np.interp(nmRaman, x, y) # interpolated y-values thru the Raman band
new_df.iloc[rowCtr,raman_nm_start_ix:raman_nm_end_ix] = ynew # replace the Raman spectrum with the interpolated values
# reassemble with original meta cols
spectra_RamanBlanked_df = pd.concat([metaCols,new_df], axis=1)
return spectra_RamanBlanked_df
# make sure dirOUT exists
import os
try:
if not os.path.exists(dirOUT): os.makedirs(dirOUT)
except:
pass
# adjust column types in place to reduce memory use
def adjustColumnTypes(df, firstLambdaColIx, nLambdaCols, cleanLambdaColHeaders=True, promoteLambdaColsToFloat64=False): # typically pass imageData in df
# promoteLambdaColsToFloat64=True if averaging to maintain highest precision, typically pass True in all avgSpectra scripts
if 'name' in df.columns: df['name'] = df['name'].astype('string') # force to str type in case names happen to be numeric
if 'classID' in df.columns: df['classID'] = df['classID'].astype('int8')
if 'T' in df.columns: df['T'] = df['T'].astype('uint8')
if 'X' in df.columns: df['X'] = df['X'].astype('uint16')
if 'Y' in df.columns: df['Y'] = df['Y'].astype('uint16')
if 'Z' in df.columns: df['Z'] = df['Z'].astype('uint8')
if 'unique4' in df.columns: df['unique4'] = df['unique4'].astype('uint16')
if 'intensity' in df.columns: df['intensity'] = df['intensity'].astype('float32')
# Convert lambda columns to float32 or 64
columns_to_convert = df.columns[firstLambdaColIx: firstLambdaColIx + nLambdaCols]
for col in columns_to_convert:
if promoteLambdaColsToFloat64:
df[col] = df[col].astype(np.float64)
else:
df[col] = df[col].astype(np.float32)
if cleanLambdaColHeaders:
df.rename(columns={col: float(col[2:]) if isinstance(col, str) and col.startswith('C#') else col for col in df.columns}, inplace=True) # strip all leading 'C#' substrings in col headers (check if col is already numeric to avoid AttributeError: 'float' object has no attribute 'startswith')
return
# restrict range of wavelengths by dropping wavelength columns from imageData_df in place
# imageData_df must be pre-processed with call to adjustColumnTypes()
# pass lambda_start = -1 to omit any modification of imageData_df
# returns the reduced nLambdas, and the truncated dataframe (all cols other than lambda cols are returned as is)
def restrictLambdaRange_OLD(imageData_df, lambda_start, lambda_end, firstLambdaCol, nLambdas):
if lambda_start < 0: return nLambdas, imageData_df
lambda_columns = imageData_df.columns[firstLambdaCol:firstLambdaCol+nLambdas] # str list of lambda headers
# Filter columns based on the numeric range
columns_to_keep = [col for col in imageData_df.columns if (col not in lambda_columns) or (lambda_start<=float(col)<=lambda_end)]
# Filter columns based on the numeric range
nColsRemoved = imageData_df.shape[1] - len(columns_to_keep)
nLambdas_new = nLambdas - nColsRemoved
filtered_df = imageData_df[columns_to_keep]
return nLambdas_new, filtered_df
# restrict range of wavelengths by dropping wavelength columns from imageData_df in place
# imageData_df must be pre-processed with call to adjustColumnTypes()
# pass lambda_start = -1 to omit any modification of imageData_df
# returns the restricted nLambdas, and the truncated dataframe (all cols other than lambda cols are returned as is)
def restrictLambdaRange(imageData_df, lambda_start, lambda_end, firstLambdaCol, nLambdas):
if lambda_start < 0:
return nLambdas, imageData_df
# Get the original lambda columns
lambda_columns = imageData_df.columns[firstLambdaCol:firstLambdaCol+nLambdas]
# Filter lambda columns based on the numeric range
lambda_columns_to_remove = [col for col in lambda_columns if not (lambda_start <= float(col) <= lambda_end)]
# Calculate how many lambda columns were removed
nLambdas_new = nLambdas - len(lambda_columns_to_remove)
# Create the filtered dataframe by dropping columns outside the range
filtered_df = imageData_df.drop(columns=lambda_columns_to_remove)
return nLambdas_new, filtered_df
def combine_arrays_with_varying_columns(coeffs):
num_rows = coeffs[0].shape[0] # Get the number of rows (consistent across arrays)
# Pre-allocate an empty array to store the combined data
combined_data = np.empty((num_rows, 0))
# Combine columns from each array
for arr in coeffs:
combined_data = np.hstack((combined_data, arr))
# Create a DataFrame from the combined data
df = pd.DataFrame(combined_data)
return df
def split_list(L, N): # for splitting wavelet_param_dicts into sublists by CPU
"""
Splits a list L into N sublists as evenly as possible.
Args:
L: The list to split.
N: the number of sublists
Returns:
A list of N sublists.
"""
n = len(L)
# Calculate the base sublist size and the number of extra elements
sublist_size = n // N
extra_elements = n % N
sublists = []
start = 0
for i in range(N):
# Adjust sublist size based on extra elements
end = start + sublist_size + (1 if i < extra_elements else 0)
sublists.append(L[start:end])
start = end
return sublists
# augment imageData by replicating names using_a, _b, etc suffixes
def augXdf(df, X):
# df is input dataframe, typically imageData, X = aug factor
if X <= 1:
uniqueNamesANDclassIDs = df.groupby('name')['classID'].first()
return df, uniqueNamesANDclassIDs
# Shuffle the modified DataFrame to randomize the assignment of subgroups
aug_df = df.sample(frac=1)
# Add a new column for the subgroup labels
aug_df['subgroup'] = aug_df.groupby('name').cumcount() % X
# Create the new name column with the subgroup suffixes
# aug_df['name'] = aug_df.apply(lambda row: f"{row['name']}_{chr(ord('a') + row['subgroup'])}", axis=1) # _a, _b, etc suffixes
aug_df['name'] = aug_df.apply(lambda row: f"{row['name']}_aug{str(row['subgroup']+1)}", axis=1) # _aug1, _aug2, etc suffixes, unlimited
# Drop the temporary 'subgroup' column
aug_df = aug_df.drop(columns='subgroup')
aug_df.sort_values(by=['classID', 'name'], inplace=True)
aug_df.reset_index(drop=True, inplace=True)
# you can check the new names and their classIDs like this: aug_df.groupby('name')['classID'].first()
uniqueNamesANDclassIDs = aug_df.groupby('name')['classID'].first()
return aug_df, uniqueNamesANDclassIDs
def augPow2df(df, target_power):
"""
Augment dataframe so final number of unique names is 2^target_power
Args:
df: Input dataframe with 'name' and 'classID' columns
target_power: Power of 2 desired for final unique name count
"""
# Get current number of unique names
current_unique = len(df['name'].unique())
target_unique = 2**target_power
if current_unique > target_unique:
return df, df.groupby('name')['classID'].first()
# Calculate required multiplication factor
X = int(np.ceil(target_unique / current_unique))
# Shuffle the DataFrame to randomize assignment
aug_df = df.sample(frac=1)
# Add subgroup column
aug_df['subgroup'] = aug_df.groupby('name').cumcount() % X
# Create new names with augmentation suffixes
aug_df['name'] = aug_df.apply(
lambda row: f"{row['name']}_aug{str(row['subgroup']+1)}",
axis=1
)
# Clean up and sort
aug_df = aug_df.drop(columns='subgroup')
aug_df.sort_values(by=['classID', 'name'], inplace=True)
aug_df.reset_index(drop=True, inplace=True)
# If we have more unique names than needed, randomly select names
if len(aug_df['name'].unique()) > target_unique:
keep_names = np.random.choice(
aug_df['name'].unique(),
size=target_unique,
replace=False
)
aug_df = aug_df[aug_df['name'].isin(keep_names)]
uniqueNamesANDclassIDs = aug_df.groupby('name')['classID'].first()
return aug_df, uniqueNamesANDclassIDs
def aug0df(df):
# df is input dataframe, typically imageData, special case for Mk6: each row is augmented to a unique instance
aug_df = df.copy()
aug_df['name'] = aug_df.groupby('name').cumcount().add(1).astype(str).radd(aug_df['name'] + '_')
return aug_df
def balancedTrainTestSplitWithAug(df, train_proportion, augX_train, augX_test=1):
# df must have 'name' and 'classID' cols, together with any number of feature cols
# returns an augmented balanced train-test split with 2 return df's ensuring that rootnames (after augmentation) appear only in one or the other set
# see also balancedTrainTestSplitWithAugFirst
import pandas as pd
from sklearn.model_selection import train_test_split
# Create a new DataFrame with unique names and their classID
name_class_df = df.groupby('name')['classID'].agg(lambda x:x.value_counts().index[0]).reset_index()
# Split unique names into train and test sets while preserving class balance
#train_names, test_names = train_test_split(name_class_df['name'], train_size=train_proportion, random_state=42, stratify=name_class_df['classID'])
train_names, test_names = train_test_split(name_class_df['name'], train_size=train_proportion, stratify=name_class_df['classID'])
# Create boolean masks for train and test sets based on names
train_mask = df['name'].isin(train_names)
test_mask = df['name'].isin(test_names)
# Split the DataFrame into train and test sets using the masks
df_train = df[train_mask]
df_test = df[test_mask]
df_train_aug, uniqueNamesANDclassIDs_train_aug = augXdf(df_train, augX_train) # data augmentation
df_test_aug, uniqueNamesANDclassIDs_test_aug = augXdf(df_test, augX_test) # data augmentation
df_train_aug.reset_index(inplace=True, drop=True) # drop means delete the old index column
df_test_aug.reset_index(inplace=True, drop=True) # drop means delete the old index column
return df_train_aug, df_test_aug, uniqueNamesANDclassIDs_train_aug, uniqueNamesANDclassIDs_test_aug
def balancedTrainTestSplitWithAugFirst(df, augX, train_proportion): # AUGMENTING FIRST IS LIKELY INCORRECT AND WILL RESULT IN DATA LEAKAGE TO TEST SET
# df must have 'name' and 'classID' cols, together with any number of feature cols
# returns an augmented balanced train-test split with 2 return df's NOT ensuring that rootnames (after augmentation) appear only in one or the other set ie. aug is done first and aug'd names are used as-is for the train-test split
# see also balancedTrainTestSplitWithAug
import pandas as pd
from sklearn.model_selection import train_test_split
# aug first:
df_aug,uniqueNamesANDclassIDs_df_aug = augXdf(df, augX)
# Create a new DataFrame with unique names and their most frequent (will all be the same here) classID
name_class_aug_df = df_aug.groupby('name')['classID'].agg(lambda x:x.value_counts().index[0]).reset_index()
# Split unique names into train and test sets while preserving class balance
train_names, test_names = train_test_split(name_class_aug_df['name'], train_size=train_proportion, random_state=42, stratify=name_class_aug_df['classID'])
# Create boolean masks for train and test sets based on names
train_mask = df_aug['name'].isin(train_names)
test_mask = df_aug['name'].isin(test_names)
# Split the DataFrame into train and test sets using the masks
df_train_aug = df_aug[train_mask]
df_test_aug = df_aug[test_mask]
df_train_aug, uniqueNamesANDclassIDs_train_aug = augXdf(df_train_aug, 1) # just to get uniqueNamesANDclassIDs_train_aug
df_test_aug, uniqueNamesANDclassIDs_test_aug = augXdf(df_test_aug, 1)
df_train_aug.reset_index(inplace=True, drop=True) # drop means delete the old index column
df_test_aug.reset_index(inplace=True, drop=True) # drop means delete the old index column
return df_train_aug, df_test_aug, uniqueNamesANDclassIDs_train_aug, uniqueNamesANDclassIDs_test_aug
def averageSpectraByName(spectra_df, lambdaColStartIx, nLambdaCols, sortByClassID=True):
# NAME-WISE AVERAGE OF ALL SPECTRA in spectra_df
# retains only name, classID and lambda cols
lambda_columns = spectra_df.columns[lambdaColStartIx:lambdaColStartIx+nLambdaCols]
lambda_columns = np.insert(lambda_columns, 0, 'classID')
# 2. Group, Aggregate, and Reset Index
df_avg = (
spectra_df.groupby('name')[lambda_columns] # Group by 'name'
.mean() # Calculate mean for each numeric column within groups
.reset_index() # Convert 'name' back into a regular column
)
df_avg['classID'] = df_avg['classID'].astype('Int64') # "mean" classID back to int
if sortByClassID:
df_avg.sort_values(by=['classID', 'name'], inplace=True)
df_avg.reset_index(drop=True, inplace=True)
new_lambdaColStartIx = 2 # because some cols were stripped this is the new lambdaColStartIx (nLambdaCols remains the same)
return df_avg, new_lambdaColStartIx
def get_column_indices(df, startChar='L'):
"""Gets the 0-based indices of columns starting with 'L' in a DataFrame (ITRK can export kernel spectra by prefixing all lambda columns with 'L')
Args:
df: The pandas DataFrame.
Returns:
A list of column indices, # of columns starting with 'L'
"""
l_columns = df.columns[df.columns.str.startswith(startChar)]
indices = [df.columns.get_loc(col) for col in l_columns]
N = len(indices)
if N == 0:
firstIndex = -1
else:
firstIndex = indices[0]
return firstIndex, N, indices
def get_numericalColumns(df):
""" Gets the 0-based indices of numerical columns
Use to detect the lambda columns in an imageDF
Args:
df: The pandas DataFrame.
Returns:
A list of column indices, # of numerical columns
"""
# Function to check if a string is numeric
def is_numeric(column_name):
try:
float(column_name)
return True
except ValueError:
return False
# Find all numeric column names and the index of the first one
numeric_column_names = [col for col in df.columns if is_numeric(col)]
if numeric_column_names:
first_numeric_column_name = numeric_column_names[0]
first_numeric_column_index = df.columns.tolist().index(first_numeric_column_name)
# print(f"The list of numeric column names: {numeric_column_names}")
# print(f"The first numeric column name is: {first_numeric_column_name}")
# print(f"The zero-based index of this column is: {first_numeric_column_index}")
lambdas_np = np.array(numeric_column_names, dtype=float)
return first_numeric_column_index, len(lambdas_np), lambdas_np # firstLamvda, lambdas as strings, convert to np float array like this: )
else:
return -1, 0, np.empty((1), dtype=float)
def custom_growth(num_points, start, end, power):
"""
Generates a monotonic custom growth curve using a power function.
Args:
start: Starting value of the sequence.
end: Ending value of the sequence.
power: Controls the steepness of the curve. Higher values make it steeper.
Must be greater than 0. power=1 results in a linear sequence.
num_points: The number of points to generate in the sequence.
Returns:
Array of output values corresponding to the custom growth curve.
start_value = 0.3
end_value = 10
power_factor = 2 # Higher values make the curve steeper
num_points = 20
y = custom_growth(start_value, end_value, power_factor, num_points)
# Visualization (with markers)
plt.figure(figsize=(10, 6))
plt.plot(y, marker='o', linestyle='-', label=f'Custom Growth (power={power_factor})')
plt.xlabel('Point Index')
plt.ylabel('Y Values')
plt.title('Custom Growth Curve (Power Function)')
plt.grid(True)
plt.legend()
plt.show()
"""
if power <= 0:
raise ValueError("Power must be greater than 0")
x = np.linspace(0, 1, num_points) # Normalized x values
return start + (end - start) * x**power
# Plotting Functions
def plot_spectra(spectra):
"""Plots a set of spectra."""
num_plots = len(spectra)
if num_plots == 1:
fig, ax = plt.subplots(1, num_plots, figsize=(7, 4)) # Changed to ax
ax.plot(spectra[0]) # .spectrum return the 1D spectral data
ax.set_title(f"Spectrum 1")
ax.set_xlabel("Wavelength")
#ax.set_ylabel("Amplitude")
ax.set_yticklabels([]) # This line removes the y-axis tick labels
else: # Multiple plots
fig, ax = plt.subplots(1, num_plots, figsize=(15, 4)) # Changed to ax
for i, spectrum in enumerate(spectra):
ax[i].plot(spectrum) # .spectrum return the 1D spectral data
#x[i].set_title(f"Spectrum {i+1}")
ax[i].set_title(spectrum) # SpectrumData class returns a custom string when referring to a class instance
ax[i].set_xlabel("Wavelength")
#ax[i].set_ylabel("Amplitude")
ax[i].set_yticklabels([]) # This line removes the y-axis tick labels
#plt.show()
# save it to a file:
plt.close()
### MP version of process_spectra
import multiprocessing as mp
from functools import partial
def process_spectra_np_MP_chunk(chunk_with_index, wavelet_name, compression, scales, convertTofloat32): # although it would make sense to put this function inside process_spectra_np_MP this does not work with MP because of a pickling error: keep it separate here
chunk, start_index = chunk_with_index
chunk_scalograms = []
for i, spectrum in enumerate(chunk):
coeffs, _ = pywt.cwt(spectrum, scales, wavelet_name, sampling_period=1.0)
coeffs = CustomWavelet.compress_complex_array(coeffs=coeffs, power_factor=compression)
if convertTofloat32:
if np.iscomplexobj(coeffs):
coeffs = coeffs.astype(np.complex64)
else:
coeffs = coeffs.astype(np.float32)
# Separate real and imaginary parts if complex
# if np.iscomplexobj(coeffs):
# coeffs = np.stack((coeffs.real, coeffs.imag), axis=-1)
# else:
# coeffs = coeffs[..., np.newaxis]
chunk_scalograms.append((start_index + i, coeffs))
return chunk_scalograms
def process_spectra_np_MP(spectra_np, wavelet_name, compression, scales, convertTofloat32=False, maxCPUs=-1, verbose=0): # verbose maybe later
""" Computes CWT scalograms for a set of spectra. Returns a 4D array (N,H,W,C)
May need to limit max CPUs because of memory
"""
num_cores = mp.cpu_count()
if maxCPUs > 0: num_cores = min(maxCPUs,num_cores)
if verbose > 0: printSB(f'process_spectra_np_MP is spawning {num_cores} processes')
total_rows = spectra_np.shape[0]
chunk_size = max(1, total_rows // num_cores)
chunks = [spectra_np[i:i+chunk_size] for i in range(0, total_rows, chunk_size)]
chunks_with_index = [(chunk, i*chunk_size) for i, chunk in enumerate(chunks)] # chunks_with_index is used to ensure that chunks are arranged in the original order of the input array
process_chunk_partial = partial(process_spectra_np_MP_chunk, wavelet_name=wavelet_name, compression=compression, scales=scales, convertTofloat32=convertTofloat32)
with mp.Pool(processes=num_cores) as pool:
results = pool.map(process_chunk_partial, chunks_with_index)
flat_results = [item for sublist in results for item in sublist]
sorted_results = sorted(flat_results, key=lambda x: x[0])
scalograms = [item[1] for item in sorted_results]
return np.array(scalograms)
### END MP version of process_spectra
def process_spectra_np(spectra_np, wavelet_name, compression, scales, convertTofloat32=False):
"""Computes CWT scalograms for a set of spectra. 1 scalogram for each row in spectra_np"""
scalograms = []
for i in range(spectra_np.shape[0]):
spectrum = spectra_np[i, :]
coeffs, freqs = pywt.cwt(spectrum, scales, wavelet_name, sampling_period=1.0)
coeffs = CustomWavelet.compress_complex_array(coeffs=coeffs, power_factor=compression) # compress dynamic range of coeffs
if convertTofloat32: # save memory
if np.iscomplexobj(coeffs): # complex
coeffs = coeffs.astype(np.complex64) # 32-bit single precision real and imag components = 64 bits total
else: # real
coeffs = coeffs.astype(np.float32) # 32-bit single precision real
scalograms.append(coeffs)
return np.array(scalograms)
def spectraToScalograms_2D(X, CWT_scales_list, CWT_type, normalizeWavelets=False, convertTofloat32=True):
# computes a 2D scalogram for each row in X which contains a spectrum. The 2D scalogram is flattened to a 1D vector (and complex scalograms have their real and complex flattened 1D vectors concatenated)
# final result returned as a 2D np array of N (=# of spectra/subjects/instances) rows of 1D vectors that are flattened scalograms
# Iterate through rows of X
scalograms = []
for index, row in enumerate(X):
coef, freqs = pywt.cwt(row, CWT_scales_list, CWT_type, sampling_period=1.0)
coef = CustomWavelet.compress_complex_array(coef) # compressed dynamic range of magnitudes
if convertTofloat32: # save memory
if np.iscomplexobj(coef): # complex
coef = coef.astype(np.complex64) # 32-bit single precision real and imag components
else: # real
coef = coef.astype(np.float32) # 32-bit single precision real
scalograms.append(coef)
scalograms_np = np.array(scalograms) # 3D array (N,CWT_N_scales,nL) of real or complex values
#printSB('scalograms_np.shape1:',scalograms_np.shape)
scalograms_np = scalograms_np.reshape(scalograms_np.shape[0],-1) # reshape to 2D array (N,CWT_N_scales*nL) of real or complex values
#printSB('scalograms_np.shape2:',scalograms_np.shape)
if np.iscomplexobj(scalograms_np): # we must handle complex results by processing mag & phase separately
#printSB(CWT_type + ': scalograms_np is complex')
scalograms_np_mag = np.abs(scalograms_np) # (N,CWT_N_scales*nL) reals
#printSB('scalograms_np_mag.shape:',scalograms_np_mag.shape)
scalograms_np_phase = np.angle(scalograms_np) # (N,CWT_N_scales*nL) imgs
if normalizeWavelets:
row_maxes = np.max(scalograms_np_mag, axis=1) # Find row maximums, handling all-zero rows
row_maxes[row_maxes == 0] = 1 # Replace zeros with 1 to avoid division by zero
scalograms_np_mag = scalograms_np_mag / row_maxes[:, np.newaxis] # Normalize each row
scalograms_np_phase = (scalograms_np_phase + np.pi) / (2 * np.pi) # Rescale phase angles from -pi..pi to to 0..1 to match normalized magnitudes
scalograms_np_combined = np.hstack([scalograms_np_mag, scalograms_np_phase]) # concat mag & phase; what is the shape here? should be (N,CWT_N_scales*nL*2) float32 reals: YES
else: # real only
#printSB(CWT_type + ': scalograms_np is real')
if normalizeWavelets:
row_maxes = np.max(scalograms_np, axis=1) # Find row maximums, handling all-zero rows
row_maxes[row_maxes == 0] = 1 # Replace zeros with 1 to avoid division by zero
scalograms_np = scalograms_np / row_maxes[:, np.newaxis] # Normalize each row
scalograms_np_combined = scalograms_np # pass through (N,CWT_N_scales*nL)
#printSB('Returning from spectraToScalograms_2D: scalograms_np_combined.shape:',scalograms_np_combined.shape)
return scalograms_np_combined # always real
def spectraToScalograms_4D(X_in, CWT_scales_list, CWT_type, compression, normalizeScalograms=False, convertTofloat32=False, verbose=0, MP=True, lambdaStartIx=-1, lambdaEndIx=-1):
# X_in is a 2D np array of features (typically wavelength bins normalized to peak 1.0), each instance in a row
# computes a 2D scalogram for each row in X_in which contains a spectrum. The 2D scalogram is NOT flattened to a 1D vector instead a 2D array/image is returned for each row (subject) in X_in, with 1 (real) or 2 (complex, depending on wavelet type) channels
# normalizeScalograms: normalizes each channel (real +/- imag) of a scalogram to its max val (not norm columns as in Standard Scaler). This is to prevent mag diffs between Re & Im channels
# final result returned as a 4D np array shape (N,H,W,C) where N is instances/names (this dimension will mirror the rows of X_in), H is scalogram height in pixles, W is scalogram width, C is channels (1 real, 2 complex)
# pass verbose>0 to produce a mod verbose log of rows being processed eg verbose=10000 will print a msg every 10000 rows/spectra
# if you ever do eg: avgSpectrum0=np.mean(features0.values, axis=0) to compute a mean spectrum from a bunch of raw spectral rows, then you want a single scalogram from this mean spectrum make you reshape X_in like this: X_in=avgSpectrum0.reshape(1,avgSpectrum0.shape[0])
# rev.2024-09-15: corrected complex normalization and float32 conversion in process_spectra_np
# rev.2025-01-26: now normalizeScalograms=False by default, but make sure you pass peak-normalized spectra in X_in
# rev.2025-03-05: lambdaStartIx, lambdaEndIx INCLUSIVE 0-based indexes to extract a subset of lambdas for computing scalograms, to restrict analysis to a certain wavelength band
if (CWT_type is None) or (CWT_type==''):
# raise ValueError('spectraToScalograms_4D: CWT_type was not defined')
printSB('FATAL ERROR: spectraToScalograms_4D: CWT_type was not defined')
sys.exit()
if (compression is None) or (compression==0):
# raise ValueError('spectraToScalograms_4D: compression was not defined or is 0')
printSB('FATAL ERROR: spectraToScalograms_4D: compression was not defined or is 0')
sys.exit()
if lambdaStartIx > -1: X_in = X_in[:, lambdaStartIx:lambdaEndIx+1] # important to extract the lambda subset BEFORE scalograms calculation?
if MP: # MP version:
if verbose > 0: printSB(f'Starting process_spectra_np_MP to transform {len(X_in)} spectra into scalograms...')
X_cwt = process_spectra_np_MP(X_in, wavelet_name=CWT_type, compression=compression, scales=CWT_scales_list, convertTofloat32=convertTofloat32, verbose=verbose) # X_cwt contains all the scalograms, 1 per mean spectrum, shape (N, height, width)
if verbose > 0: printSB(' DONE')
else: # non-MP version:
if verbose > 0: printSB(f'Starting process_spectra_np to transform {len(X_in)} spectra into scalograms...')
X_cwt = process_spectra_np(X_in, wavelet_name=CWT_type, compression=compression, scales=CWT_scales_list, convertTofloat32=convertTofloat32) # X_cwt contains all the scalograms, 1 per mean spectrum, shape (N, height, width)
if verbose > 0: printSB(' DONE')
# some wavelets return complex data:We separate the real (X_train.real) and imaginary (X_train.imag) components into two separate arrays, then stack them along the last axis to create a 3D with shape (N, height, width) for both real and complex wavelets, where N is instances, H is scales, W is wavelength bins (4D tensor with extra real & im channels is created below)
waveletIsComplex = np.iscomplexobj(X_cwt)
# printSB(f'DEBUG: lambdaStartIx={lambdaStartIx}. lambdaEndIx={lambdaEndIx}')
if waveletIsComplex:
#printSB('spectraToScalograms_4D is processing wavelet: ' + CWT_type + ' (complex)...')
X_cwt_mag = np.abs(X_cwt)
X_cwt_phase = np.angle(X_cwt)
if normalizeScalograms:
max_values = np.max(np.abs(X_cwt_mag), axis=(1, 2)) # Compute max values for each image (along the 2nd and 3rd axes)
max_values_reshaped = max_values[:, np.newaxis, np.newaxis] # Reshape max_values to enable broadcasting
X_cwt_mag = X_cwt_mag / max_values_reshaped # Normalize each image (not column as in StandardScaler) by its max value
# Verify the normalization (optional)
# max_after_normalization = np.max(X_cwt_mag, axis=(1, 2))
# printSB("Max values after normalization:", max_after_normalization) # Should be all 1.0
X_cwt_phase = (X_cwt_phase + np.pi) / (2 * np.pi) # Rescale phase angles from -pi..pi to 0..1 to match normalized magnitudes
X_cwt_stacked = np.stack((X_cwt_mag, X_cwt_phase), axis=-1) # Shape: (N, height, width, 2 channels)
X_out = X_cwt_stacked # this is the new input type
else: # real wavelet coefs
#printSB('spectraToScalograms_4D is processing wavelet: ' + CWT_type + '...')
if normalizeScalograms:
max_values = np.max(np.abs(X_cwt), axis=(1, 2)) # Compute max values for each image (along the 2nd and 3rd axes)
max_values_reshaped = max_values[:, np.newaxis, np.newaxis] # Reshape max_values to enable broadcasting
X_cwt = X_cwt / max_values_reshaped # Normalize each image by its max value
X_out = np.expand_dims(X_cwt, axis=-1) # Reshape to include channel dimension for Conv2D, now shape (N,H,W,C)
# independently standardize the channels (scales of X_mag & X_phase may be very different):
#X_train, train_mean, train_std = standardize_channels(X_train) THIS WAS DONE ABOVE
if verbose > 0: printSB('spectraToScalograms_4D: X_in.shape:',X_in.shape,' X_out.shape:',X_out.shape) # debug
return X_out, waveletIsComplex # shape (N,H,W,C), always real
def spectraToScalograms_3Dcomplex(X_in, CWT_scales_list, CWT_type, compression, convertTofloat32=False, verbose=0, MP=True, lambdaStartIx=-1, lambdaEndIx=-1):
# X_in is a 2D np array of features (typically wavelength bins normalized to peak 1.0), each instance in a row
# computes a 2D scalogram for each row in X_in which contains a spectrum. The 2D scalogram is NOT flattened to a 1D vector instead a 2D array/image is returned for each row (subject) in X_in, with 1 (real) or 2 (complex, depending on wavelet type) channels
# like spectraToScalograms_4D except returns a 3D array shape (N,H,W) of real or complex numbers, depending on the wavelet
if (CWT_type is None) or (CWT_type==''):
raise ValueError('FATAL ERROR: spectraToScalograms_3Dcomplex: CWT_type was not defined')
if (compression is None) or (compression==0):
# raise ValueError('spectraToScalograms_4D: compression was not defined or is 0')
raise ValueError('FATAL ERROR: spectraToScalograms_3Dcomplex: compression was not defined or is 0')
if lambdaStartIx > -1: X_in = X_in[:, lambdaStartIx:lambdaEndIx+1] # important to extract the lambda subset BEFORE scalograms calculation?
if MP: # MP version:
if verbose > 0: printSB(f'Starting process_spectra_np_MP to transform {len(X_in)} spectra into scalograms...')
X_cwt = process_spectra_np_MP(X_in, wavelet_name=CWT_type, compression=compression, scales=CWT_scales_list, convertTofloat32=convertTofloat32, verbose=verbose) # X_cwt contains all the scalograms, 1 per mean spectrum, shape (N, height, width)
if verbose > 0: printSB(' DONE')
else: # non-MP version:
if verbose > 0: printSB(f'Starting process_spectra_np to transform {len(X_in)} spectra into scalograms...')
X_cwt = process_spectra_np(X_in, wavelet_name=CWT_type, compression=compression, scales=CWT_scales_list, convertTofloat32=convertTofloat32) # X_cwt contains all the scalograms, 1 per mean spectrum, shape (N, height, width)
if verbose > 0: printSB(' DONE')
# some wavelets return complex data:We separate the real (X_train.real) and imaginary (X_train.imag) components into two separate arrays, then stack them along the last axis to create a 3D with shape (N, height, width) for both real and complex wavelets, where N is instances, H is scales, W is wavelength bins (4D tensor with extra real & im channels is created below)
waveletIsComplex = np.iscomplexobj(X_cwt)
X_out = X_cwt
if verbose > 0: printSB('spectraToScalograms_3Dcomplex: X_in.shape:',X_in.shape,' X_out.shape:',X_out.shape) # debug
return X_out, waveletIsComplex # shape (N,H,W), real or complex numbers
def plot_scalograms(scalograms):
"""Plots a set of scalograms."""
num_plots = len(scalograms)
if num_plots == 1:
fig, ax = plt.subplots(1, num_plots, figsize=(4, 4))
ax.imshow(np.abs(scalograms[0]), extent=[0, N_points-1, CWT_scales_end, CWT_scales_start], cmap='viridis', aspect='auto') # Adjust extent based on scales
ax.set_title(f"Scalogram 1")
ax.set_ylabel("Scale")
ax.set_xlabel("Wavelength (nm)")
else: # Multiple plots
fig, ax = plt.subplots(1, num_plots, figsize=(15, 4))
for i, scalogram in enumerate(scalograms):
ax[i].imshow(np.abs(scalogram), extent=[0, N_points-1, CWT_scales_end, CWT_scales_start], cmap='viridis', aspect='auto') # Adjust extent based on scales, we use np.abs() in case scalograms are complex
ax[i].set_title(f"Scalogram {i+1}")
ax[i].set_ylabel("Scale")
ax[i].set_xlabel("Wavelength (nm)")
#plt.show()
# save it to a file:
plt.close()
def standardize_channels(X, mean=None, std=None):
"""Standardizes the channels of a 2D array of images.
Args:
X: A NumPy array with shape (num_samples, height, width, channels).
mean: (Optional) Precomputed mean for each channel.
std: (Optional) Precomputed standard deviation for each channel.
Returns:
The standardized array with the same shape as X, along with
the calculated mean and std if not provided.
Example usage:
X_train_standardized, train_mean, train_std = standardize_channels(X_train)
Use the SAME mean and std to standardize X_test
X_test_standardized, _, _ = standardize_channels(X_test, mean=train_mean, std=train_std)
"""
num_channels = X.shape[-1]
if mean is None:
mean = np.mean(X, axis=(0, 1, 2)) # Calculate mean per channel
if std is None:
std = np.std(X, axis=(0, 1, 2)) # Calculate std per channel
for channel in range(num_channels):
X[:, :, :, channel] = (X[:, :, :, channel] - mean[channel]) / std[channel]
return X, mean, std
### plot class0 vs class1 average scalograms (we do this one first else this is the one [most recent] that is displayed in SB)
def plotClassAvgSpectraAndScalograms(imageData, firstLambdaCol, nLambdas, normMode, D, CWT_scales_list, CWT_type, fBaseNoExt, outPDF='', pad=0.5):
# fBaseNoExt: csv input fname, for graphtitle
# fig.tight_layout(pad=pad): pad value may need adjustment
# class-wise avg spectra (or their derivatives)
classCol = imageData['classID']
lambdaCols = imageData.iloc[:,firstLambdaCol:firstLambdaCol+nLambdas]
# lambdas = imageData_avg.columns[firstLambdaCol:].to_numpy().astype(np.float64)
#lambdas = df_total_NONaug_avg.columns[firstLambdaCol_avg:firstLambdaCol_avg+nLambdas].to_numpy().astype(np.float64)
lambdas = imageData.columns[firstLambdaCol:firstLambdaCol+nLambdas].to_numpy().astype(np.float64)
tick_locations = np.arange(400, lambdas.max() + 50, 50) # Start at 400, end at max value + 50, with a step of 50
# Filter lambdaCols based on classCol=0 values
filtered_lambdaCols = lambdaCols[classCol.iloc[:] == 0]
# Calculate the mean across rows
class0avg = filtered_lambdaCols.mean(axis=0)
# Filter lambdaCols based on classCol=1 values
filtered_lambdaCols = lambdaCols[classCol.iloc[:] == 1]
# Calculate the mean across rows
class1avg = filtered_lambdaCols.mean(axis=0)
if normMode == 1:
class0avg /= class0avg.max() # normalize to 1.0
class1avg /= class1avg.max()
# derivatives?
if D==1: # 1st derivative
class0avg = np.pad(np.diff(class0avg), pad_width=(0,1), mode='edge') # pad to retain same num of elements else 1 fewer than original
class1avg = np.pad(np.diff(class1avg), pad_width=(0,1), mode='edge')
elif D==2: # 2nd derivative
class0avg = np.pad(np.diff(class0avg), pad_width=(0,1), mode='edge') # pad to retain same num of elements else 1 fewer than original
class0avg = np.pad(np.diff(class0avg), pad_width=(0,1), mode='edge')
class1avg = np.pad(np.diff(class1avg), pad_width=(0,1), mode='edge')
class1avg = np.pad(np.diff(class1avg), pad_width=(0,1), mode='edge')
### 3-panel graph
# plot formats:
nSubplotRows = 1
nSubplotCols = 3
figHeight_inches = 5
figWidth_inches = 2.07*figHeight_inches # match aspect ratio of canvas for nice display
plotTitleFontSize = 12
axisLabelFontSize = 13
tickLabelFontSize = 11
statsStrFontSize = 10
nameLabelsFontSize = 7
fig, axs = plt.subplots(nSubplotRows,nSubplotCols,figsize=(figWidth_inches, figHeight_inches))
fig.tight_layout(pad=pad) # adjust spacing between subplots; https://www.geeksforgeeks.org/how-to-set-the-spacing-between-subplots-in-matplotlib-in-python/
### panel 1: class-wise avg spectra (or their derivatives)
ax=plt.subplot(nSubplotRows, nSubplotCols, 1) # left panel
graphTitle = fBaseNoExt + '\nAveraged Spectra'
if D==0:
ylabel = 'Intensity'
elif D==1: # 1st derivative
graphTitle = graphTitle + ' (D1)'
ylabel = 'D1'
elif D==2: # 2nd derivative
graphTitle = graphTitle + ' (D2)'
ylabel = 'D2'
ax.set_title(graphTitle,size=plotTitleFontSize-1)
ax.set_xlabel('Wavelength (nm)', size= axisLabelFontSize)
if kNormalizeSpectralOverlay:
ax.set_ylabel('Normalized ' + ylabel, size= axisLabelFontSize)
else:
ax.set_ylabel(ylabel, size= axisLabelFontSize)
ax.tick_params(axis='both', labelsize=tickLabelFontSize, direction='out', length=10, bottom=True, left=True, width=2) # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.tick_params.html
#ax.set_xticks([0,1],['0 (n=' + str(score_df[score_df['classID'] == 0].shape[0]) + ')','1 (n=' + str(score_df[score_df['classID'] == 1].shape[0]) + ')']) # classID labels and n's
# Define the tick locations
# Get the numerical columns for plotting
ax.grid(axis='y', color='0.8', linewidth=1.0)
#ax.axhline(0, color='blue')
#ax.legend(['Class 0', 'Class 1'])
#lambdas = imageData_avg.columns[firstLambdaCol:].to_numpy().astype(np.float64)
#tick_locations = np.arange(400, lambdas.max() + 50, 50) # Start at 400, end at max value + 50, with a step of 50
# Set the tick locator
ax.set_xlim([tick_locations.min(), tick_locations.max()]) # Set limits to encompass ticks
ax.plot(lambdas, class0avg, color='g',linewidth=0.5)
ax.plot(lambdas, class1avg, color='r',linewidth=0.5)
# Add the legend AFTER plotting
ax.legend(['class0', 'class1'], loc='upper right', fontsize=11)
### panel 2: mean scalogram
mean_scalogram, freqs = pywt.cwt((class0avg + class1avg) / 2, CWT_scales_list, CWT_type, sampling_period=1.0) # mean scalogram from both classes
mean_scalogram = CustomWavelet.compress_complex_array(mean_scalogram) # compressed dynamic range of magnitudes
ax=plt.subplot(nSubplotRows, nSubplotCols, 2) # middle panel
graphTitle = 'Mean scalogram (' + CWT_type
if D==0:
graphTitle = graphTitle + ')'
elif D==1: # 1st derivative
graphTitle = graphTitle + ', D1)'
elif D==2: # 2nd derivative
graphTitle = graphTitle + ', D2)'
ax.set_title(graphTitle,size=plotTitleFontSize)
ax.tick_params(axis='both', labelsize=tickLabelFontSize, direction='out', length=10, bottom=True, left=True, width=2) # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.tick_params.html
# ax.imshow(np.abs(class0scalogram), extent=[lambdas[0], lambdas[-1], CWT_scales_list[-1], CWT_scales_list[0]], cmap='viridis', aspect='auto') # Adjust extent based on scales
ax.imshow(np.abs(mean_scalogram), extent=[lambdas[0], lambdas[-1], len(CWT_scales_list)-1, 0], cmap='viridis', aspect='auto') # extent is set linearly to the # of scales: we re-label the ticks with the real non-linear values below
ax.set_ylabel("Scale (higher values=lower frequencies)", size= axisLabelFontSize-1)
ax.set_xlabel("Wavelength (nm)", size= axisLabelFontSize)
num_ticks = 10
y_min, y_max = 0, len(CWT_scales_list)-1
y_ticks_positions = np.linspace(y_min, y_max, num_ticks, dtype=int)
y_ticks_values = CWT_scales_list[y_ticks_positions]
# Format tick labels based on value
y_tick_labels = [f"{val:.1f}" if val < 10 else f"{val:.0f}" for val in y_ticks_values]
# Set y-ticks and labels with fontsize
ax.set_yticks(y_ticks_positions) # this sets the even tick spacing along the y-axis
ax.set_yticklabels(y_tick_labels,fontsize=tickLabelFontSize-1) # here we apply the corrrect numerical labels for the non-linear CWT_scales_list
# difference scalogram
class0scalogram, freqs = pywt.cwt(class0avg, CWT_scales_list, CWT_type, sampling_period=1.0)
class1scalogram, freqs = pywt.cwt(class1avg, CWT_scales_list, CWT_type, sampling_period=1.0)
class0scalogram = CustomWavelet.compress_complex_array(class0scalogram) # compressed dynamic range of magnitudes
class1scalogram = CustomWavelet.compress_complex_array(class1scalogram) # compressed dynamic range of magnitudes
# Calculate difference image
diff_image = np.abs(class1scalogram) - np.abs(class0scalogram)
mean_absolute_difference_per_pixel = np.mean(np.abs(diff_image)) / diff_image.size
#printSB('diff_image mean_absolute_difference_per_pixel:',mean_absolute_difference_per_pixel) # later could be use this to find the wavelet that generates the largest difference between class0 & 1?
# Create custom colormap (blue for negative, red for positive)
from matplotlib.colors import LinearSegmentedColormap
from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar
import matplotlib.font_manager as fm
from matplotlib.transforms import Bbox
import matplotlib as mpl
colors = [(0, 0, 1), (1, 0, 0)] # Blue to Red
cm = LinearSegmentedColormap.from_list("BlueRed", colors, N=256)
cm = mpl.colormaps['coolwarm'] # canned LUTs, eg coolwarm, seismic, RdBu_r, bwr
ax=plt.subplot(nSubplotRows, nSubplotCols, 3) # right panel
graphTitle = 'Mean difference scalogram'
ax.set_title(graphTitle,size=plotTitleFontSize)
ax.tick_params(axis='both', labelsize=tickLabelFontSize, direction='out', length=10, bottom=True, left=True, width=2) # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.tick_params.html
im = ax.imshow(diff_image, extent=[lambdas[0], lambdas[-1], len(CWT_scales_list)-1, 0], cmap =cm, aspect='auto') # Adjust extent based on scales
ax.set_ylabel("Scale (higher values=lower frequencies)", size= axisLabelFontSize-1)
ax.set_xlabel("Wavelength (nm)", size= axisLabelFontSize)
# Add scale bar (adjusted for subplots)
# Get the position of the last subplot
pos = ax.get_position()
#cbar = plt.colorbar(im, ax=ax, label='Diff')
#cbar = plt.colorbar(im, ax=ax)
# Define the colorbar position and width (adjust 'width' as needed)
cax = fig.add_axes([pos.x1 + 0.01, pos.y0, 0.02, pos.height]) # The pos.x1 + 0.02 in fig.add_axes() places the colorbar slightly to the right of the last subplot. The 0.02 controls the width of the colorbar (you can adjust this to your preference). The pos.height ensures the colorbar has the same height as the last subplot.
cbar = plt.colorbar(im, cax=cax)
cbar.ax.tick_params(labelsize=8, width=0.75, pad=1) # Colorbar tick labels and width, pad=2 moves the labels closer to the bar
cbar.ax.set_ylabel(cbar.ax.get_ylabel(), fontsize=8) # Colorbar label (title)
cbar.outline.set_linewidth(0.75) # Colorbar frame linewidth
# Shorten colorbar ticks
for tick in cbar.ax.get_yticklines():
tick.set_markersize(3) # Adjust the markersize (length) of the ticks
# Set y-ticks and labels
ax.set_yticks(y_ticks_positions) # this sets the even tick spacing along the y-axis
ax.set_yticklabels(y_tick_labels,fontsize=tickLabelFontSize-1) # here we apply the corrrect numerical labels for the non-linear CWT_scales_list
if outPDF == '':
plt.show() # presumably running interactively in jupyterlab
else:
plt.savefig(outPDF, bbox_inches='tight')
printSB('Avg scalograms saved to: ' + outPDF)
plt.close() # clear for next plot else we have old annotations persisting
return lambdas
### save average scalograms to avg_scalograms_outPDF (we do this one first else this is the one [most recent] that is displayed in SB)
def saveMeanAndDiffScalograms(class_col_np, featureCols_avg_np, lambdas, normalize, D, CWT_scales_list, CWT_type, outPDF):
# featureCols_avg_np: 2D np array of class-wise mean spectra (±augmentation, ±differentiation)
# class_col_np: 1D np vector of classIDs (can convert a df to an np array like this: class_col_df.to_numpy(copy=False).astype(int) )
# lambdas: np array of wavelengths for featureCols_avg_np matrix
# normalize: boolean, normalize mean spectra to 1.0 peak
# D: 0,1,2, derivative, just for adjusting plot labels, differentiation of featureCols_avg_np must be done prior to the call
kAvgSpectraLinewidth = 1.0
tick_min = (lambdas[0] // 50) * 50 # next lower multiple of 50
tick_max = math.ceil(lambdas[-1] / 50) * 50
tick_locations = np.arange(tick_min, tick_max+50, 50) # steps of 50
# Filter rows based on class0
mask = (class_col_np == 0)
filtered_rows = featureCols_avg_np[mask]
# Calculate the mean
class0avg = np.mean(filtered_rows, axis=0)
# Filter rows based on class1
mask = (class_col_np == 1)
filtered_rows = featureCols_avg_np[mask]
class1avg = np.mean(filtered_rows, axis=0)
if normalize:
class0avg /= class0avg.max() # normalize to 1.0
class1avg /= class1avg.max()
### 3-panel graph: mean subject-wise spectra, mean scalogram, class-diff scalogram
# plot formats:
nSubplotRows = 1
nSubplotCols = 3
figHeight_inches = 5
figWidth_inches = 2.07*figHeight_inches # match aspect ratio of canvas for nice display
plotTitleFontSize = 12
axisLabelFontSize = 14
tickLabelFontSize = 12
statsStrFontSize = 10
nameLabelsFontSize = 7
fig, axs = plt.subplots(nSubplotRows,nSubplotCols,figsize=(figWidth_inches, figHeight_inches))
fig.tight_layout(pad=1.0) # adjust spacing between subplots; https://www.geeksforgeeks.org/how-to-set-the-spacing-between-subplots-in-matplotlib-in-python/
### panel1: class-wise avg spectra (or their derivatives)
graphTitle = 'Averaged Spectra'
if D==0:
ylabel = 'Intensity'
elif D==1: # 1st derivative
graphTitle = graphTitle + ' (D1)'
ylabel = 'D1'
elif D==2: # 2nd derivative
graphTitle = graphTitle + ' (D2)'
ylabel = 'D2'
ax=plt.subplot(nSubplotRows, nSubplotCols, 1) # left panel
ax.set_title(graphTitle,size=plotTitleFontSize)
ax.set_xlabel('Wavelength (nm)', size= axisLabelFontSize)
if normalize:
ax.set_ylabel('Normalized ' + ylabel, size=axisLabelFontSize)
else:
ax.set_ylabel(ylabel, size=axisLabelFontSize)
ax.tick_params(axis='both', labelsize=tickLabelFontSize, direction='out', length=10, bottom=True, left=True, width=2) # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.tick_params.html
#ax.set_xticks([0,1],['0 (n=' + str(score_df[score_df['classID'] == 0].shape[0]) + ')','1 (n=' + str(score_df[score_df['classID'] == 1].shape[0]) + ')']) # classID labels and n's
# Define the tick locations
# Get the numerical columns for plotting
ax.grid(axis='y', color='0.8', linewidth=1.0)
#ax.axhline(0, color='blue')
#ax.legend(['Class 0', 'Class 1'])
#columns_to_plot = imageData_avg.columns[firstLambdaCol:].to_numpy().astype(np.float64)
#tick_locations = np.arange(400, columns_to_plot.max() + 50, 50) # Start at 400, end at max value + 50, with a step of 50
# Set the tick locator
ax.set_xlim([tick_locations.min(), tick_locations.max()]) # Set limits to encompass ticks
ax.plot(lambdas, class0avg, color='g',linewidth=kAvgSpectraLinewidth)
ax.plot(lambdas, class1avg, color='r',linewidth=kAvgSpectraLinewidth)
### panel 2: mean scalogram
mean_scalogram, freqs = pywt.cwt((class0avg + class1avg) / 2, CWT_scales_list, CWT_type, sampling_period=1.0) # mean scalogram from both classes
mean_scalogram = CustomWavelet.compress_complex_array(mean_scalogram) # compressed dynamic range of magnitudes
ax=plt.subplot(nSubplotRows, nSubplotCols, 2) # middle panel
graphTitle = 'Mean scalogram (' + CWT_type
if D==0:
graphTitle = graphTitle + ')'
elif D==1: # 1st derivative
graphTitle = graphTitle + ', D1)'
elif D==2: # 2nd derivative
graphTitle = graphTitle + ', D2)'
ax.set_title(graphTitle,size=plotTitleFontSize)
ax.tick_params(axis='both', labelsize=tickLabelFontSize, direction='out', length=10, bottom=True, left=True, width=2) # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.tick_params.html
# ax.imshow(np.abs(class0scalogram), extent=[lambdas[0], lambdas[-1], CWT_scales_list[-1], CWT_scales_list[0]], cmap='viridis', aspect='auto') # Adjust extent based on scales
ax.imshow(np.abs(mean_scalogram), extent=[lambdas[0], lambdas[-1], len(CWT_scales_list)-1, 0], cmap='viridis', aspect='auto') # extent is set linearly to the # of scales: we re-label the ticks with the real non-linear values below
ax.set_ylabel("Scale (higher values=lower frequencies)", size= axisLabelFontSize-1)
ax.set_xlabel("Wavelength (nm)", size= axisLabelFontSize)
num_ticks = 10
y_min, y_max = 0, len(CWT_scales_list)-1
y_ticks_positions = np.linspace(y_min, y_max, num_ticks, dtype=int)
y_ticks_values = CWT_scales_list[y_ticks_positions]
# Format tick labels based on value
y_tick_labels = [f"{val:.1f}" if val < 10 else f"{val:.0f}" for val in y_ticks_values]
# Set y-ticks and labels with fontsize
ax.set_yticks(y_ticks_positions) # this sets the even tick spacing along the y-axis
ax.set_yticklabels(y_tick_labels,fontsize=tickLabelFontSize-1) # here we apply the corrrect numerical labels for the non-linear CWT_scales_list
# difference scalogram
class0scalogram, freqs = pywt.cwt(class0avg, CWT_scales_list, CWT_type, sampling_period=1.0)
class1scalogram, freqs = pywt.cwt(class1avg, CWT_scales_list, CWT_type, sampling_period=1.0)
class0scalogram = CustomWavelet.compress_complex_array(class0scalogram) # compressed dynamic range of magnitudes
class1scalogram = CustomWavelet.compress_complex_array(class1scalogram) # compressed dynamic range of magnitudes
# Calculate difference image
diff_image = np.abs(class1scalogram) - np.abs(class0scalogram)
#mean_absolute_difference_per_pixel = np.mean(np.abs(diff_image)) / diff_image.size
#printSB('diff_image mean_absolute_difference_per_pixel:',mean_absolute_difference_per_pixel) # later could be use this to find the wavelet that generates the largest difference between class0 & 1?
# Create custom colormap (blue for negative, red for positive)
from matplotlib.colors import LinearSegmentedColormap
from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar
import matplotlib.font_manager as fm
from matplotlib.transforms import Bbox
import matplotlib as mpl
colors = [(0, 0, 1), (1, 0, 0)] # Blue to Red
cm = LinearSegmentedColormap.from_list("BlueRed", colors, N=256)
cm = mpl.colormaps['coolwarm'] # canned LUTs, eg coolwarm, seismic, RdBu_r, bwr
ax=plt.subplot(nSubplotRows, nSubplotCols, 3) # right panel
graphTitle = 'Difference scalogram'
ax.set_title(graphTitle,size=plotTitleFontSize)
ax.tick_params(axis='both', labelsize=tickLabelFontSize, direction='out', length=10, bottom=True, left=True, width=2) # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.tick_params.html
im = ax.imshow(diff_image, extent=[lambdas[0], lambdas[-1], len(CWT_scales_list)-1, 0], cmap =cm, aspect='auto') # Adjust extent based on scales
ax.set_ylabel("Scale (higher values=lower frequencies)", size= axisLabelFontSize-1)
ax.set_xlabel("Wavelength (nm)", size= axisLabelFontSize)
# Add colorbar bar (adjusted for subplots)
# Get the position of the last subplot
pos = ax.get_position()
#cbar = plt.colorbar(im, ax=ax, label='Diff')
#cbar = plt.colorbar(im, ax=ax)
# Define the colorbar position and width (adjust 'width' as needed)
cax = fig.add_axes([pos.x1 + 0.02, pos.y0, 0.02, pos.height]) # The pos.x1 + 0.02 in fig.add_axes() places the colorbar slightly to the right of the last subplot. The 0.02 controls the width of the colorbar (you can adjust this to your preference). The pos.height ensures the colorbar has the same height as the last subplot.
cbar = plt.colorbar(im, cax=cax)
cbar.ax.tick_params(labelsize=8, width=0.75, pad=-5) # Colorbar tick labels and width, pad=2 moves the labels closer to the bar
cbar.ax.set_ylabel(cbar.ax.get_ylabel(), fontsize=8) # Colorbar label (title)
cbar.outline.set_linewidth(0.75) # Colorbar frame linewidth
# Shorten colorbar ticks
for tick in cbar.ax.get_yticklines():
tick.set_markersize(3) # Adjust the markersize (length) of the ticks
# Set y-ticks and labels
ax.set_yticks(y_ticks_positions) # this sets the even tick spacing along the y-axis
ax.set_yticklabels(y_tick_labels,fontsize=tickLabelFontSize-1) # here we apply the corrrect numerical labels for the non-linear CWT_scales_list
if outPDF == '':
plt.show() # presumably running interactively in jupyterlab
else:
plt.savefig(outPDF, bbox_inches='tight')
printSB('Avg scalograms saved to: ' + outPDF)
plt.close() # clear for next plot else we have old annotations persisting
return
def generate2classNoisyGaussians(N, num_points, x0, x1, width, noise_level, kRand, amplitude=1.0, start=400, end=750, includePlot=True):
"""
Generate 2 sets of noisy gaussians (total N instances) centered @ x0 (class0) and x1 (class1)
Set kRand to True to randomize labels
Example usage:
N = int(100) # Number of samples
nL = 32 # Original number of features
center0=550
center1=560
width=50
noise_level=0.1
start=400
end=750
amplitude=1.0
kRand = False # Set to True to randomize labels
lambdas, X, y, X0, X1, y0, y1 = generate2classNoisyGaussians(N, nL, center0, center1, width, noise_level, kRand, amplitude, start, end)
plt.figure(figsize=(10, 6))
plt.plot(lambdas, X0[0], label='class0 Gaussian')
plt.plot(lambdas, X1[0], label='class1 Gaussian')
plt.xlabel('lambdas')
plt.ylabel('Intensity')
plt.title('2-class Noisy Gaussian Data')
plt.legend()
plt.grid(True)
plt.show()
"""
import matplotlib.pyplot as plt
import numpy as np
def generate_noisy_gaussians(N, num_points, center, width, noise_level, amplitude=1.0, start=400, end=750):
"""
Generates N 1D noisy Gaussian arrays and returns them as a 2D numpy array together with 0 or 1 labels
Args:
start (int): The starting value of the x-range.
end (int): The ending value of the x-range.
center (int): The center of the Gaussian distribution.
width (int): The width (standard deviation) of the Gaussian distribution.
amplitude (float): The peak amplitude of the Gaussian distribution.
noise_level (float): The percentage of noise to add (between 0 and 1).
num_points (int): The number of points in each generated array.
N (int): The number of noisy Gaussian arrays to generate.
Returns:
lambdas
np.ndarray: A 2D numpy array containing N noisy Gaussian arrays.
class labels
"""
x = np.linspace(start, end, num_points)
def gaussian(x, mu, sigma, amplitude):
return amplitude * np.exp(-((x - mu) / sigma) ** 2 / 2)
y = gaussian(x, center, width, amplitude)
noisy_gaussians = []
for _ in range(N):
noise = np.random.normal(0, noise_level, num_points)
y_noisy = np.clip(y + noise, 0, None) # avoid neg values, 'None' means no upper bound for clipping
noisy_gaussians.append(y_noisy)
return x, np.array(noisy_gaussians)
# Set a seed for reproducibility
#np.random.seed(42)
# Generate X with two sets of Gaussians
lambdas, X0 = generate_noisy_gaussians(N // 2, num_points, x0, width, noise_level, amplitude, start, end)
_, X1 = generate_noisy_gaussians(N // 2, num_points, x1, width, noise_level, amplitude, start, end)
X = np.concatenate((X0, X1), axis=0)
# generate clean gaussians for the plot overlay
num_points_clean = 200
lambdas_clean, X0_clean = generate_noisy_gaussians(1, num_points_clean, x0, width, 0, amplitude, start, end)
_, X1_clean = generate_noisy_gaussians(1, num_points_clean, x1, width, 0, amplitude, start, end)
# Generate initial y ie the 2 classes
y0 = np.zeros(N // 2, dtype=int)
y1 = np.ones(N // 2, dtype=int)
y_initial = np.concatenate((y0, y1), axis=0)
# Randomize labels if kRand is True
if kRand:
y = shuffle(y_initial)
else:
y = y_initial
if includePlot:
plt.figure(figsize=(13, 6))
# plot 1st noisy gaussians from each class
plt.plot(lambdas, X0[0], label='class0', color='g')
plt.plot(lambdas, X1[0], label='class1', color='r')
# clean base gaussians
plt.plot(lambdas_clean, X0_clean[0], label='class0 base', linestyle='--', dashes=(6, 7), color='g', linewidth=0.7)
plt.plot(lambdas_clean, X1_clean[0], label='class1 base', linestyle='--', dashes=(6, 7), color='r', linewidth=0.7)
plt.xlabel('Lambda')
plt.ylabel('Intensity')
plt.title('Two-class Noisy Gaussians (1st from each class plotted)')
plt.legend()
plt.grid(True)
#plt.show()
else:
plt = None
return lambdas, X, y, X0, X1, y0, y1, plt
def generate2classGaussianImageData(nRowsPerSubject, N, nL, x0, x1, width, noise_level, outputCSV='', kRand=False, amplitude=1.0, start=400, end=750):
"""
Generate a 2-class multirow synthetic gaussian imageData df for testing
nRowsPerSubject is # rows/subject
If outputCSV != '' writes the df to a csv as well
Example usage:
nRowsPerSubject = int(1e2)
N = 5+5 # total number of unique names split among 2 classes
nL = 32
x0 = 550 # center wavelength for class0
x1 = 560 # center wavelength for class1
width = 50
noise_level = 0.05
kRand = False
outputCSV = '/Users/pstys/Documents/imageDataGaussians_df.csv'
imageDataGaussians_df, plt = generate2classGaussianImageData(nRows, nSubjects, nL, x0, x1, width, noise_level, outputCSV, kRand)
plt.show
"""
import pandas as pd
import numpy as np
def generate_name_array(N, nSubjects, prefix): # prefix like 'X0_'
"""
Generates a 1D array of N strings with format 'X0_xxx'
Args:
N: Total number of strings to generate.
nSubjects: Number of groups to divide the strings into.
Returns:
A 1D numpy array of strings.
"""
strings = np.empty(shape=N, dtype=object)
group_size = N // nSubjects
for i in range(N):
group_index = i // group_size
strings[i] = prefix + f"{group_index:03d}" # Pad group index with zeros
return strings
lambdas, X, y, X0, X1, y0, y1, plt = generate2classNoisyGaussians(nRowsPerSubject * N, nL, x0, x1, width, noise_level, kRand, amplitude, start, end)
# dummy cols with simulated Orange headers for import
dummy_data = np.full(nRowsPerSubject * N // 2, 100)
intensity_df = pd.DataFrame({'C#intensity': dummy_data})
X_df = pd.DataFrame({'C#X': dummy_data})
Y_df = pd.DataFrame({'C#Y': dummy_data})
Z_df = pd.DataFrame({'C#Z': dummy_data})
unique4_df = pd.DataFrame({'S#unique4': dummy_data})
# lambda col names
lambdas_col_names = np.array(["L" + format(x, "003.1f") for x in lambdas]) # like 'L411.3' etc
# generate the class0 sub-df
class0names = generate_name_array(nRowsPerSubject * N // 2, nSubjects // 2, 'X0_') # like 'X0_000', 'X0_001', etc
# Create a dictionary to map column names to NumPy arrays
class0data = {'mS#name': class0names, 'cD#classID': y0}
class0df = pd.DataFrame(class0data)
spectra0df = pd.DataFrame(X0, columns=lambdas_col_names)
# complete class0 table
class0df = pd.concat([class0df,intensity_df,spectra0df,X_df,Y_df,Z_df,unique4_df],axis=1)
# generate the class1 sub-df
class1names = generate_name_array(nRowsPerSubject * N // 2, nSubjects // 2, 'X1_') # like 'X1_000', 'X1_001', etc
# Create a dictionary to map column names to NumPy arrays
class1data = {'mS#name': class1names, 'cD#classID': y1}
class1df = pd.DataFrame(class1data)
spectra1df = pd.DataFrame(X1, columns=lambdas_col_names)
# complete class0 table
class1df = pd.concat([class1df,intensity_df,spectra1df,X_df,Y_df,Z_df,unique4_df],axis=1)
imageDataGaussians_df = pd.concat([class0df,class1df],axis=0)
if outputCSV != '': imageDataGaussians_df.to_csv(outputCSV,index = False, mode='w') # append df data to line1, if any, omitting the index column
return imageDataGaussians_df, plt
def generate2classJitteredGaussianImageData(n_root_subjects_per_class = 5,
n_aug = 100, # number of "augmented" spectra (ie averaged subgroup of spectra) per root subject; all n_aug spectra will have the same center wavelength
nL = 32, # number of points per gaussin (wavelength bins per spectrum)
noise_level = 0.05,
width = 30, # FWHM
gauss_ampl = 1.0, # ampl of gaussians
center_wavelength_class0 = 550, # all class0 gaussians +/-jitter
center_wavelength_class1 = 555, # all class5 gaussians +/-jitter
center_wavelength_jitter_per_class = 5, # center wavelengths of n_root_subjects_per_class in each class will be jittered by this amount; all n_aug spectra per rootname will have the same jittered centered wavelength, but will differ only by noise
outputCSV='',
includePlot=True
):
gauss_list=[]
name_list=[]
classID_list=[]
for n_root_subjects_per_class_ctr in range(n_root_subjects_per_class):
jitter0 = random.uniform(-center_wavelength_jitter_per_class/2 , center_wavelength_jitter_per_class/2)
jitter1 = random.uniform(-center_wavelength_jitter_per_class/2 , center_wavelength_jitter_per_class/2)
for n_aug_ctr in range(n_aug):
# generate a jittered pair of class0 and class1 gaussians. The class jitter will be the same for all n_aug spectra for this root/class, but they will differ by noise only
lambdas, _, _, X0, X1, y0, y1, _ = generate2classNoisyGaussians(N=2,
num_points=nL,
x0=center_wavelength_class0+jitter0,
x1=center_wavelength_class1+jitter1,
width=width,
noise_level=noise_level,
kRand=False,
amplitude=1.0,
start=400,
end=750,
includePlot=False)
name_list.append(f'class0_{n_root_subjects_per_class_ctr+1}')
classID_list.append(y0[0])
gauss_list.append(X0[0])
name_list.append(f'class1_{n_root_subjects_per_class_ctr+1}')
classID_list.append(y1[0])
gauss_list.append(X1[0])
# combine into a df
# Create a list of dictionaries
lambdas = round_to_N_sig(lambdas, N=3)
data = []
for name, classID, gauss in zip(name_list, classID_list, gauss_list):
row = {
'mS#name': name,
'cD#classID': classID,
'C#intensity': 100, # dummy
}
# Add the gauss array values with lambda values as column names
for i, lambda_val in enumerate(lambdas):
row[f'L{lambda_val:g}'] = gauss[i]
row['C#X']=100 # dummy
row['C#Y']=100 # dummy
row['C#Z']=100 # dummy
row['S#unique4']=100 # dummy
data.append(row)
# Create the DataFrame
df = pd.DataFrame(data)
if outputCSV != '': df.to_csv(outputCSV,index = False, mode='w') # append df data to line1, if any, omitting the index column
if includePlot:
import itertools
# Get the lambda columns (excluding 'name' and 'classID')
#lambda_columns = [col.lstrip('L') for col in df.columns if col not in ['mS#name', 'cD#classID']]
lambda_columns = df.columns[3:3+len(lambdas)]
# Create a color map
color_map = {0: 'green', 1: 'red'}
# Create the plot
plt.figure(figsize=(12, 6))
max_rows = 2 # Specify the maximum number of rows you want to iterate over
for index, row in itertools.islice(df.iterrows(), max_rows):
color = color_map[row['cD#classID']]
# Convert lambda values to float and plot
x = np.array([float(col.lstrip('L')) for col in lambda_columns])
y = row[lambda_columns].values
plt.plot(x, y, color=color, alpha=0.5)
# Customize the plot
plt.xlabel('Lambda')
plt.ylabel('Feature Value')
plt.title('Sample Class0 vs Class1 Gaussians')
# Add a legend
plt.plot([], [], color='green', label='Class 0')
plt.plot([], [], color='red', label='Class 1')
plt.legend()
# Show the plot
plt.grid(True)
#plt.show()
return df, plt
else:
return df, None
import psutil
def get_ram_infoinGB():
mem = psutil.virtual_memory()
total_ram = mem.total / (1024**3)
available_ram = mem.available / (1024**3)
used_ram = mem.used / (1024**3)
percent_used = mem.percent
return total_ram, available_ram, used_ram, percent_used
def adaptive_multi_split(X_df, y_df, N, D=2, test_size=0.2):
# X_df has a 'name' and many features cols
# y_df has a 0/1 target labels col
# returns a list of train-test splits so that taken together each name appears at most N times in all test sets, and no less than N-D times
from collections import Counter
def check_difference(report_df, D):
"""
Check if the difference between the smallest and greatest value in 'test_set_count' is <= D.
Parameters:
- report_df: DataFrame with 'test_set_count' column
- D: Difference threshold
Returns:
- True if difference between smallest and greatest value is <= D, else False
"""
min_count = report_df['test_set_count'].min()
max_count = report_df['test_set_count'].max()
return (max_count - min_count) <= D
def multi_split_and_select_subset(X_df, y_df, N, test_size=0.2, random_state=None):
"""
Performs multiple train_test_splits and selects a subset of test sets to ensure each 'name'
appears exactly N times in the combined test sets.
Args:
X_df: Features DataFrame with 'name' column.
y_df: Target variable DataFrame.
N: Exact number of times each 'name' should appear in the final combined test sets.
test_size: Proportion of data to include in each test split.
random_state: Seed for reproducibility.
Returns:
A list of tuples, where each tuple contains the train and test sets from the selected splits.
A DataFrame containing 'name' and the count of their appearances in the final combined test sets (which should all be N).
"""
name_counts = X_df['name'].value_counts().to_dict()
name_test_counts = {name: 0 for name in name_counts}
all_splits = []
# Perform splits until we have enough data to select the desired subset
while any(count < N for count in name_test_counts.values()):
X_train, X_test, y_train, y_test = train_test_split(
X_df, y_df, test_size=test_size, random_state=random_state, stratify=y_df
)
all_splits.append((X_train, X_test, y_train, y_test))
# Update counts for names in this test set
for name in X_test['name']:
name_test_counts[name] += 1
# Select a subset of splits to achieve the exact N count for each name
selected_splits = []
final_name_counts = Counter()
for X_train, X_test, y_train, y_test in all_splits:
temp_counts = final_name_counts + Counter(X_test['name'])
if all(count <= N for count in temp_counts.values()):
selected_splits.append((X_train, X_test, y_train, y_test))
final_name_counts = temp_counts
report_df = pd.DataFrame(list(final_name_counts.items()), columns=['name', 'test_set_count'])
return selected_splits, report_df
nIters=0
while True:
# Code to be executed repeatedly
nIters += 1
selected_splits, report_df = multi_split_and_select_subset(X_df, y_df, N, test_size)
if check_difference(report_df, D): # Check the "until" condition
break # Exit the loop if the condition is met
printSB('nIters:',nIters)
return selected_splits, report_df
def RAND(kRANDiters, dimReductionMode, bestDict):
########################### RANDOMIZE USING BEST HYPERPARAMS ###########################
# fetch the best estimator from the grid search
grid_search = bestDict['grid_search']
best_pipe = bestDict['best_pipe']
cv = bestDict['cv'] # cross-validation object
grid_search_best_params_ = bestDict['grid_search_best_params_']
#best_params = grid_search.best_params_
best_pipe.set_params(**grid_search_best_params_) # Set the best hyperparameters
# fetch the original input data
X_combined_df = bestDict['X_combined_df'] # name, raw input features, score
X_featuresOnly = X_combined_df.drop(columns=['name','score']) # features only for fit
nameCol = X_combined_df['name']
y_df = bestDict['y_df']
nameAndClassID_df = pd.concat([nameCol,y_df],axis=1)
nameAndClassID_df.columns = ['name','classID']
name_classID_features_df = pd.concat([nameAndClassID_df,X_featuresOnly],axis=1) # name, classID, all feature cols
unique_names_classIDs = name_classID_features_df.drop_duplicates(subset=['name', 'classID']).copy().iloc[:,0:2] # unique names and matching classID cols
y = y_df.to_numpy(copy=False)
printSB('\nStarting ' + str(kRANDiters) + ' RAND iterations...')
d_RANDlist = [] # list of dicts that we will average into a aggregate RAND result
for RANDctr in range(kRANDiters):
if (RANDctr>0) and ((RANDctr+1)%5 == 0): printSB(' RAND ' + str(RANDctr+1) + '...')
temp_RAND, randSuffix = randomizeImageData(name_classID_features_df, unique_names_classIDs, randflag=1, verbose=0) # classID-randomized name, classID, all feature cols
nameAndClassID_RAND = temp_RAND[['name', 'classID']]
X_featuresOnly_RAND = temp_RAND.drop(columns=['name','classID']) # features only for fit
y_RAND = temp_RAND['classID'].to_numpy(copy=False)
grid_search.fit(X_featuresOnly_RAND, y_RAND) # repeat the gridsearch on target-randomized sets
accuracy = grid_search.best_score_ # test_accuracy
best_pipe = grid_search.best_estimator_
best_pipe.set_params(**grid_search.best_params_) # Set the best hyperparameters
distances = best_pipe.decision_function(X_featuresOnly_RAND) # by calling best_pipe's decision_function, X_featuresOnly will flow thru the entire scale->PCA->LDA pipeline
distances_series = pd.Series(distances, name='score') # Convert the NumPy array to a Series for seamless appending
distances_df = pd.DataFrame(distances_series)
d_RAND = {}
d_RAND['nameAndClassID_RAND'] = nameAndClassID_RAND
d_RAND['distances_df'] = distances_df
d_RAND['accuracy'] = accuracy
d_RANDlist.append(d_RAND)
#printSB('RAND:', d_best_RAND['Pstr'], d_best_RAND['AUCstr'], d_best_RAND['CM_logbase10str'])
### here we have several d_RANDs in the d_RANDlist: average the results into a new d_best_RAND dict
# average scores by name, then compute stats on these aggregate avg scores
nameColRAND = d_RANDlist[0]['nameAndClassID_RAND']['name']
name_and_distance_RAND = pd.concat([nameColRAND,d_RANDlist[0]['distances_df']], axis=1) # first one
distances_df_RANDconcat = name_and_distance_RAND
for j in range(1,len(d_RANDlist)):
nameColRAND = d_RANDlist[j]['nameAndClassID_RAND']['name']
distances = d_RANDlist[j]['distances_df']
name_and_distances_df_RAND = pd.concat([nameColRAND,distances], axis=1)
distances_df_RANDconcat = pd.concat([distances_df_RANDconcat,name_and_distances_df_RAND],axis=0) # append new rows
distances_df_RANDconcat.reset_index(inplace = True, drop = True)
# group all scores by name, then calculate the mean of these grouped scores
grouped = distances_df_RANDconcat.groupby('name')['score'].mean()
distances_df_RAND_avg = pd.DataFrame(grouped).reset_index()
# for some reason distances on RAND data often have a large range: do we want to scale these for comparison with non-RAND?
scaler = StandardScaler()
distances_df_RAND_avg['score'] = scaler.fit_transform(distances_df_RAND_avg[['score']]) # std scale only the 'score' col
# stats on distances_df_RAND_avg (now original classes need not be 0 or 1, the classID column was reset to 0 & 1 for simplicity)
class0 = distances_df_RAND_avg.loc[y == 0]['score']
class1 = distances_df_RAND_avg.loc[y == 1]['score']
_, pValue_RAND = stats.ttest_ind(class0, class1, equal_var = True, alternative='less') #run independent 2 sample T-Test
Pstr_RAND = formatPstring(pValue_RAND)
# AUC
y_pred_RAND = np.where(distances_df_RAND_avg['score'] > 0, 1, 0) # so convert it to a predicted 0 or 1 class (need not correspond to classPair), depending on <0 vs >0
# y_actual = names_classIDs_X_reduced_df['classID']
y_actual = y
AUC_RAND = metrics.roc_auc_score(y_actual, y_pred_RAND)
if AUC_RAND == 1.0:
AUCstr_RAND = 'AUC = 1.0'
else:
AUCstr_RAND = 'AUC = {:.2f}'.format(AUC_RAND)
# accuracy
accuracy_RAND = metrics.accuracy_score(y_actual, y_pred_RAND)
if accuracy_RAND == 1.0:
accuracyStr_RAND = 'accuracy=1.0'
else:
accuracyStr_RAND = 'accuracy={:.2f}'.format(accuracy_RAND)
# f1 score
f1_RAND = metrics.f1_score(y_actual, y_pred_RAND)
# composite metric
CM_RAND, CMstr_RAND = calcCM(pValue_RAND, AUC_RAND, accuracy_RAND, p_value_logbase=10) # Higher values of p_value_logbase emphasize AUC over P value; p_value_logbase=0 uses AUC only
# Merging the DataFrames based on the 'name' column
out_df_RAND_avg = nameAndClassID_df.merge(distances_df_RAND_avg[['name', 'score']], on='name', how='left')
d = {}
d['pValue_RAND'] = pValue_RAND
d['Pstr_RAND'] = Pstr_RAND
d['AUC_RAND'] = AUC_RAND
d['AUCstr_RAND'] = AUCstr_RAND
d['accuracy_RAND'] = accuracy_RAND
d['accuracyStr_RAND'] = accuracyStr_RAND
d['CM_RAND'] = CM_RAND
d['CMstr_RAND'] = CMstr_RAND
d['out_df_RAND_avg'] = out_df_RAND_avg
return d
# custom print that wraps text in tokens that SB detects allows it to pass through for printing when soft-suspend is on. Useful to suppress unwanted warning messages from deep within certain libraries
def printSB_OLD(*args):
# must match in SB:
kSoftSuspendToken_start = "|||SB"
#kSoftSuspendToken_end = "SB|||"
if not args:
print() # Print an empty line if no arguments are provided
return
modified_args = list(args) # Convert tuple to list for easy modification
# Modify the first argument
modified_args[0] = f"{kSoftSuspendToken_start}{modified_args[0]}"
# Modify the last argument
#modified_args[-1] = f"{modified_args[-1]}{kSoftSuspendToken_end}"
# Print all arguments
print(*modified_args)
def printSB(*args, **kwargs):
kSoftSuspendToken_start = "|||SB"
#kSoftSuspendToken_end = "SB|||"
end = kwargs.get('end', '')
# Handle the 'sep' parameter if provided in kwargs
sep = kwargs.get('sep', ' ')
# Convert all positional arguments to strings and join them with spaces
args_str = ' '.join(str(arg) for arg in args)
args_str = sep.join(str(arg) for arg in args)
# Add the end character(s)
result = args_str + end # can write this to disk
print(kSoftSuspendToken_start + result)
return
"""
if not args:
print(**kwargs) # Print an empty line if no arguments are provided, respecting kwargs
# printStr = **kwargs
else:
modified_args = list(args) # Convert tuple to list for easy modification
modified_args[0] = f"{kSoftSuspendToken_start}{modified_args[0]}"
# Extract 'end' from kwargs if present, otherwise use default '\n'
end = kwargs.pop('end', '\n')
# Print all arguments, using the specified or default 'end'
print(*modified_args, end=end, **kwargs)
#if end == '\n': print(kSoftSuspendToken_end) # close bracket, this traps end=??? forms
"""
def printARC(*args, **kwargs):
# if printARC_fOut is a valid path will also write the string(s) to disk
end = kwargs.get('end', '')
# Handle the 'sep' parameter if provided in kwargs
sep = kwargs.get('sep', ' ')
# Convert all positional arguments to strings and join them with spaces
args_str = ' '.join(str(arg) for arg in args)
args_str = sep.join(str(arg) for arg in args)
# Add the end character(s)
result = args_str + end # can write this to disk
print(result)
try:
if printARC_fOut != '':
# Append mode - adds to existing content
with open(printARC_fOut, 'a') as file:
file.write(result+'\n')
except NameError:
print(f"DEBUG: printARC_fOut ({printARC_fOut}) does not exist")
return
def twoComponentScatterplot(names_classIDs_X_reduced_df, classPair, lambdas, normalize, axisLabels, CWT_type, graph_outPDF, subjectEmphasisList=[''], clf=None, ax=None, extraTitleLine=''):
# names_classIDs_X_reduced_df is a >=2 col df of the 2 components eg PC1,PC2, etc. Can have more than 2 cols, only the 1st 2 are plotted
# eg axisLabels = ['PC1',PC2']
# if an SVC was done pass clf to plot the decision boundary
# pass an existing subplot in ax to add to a graph instead
#printSB('twoComponentScatterplot:names_classIDs_X_reduced_df:\n',names_classIDs_X_reduced_df)
### 2-panel graph
axIsNone = (ax==None)
targets = 2
colors = ['g','r']
legend = ['class ' + str(classPair[0]), 'class ' + str(classPair[1])]
xLabel = axisLabels[0]
yLabel = axisLabels[1]
plotTitleFontSize = 14
axisLabelFontSize = 13
tickLabelFontSize = 11
statsStrFontSize = 12
nameLabelsFontSize = 7
legendFontSize = 12
dotplotMarkerSize = 40
# MANOVA P
# Filter data based on classID (previously reset to 0 and 1 regardless of actual) and select the first two numerical columns
group_0_data = names_classIDs_X_reduced_df[names_classIDs_X_reduced_df['classID'] == classPair[0]].iloc[:, 2:4].values
group_1_data = names_classIDs_X_reduced_df[names_classIDs_X_reduced_df['classID'] == classPair[1]].iloc[:, 2:4].values
# Convert to numpy arrays (if not already)
grp0 = np.array(group_0_data)
grp1 = np.array(group_1_data)
_,MANOVA_Pstr = stats_MANOVA_np(grp0,grp1)
if ax == None:
# plot formats:
nSubplotRows = 1
nSubplotCols = 2
figHeight_inches = 5
figWidth_inches = 2.07*figHeight_inches # match aspect ratio of canvas for nice display
axIndex = 1
fig, axs = plt.subplots(nSubplotRows,nSubplotCols,figsize=(figWidth_inches, figHeight_inches))
fig.tight_layout(pad=1.0) # adjust spacing between subplots; https://www.geeksforgeeks.org/how-to-set-the-spacing-between-subplots-in-matplotlib-in-python/
### avg spectra
ax=plt.subplot(nSubplotRows, nSubplotCols, 1) # left panel
ax.set_title('Averaged Spectra',size=plotTitleFontSize)
ax.set_xlabel('Wavelength (nm)', size= axisLabelFontSize)
if kNormalizeSpectralOverlay:
ax.set_ylabel('Normalized Intensity', size= axisLabelFontSize)
else:
ax.set_ylabel('Intensity', size= axisLabelFontSize)
ax.tick_params(axis='both', labelsize=tickLabelFontSize, direction='out', length=10, bottom=True, left=True) # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.tick_params.html
#ax.set_xticks([0,1],['0 (n=' + str(score_df[score_df['classID'] == 0].shape[0]) + ')','1 (n=' + str(score_df[score_df['classID'] == 1].shape[0]) + ')']) # classID labels and n's
# Define the tick locations
# Get the numerical columns for plotting
ax.grid(axis='y', color='0.8', linewidth=1.0)
#ax.axhline(0, color='blue')
#ax.legend(['Class 0', 'Class 1'])
columns_to_plot = imageData_avg.columns[firstLambdaCol:].to_numpy().astype(np.float64)
tick_locations = np.arange(400, columns_to_plot.max() + 50, 50) # Start at 400, end at max value + 50, with a step of 50
# Set the tick locator
ax.set_xlim([tick_locations.min(), tick_locations.max()]) # Set limits to encompass ticks
# Iterate over each unique name in the DataFrame
for name in imageData_avg['name'].unique():
# Filter the DataFrame for the current name
name_data = imageData_avg[imageData_avg['name'] == name]
spectrum = name_data.iloc[0,firstLambdaCol:firstLambdaCol+nLambdas].to_numpy().astype(np.float64)
if normalize: spectrum /= spectrum.max()
classID = int(name_data['classID'].iloc[0])
linewidth=1.0 # default
if name in subjectEmphasisList: linewidth+=1.5 # plot with heavier stroke
ax.plot(lambdas, spectrum, color=colors[classID],linewidth=linewidth)
leg = ax.legend([legend[0], legend[1]],
loc='upper right',
fontsize=11,
facecolor='white',
edgecolor='lightgray', # Remove legend outline
handlelength=0,
handletextpad=0, # Remove colored lines/handles
labelcolor=[colors[0],colors[1]])
# this is the only way to suppress the tiny handles in front of the legend labels:
for line in leg.get_lines():
line.set_linewidth(0.0)
axIndex += 1
ax = plt.subplot(nSubplotRows, nSubplotCols, axIndex)
### 2 component scatterplot
if extraTitleLine != '': extraTitleLine = '\n' + extraTitleLine
if clf == None:
graphTitle = 'First 2 components (' + CWT_type + ')' + extraTitleLine
ax.set_title(graphTitle, fontsize = plotTitleFontSize-1)
else:
graphTitle = 'First 2 components (' + CWT_type + ')\n(+nuSVC decision boundary)'
ax.set_title(graphTitle, fontsize = plotTitleFontSize-2) + extraTitleLine
ax.set_xlabel(xLabel, fontsize = axisLabelFontSize)
ax.set_ylabel(yLabel, fontsize = axisLabelFontSize)
ax.tick_params(axis='both', labelsize=tickLabelFontSize, direction='out', length=10, bottom=True, left=True) # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.tick_params.html
X_reduced_df = names_classIDs_X_reduced_df.iloc[:, 2:4] # 2 reduced components
X_reduced_np = X_reduced_df.values
for n in range(2): # classID's were reset to 0 & 1 for simplicity
# indicesToKeep = (imageData_avg_df['classID'] == n)
indicesToKeep = (names_classIDs_X_reduced_df['classID'] == n)
#printSB('indicesToKeep @ ' + str(n),'\n',indicesToKeep)
#printSB('X_reduced_df:\n',X_reduced_df)
color = colors[n]
ax.scatter(X_reduced_df.loc[indicesToKeep, xLabel]
, X_reduced_df.loc[indicesToKeep, yLabel]
, c = color
, s = dotplotMarkerSize
, alpha=1)
# marker emphasis
# for index, row in imageData_avg_df.iterrows():
for index, row in names_classIDs_X_reduced_df.iterrows():
if row['name'] in subjectEmphasisList: ax.scatter(X_reduced_df.iloc[index,0], X_reduced_df.iloc[index,1], s=110, c=colors[0] if row['classID'] == 0 else colors[1])
# Auto-adjust x and y axis limits
x_margin = (X_reduced_np[:, 0].max() - X_reduced_np[:, 0].min()) * 0.05 # 5% margin on each side
y_margin = (X_reduced_np[:, 1].max() - X_reduced_np[:, 1].min()) * 0.05
ax.set_xlim(X_reduced_np[:, 0].min() - x_margin, X_reduced_np[:, 0].max() + x_margin)
ax.set_ylim(X_reduced_np[:, 1].min() - y_margin, X_reduced_np[:, 1].max() + y_margin)
#ax.legend(legend,fontsize = legendFontSize)
# all the options are to make the box smaller
ax.legend(
legend,
fontsize = legendFontSize-1,
loc='upper left',
#bbox_to_anchor=(0, 1),
#bbox_transform=ax.transAxes,
markerscale=0.5, # Smaller marker (0.5 is half the default size)
handlelength=1.0, # Shorter connecting line
labelspacing=0.2, # Tighter vertical spacing
borderpad=0.2, # Tighter border padding
handletextpad=0.2 # Less padding between marker/line and text
)
if clf != None: # plot decision boundary
# Create a meshgrid to visualize the decision boundary
h = .02 # Step size in the mesh
xlim = ax.get_xlim()
ylim = ax.get_ylim()
# rand is
x_min, x_max = X_reduced_np[:, 0].min() - 1, X_reduced_np[:, 0].max() + 1
y_min, y_max = X_reduced_np[:, 1].min() - 1, X_reduced_np[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, h),np.arange(y_min, y_max, h))
# Predict the class labels for the meshgrid points
Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
# Plot the decision boundary and data points
#ax.contourf(xx, yy, Z, cmap=plt.cm.coolwarm, alpha=0.3)
from matplotlib.colors import ListedColormap
cmap = ListedColormap(['g', 'r']) # Green for class 0, red for class 1
ax.contourf(xx, yy, Z, cmap=cmap, alpha=0.1)
# print stats
x_min, x_max = ax.get_xlim()
y_min, y_max = ax.get_ylim()
x_range = x_max - x_min
y_range = y_max - y_min
ax.annotate(text=MANOVA_Pstr, xy=(x_min+x_range/2, y_min+y_range*0.5), xycoords='data', textcoords='offset points',fontsize=12, horizontalalignment='center')
if kLabelDatapoints:
# xRange = plotDf[xLabel].max() - plotDf[xLabel].min()
# for i in range(len(imageData_avg_df)):
for i in range(len(names_classIDs_X_reduced_df)):
# classID = imageData_avg_df.iloc[i].at['classID']
# name = imageData_avg_df.iloc[i].at['name']
classID = names_classIDs_X_reduced_df.iloc[i].at['classID']
name = names_classIDs_X_reduced_df.iloc[i].at['name']
PCx = X_reduced_df.iloc[i,0]
PCy = X_reduced_df.iloc[i,1]
xytext_xoffset = 4
if name in subjectEmphasisList: xytext_xoffset += 2
ax.annotate(name, xy=(PCx, PCy), xycoords='data',xytext=(xytext_xoffset,-2), textcoords='offset points',fontsize=nameLabelsFontSize) # https://matplotlib.org/stable/users/explain/text/annotations.html
if axIsNone: # we only want to save if not part of an existing graph passed in ax
plt.savefig(graph_outPDF, bbox_inches='tight')
# printSB('\n' + xLabel + ' vs ' + yLabel +' subject means graph saved to: ' + graph_outPDF)
plt.close() # clear for next plot else we have old annotations persisting
def binarize01(y): # converts all values in y <= 0 to 0, otherwise 1, returns a copy
y01 = y.copy()
# convert to 0 v 1
y01[y01 <= 0] = 0
y01[y01 > 0] = 1
y01 = y01.astype(int)
return y01
def check_type(obj):
if isinstance(obj, np.ndarray):
return "NumPy Array"
elif isinstance(obj, pd.DataFrame):
return "pandas DataFrame"
elif isinstance(obj, pd.Series):
return "pandas Series"
elif isinstance(obj, list):
return "Python List"
elif isinstance(obj, dict):
return "Python Dictionary"
else:
return type(obj).__name__
def round_to_N_sig(x, N=2):
"""
Round the input to N significant digits.
Parameters:
x : numpy array or scalar
The input value(s) to be rounded
N : int
The number of significant digits to round to
Returns:
numpy array or scalar
The input rounded to N significant digits eg.
"""
sign_x = np.sign(x)
x = np.abs(x) # we'll operate on +ve values only then adjust
with np.errstate(divide='ignore', invalid='ignore'):
# exponent = np.floor(np.log10(np.abs(x)))
exponent = np.floor(np.log10(x)) # x will now always be +ve
mantissa = x / (10 ** exponent)
rounded_mantissa = np.round(mantissa, N - 1)
r = rounded_mantissa * (10 ** exponent)
return sign_x * r
def round_to_N_sig_str(x, N=2):
"""
as above but returns a string representation
"""
r = round_to_N_sig(x=x, N=N)
return f'{r:g}'
def format_tick_labels(ax):
# formats numerical tick labels for subplots using the 'g' auto-format
# For y-axis
yticks = ax.get_yticks()
ax.yaxis.set_major_formatter(plt.FormatStrFormatter('%g'))
# For x-axis
xticks = ax.get_xticks()
ax.xaxis.set_major_formatter(plt.FormatStrFormatter('%g'))
### dimensionality reduction using a CNN2D auto encoder
# Define the autoencoder model
def build_autoencoder(input_shape, n_components):
# Encoder
encoder_input = layers.Input(shape=input_shape)
x = layers.Conv2D(32, (3, 3), activation='relu', padding='same')(encoder_input)
x = layers.MaxPooling2D((2, 2), padding='same')(x)
x = layers.Conv2D(16, (3, 3), activation='relu', padding='same')(x)
x = layers.MaxPooling2D((2, 2), padding='same')(x)
x = layers.Conv2D(8, (3, 3), activation='relu', padding='same')(x)
x = layers.MaxPooling2D((2, 2), padding='same')(x)
x = layers.Flatten()(x)
encoded = layers.Dense(n_components, activation='linear')(x)
# Decoder
x = layers.Dense(8 * (input_shape[0] // 8) * (input_shape[1] // 8), activation='relu')(encoded)
x = layers.Reshape((input_shape[0] // 8, input_shape[1] // 8, 8))(x)
x = layers.Conv2D(8, (3, 3), activation='relu', padding='same')(x)
x = layers.UpSampling2D((2, 2))(x)
x = layers.Conv2D(16, (3, 3), activation='relu', padding='same')(x)
x = layers.UpSampling2D((2, 2))(x)
x = layers.Conv2D(32, (3, 3), activation='relu', padding='same')(x)
x = layers.UpSampling2D((2, 2))(x)
decoded = layers.Conv2D(input_shape[2], (3, 3), activation='sigmoid', padding='same')(x)
# Autoencoder
autoencoder = models.Model(encoder_input, decoded)
encoder = models.Model(encoder_input, encoded)
return autoencoder, encoder
# Function to perform dimensionality reduction
def reduce_dimensionality_2DCNN_AE(data, n_components, epochs=50, batch_size=32, verbose=1):
N, H, W, C = data.shape
input_shape = (H, W, C)
# Build and compile the autoencoder
autoencoder, encoder = build_autoencoder(input_shape, n_components)
autoencoder.compile(optimizer='adam', loss='mse')
# Split the data into training and validation sets
X_train, X_val = train_test_split(data, test_size=0.2, random_state=42)
# Train the autoencoder
history = autoencoder.fit(X_train, X_train, epochs=epochs, batch_size=batch_size,
shuffle=True, validation_data=(X_val, X_val), verbose=verbose)
# Use the encoder to get the reduced representation
reduced_data = encoder.predict(data)
# Evaluate the model
reconstructed_data = autoencoder.predict(data)
mse = mean_squared_error(data.flatten(), reconstructed_data.flatten())
return reduced_data, mse, history
def estimatedTimeRemaining(start_time, nIters, totalIters):
# save start_time = datetime.now() at start of long operatino and pass here
if (nIters <= 0) or (nIters > totalIters): return '?'
current_time = datetime.now()
elapsed_time = current_time - start_time
elapsed_seconds = elapsed_time.total_seconds()
proportion_complete = nIters / totalIters
# total_seconds = elapsed_seconds / proportion_complete # estimated total seconds
total_seconds = (totalIters/nIters)*elapsed_seconds # estimated total run time
# seconds_remaining = (totalIters-nIters)*total_seconds/totalIters
seconds_remaining = total_seconds-elapsed_seconds
if seconds_remaining < 60: # < 1hr
return '<1 min'
elif seconds_remaining < 3600: # < 1hr
return f'{int(seconds_remaining/60+0.5)} min'
elif seconds_remaining < 3600*10: # < 10hr
return f'{seconds_remaining/3600:.1f} hrs'
else: # >10hr
return f'{seconds_remaining/3600:.0f} hrs'
def remove_proportional_baseline(S,B):
# S is 1D array containing unknown spectrum X + p*baseline spectrum B
# finds p and recovers X
def objective(p):
return np.sum((S - p * B) ** 2)
def find_optimal_proportion(S, B):
result = minimize( objective,
x0=0.5,
bounds=[(0, None)],
method='SLSQP', # https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
options={
'ftol': 1e-30,
'eps': 1e-30,
'maxiter': 10000,
'disp': False
}
)
return result.x[0]
optimal_proportion = find_optimal_proportion(S, B)
X_recovered = S - optimal_proportion * B
return X_recovered
def RandomForestClassifier_loo_grid(feature_matrix_np, target_np, n_jobs=-1, swarmplot=False, q=15):
# feature_matrix_np shape (N,P) where N is instances, P is predictors eg PC1, PC2, etc
# target_np: class labels eg 0, 1
# maybe need to pass n_jobs=1 if already in an MP enivronment?
# q: Number of instances above which we switch from LeaveOneOut to RepeatedStratifiedKFold cross-val
# returns instance_probas = class1 probability
# Define parameter grid
param_grid = {
'n_estimators': [50, 100, 200],
'max_depth': [None, 10, 20, 30],
'min_samples_split': [2, 5, 10],
'min_samples_leaf': [1, 2, 4],
'max_features': ['sqrt', 'log2', None] # Remove 'auto', use valid options
}
# Choose cross-validation strategy for GridSearch
if len(target_np) < q:
cv_grid = LeaveOneOut()
cv_pred = LeaveOneOut()
else:
cv_grid = RepeatedStratifiedKFold(n_splits=5, n_repeats=5, random_state=42)
cv_pred = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
# Initialize components and perform GridSearch with LOO
base_clf = RandomForestClassifier(random_state=42)
# Initialize GridSearchCV with error handling
grid_search = GridSearchCV(
estimator=base_clf,
param_grid=param_grid,
scoring='accuracy',
cv=cv_grid,
n_jobs=n_jobs,
error_score=np.nan
)
# Fit and get predictions
with np.errstate(invalid='ignore'): # ignore failed fits
grid_search.fit(feature_matrix_np, target_np)
# Get cross-validated probabilities
cv_probas = cross_val_predict(
grid_search.best_estimator_,
feature_matrix_np,
target_np,
cv=cv_pred,
method='predict_proba'
)
instance_probas = cv_probas[:, 1] # class1 probas
predicted_classes = (instance_probas > 0.5).astype(int)
# predicted_accuracy = accuracy_score(target_np, predicted_classes) # More likely to be optimistically biased (overfit)
mean_accuracy = grid_search.best_score_ # Obtained from cross-validation during grid search • Averages the accuracy scores across all CV folds • Represents the mean performance across different train-test splits • Generally more realistic estimate of model performance
if swarmplot:
# Create swarm plot
plt.figure(figsize=(7, 4))
sns.swarmplot(x=target_np, y=instance_probas, hue=target_np, palette={0: 'green', 1: 'red'}, legend=False)
plt.axhline(y=0.5, color='gray', linestyle='--', alpha=0.5)
plt.xlabel('Class')
plt.ylabel('Predicted Probability')
plt.title(f'Predicted Probabilities by Class (mean accuracy {mean_accuracy:.2f})')
#plt.show()
return mean_accuracy, predicted_classes, instance_probas, grid_search.best_params_, plt
else:
return mean_accuracy, predicted_classes, instance_probas, grid_search.best_params_, None
def array_to_line_distances(X, w, b):
"""Calculate signed Euclidean distances from points to decision boundary.
Points where (X·w + b) > 0 get positive distances (class 1 side)
Points where (X·w + b) < 0 get negative distances (class 0 side)
"""
# Calculate the decision function values
decision_values = np.dot(X, w) + b
# Convert to distances by normalizing by the weight vector norm
distances = decision_values / np.linalg.norm(w)
# Optionally reverse signs if you want class 0 side to be positive:
# distances = -distances
return distances
"""
# sample code for scatterplot and perpendiculars with distances:
def plot_with_perpendiculars(X, y, w, b, distances):
plt.figure(figsize=(8, 8)) # Square figure
# Plot points
plt.scatter(X[y == 0][:, 0], X[y == 0][:, 1],

stys_WTv5xFAD.py, under CC-BY-4.0 · at the source

Overview

Authors: Shigeki Tsutsui^1,2, Anastasiia A Stepanchuk¹, Julian P Stys^1,3, Stefanie A G Black^1,2, George W Templeton², Russell Greiner⁴, Peter K Stys^1,2

ORCID iDs: Anastasiia A Stepanchuk, Stefanie A G Black, George W Templeton, Russell Greiner, Peter K Stys

Hotchkiss Brain Institute, Department of Clinical Neurosciences, Cumming School of Medicine, University of Calgary, Calgary, Alberta, Canada
Amira Medical Technologies Inc., Calgary, Alberta, Canada
NovaSoft Interactive, Calgary, Alberta, Canada
Department of Computing Science, Alberta Machine Intelligence Institute, University of Alberta, Edmonton, Canada

Institutions: University of Calgary (Canada); Alberta University of the Arts (Canada); University of Alberta (Canada)

Journal: Journal of Alzheimer's disease : JAD, volume 112, issue 2, pages 933-946

Dates: received 2 October 2025; accepted 19 April 2026; published online 3 June 2026; in print July 2026

Type: Research article · Language: English

License: CC BY-NC

Identifiers: DOI 10.1177/13872877261453512 · PMID 42231859 · PMCID PMC13334060 · OpenAlex W7163317707

Open access: hybrid, a free copy (OpenAlex)

Status: code verified

Categories: histology / microscopy (modality), human (organism), mouse (organism), Alzheimer's / dementia (population), methods / tools (subfield)

Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, fMRI & imaging

Keywords: Alzheimer's disease, conformationally-sensitive amyloid probes, machine learning, peripheral blood mononuclear cells, spectral confocal microscopy, wavelet

MeSH: Alzheimer Disease*, Leukocytes, Mononuclear*, Machine Learning*, Amyloid beta-Peptides, Animals, Female, Humans, Mice, Mice, Transgenic, Microscopy, Confocal, Spectrometry, Fluorescence, Support Vector Machine (* major topic)

Topic: Alzheimer's disease research and treatments (Physiology, Medicine), according to OpenAlex

Funding: Krembil Foundation; W. Garfield Weston Foundation; Brain Canada/Alzheimer’s Society Research Program; Accelerating Innovations into CarE; Alzheimer's Drug Discovery Foundation

Citations: not cited yet (Europe PMC); 68 references in the paper

Abstract

Background: Alzheimer's disease (AD) is the most common cause of dementia whose prevalence is projected to increase significantly in the coming decades. The recent advent of disease modifying therapies is a welcome development; however, it is also now apparent that early treatment maximizes the benefits of these drugs. Therefore, it is important to develop reliable methods of disease detection, preferably from an easily accessible matrix such as blood.

Objective: To develop a method for detecting AD from circulating white blood cells using spectral confocal microscopy.

Methods: Using K114-stained wild type and 5xFAD transgenic mouse cortical sections as proof-of-principle, spectral imaging of K114 fluorescence coupled with a signal processing/machine learning pipeline (spectral wavelet decomposition, dimensionality reduction, support vector machine classifier) can reliably distinguish non-plaque background parenchyma in the two strains. We then performed immunoprecipitation of Aβ from peripheral blood mononuclear cells (PBMCs) obtained from non-neurological controls and histopathologically-proven AD cases. We spectrally imaged the immunobeads labeled with K114, then used similar machine learning methods to classify control versus AD samples.

Results: Normal-appearing non-plaque 5xFAD background was reliably distinguished from wild type mouse brain. We could also classify AD with a high degree of reliability (area under the receiver operating curve = 0.95, p = 6.1e-5) and predict neuropathological scores from these blood elements (R = 0.89).

Conclusions: Our spectral imaging method, together with automated machine learning analysis of spectral micrographs, using readily obtainable PBMCs from blood, represents a potentially useful approach for detection of AD in living subjects.

Reproduced under the paper's license (CC BY-NC), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 4 matches between paragraphs and lines of code.

Zenodo 18217440

License: CC-BY-4.0

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Languages: Python (1)

Size: 1 file, 1 script

Software Heritage: not checked

Found in: the text, “Machine learning methods and statistics for spec”

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: Keras (1 file), Matplotlib (1 file), NumPy (1 file), pandas (1 file), PyWavelets (1 file), scikit-learn (1 file), SciPy (1 file), seaborn (1 file), statsmodels (1 file), TensorFlow (1 file), UMAP (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)

27 September 2026: the link answers (HTTP 200)

1 file

stys_WTv5xFAD.py, Python, 4,178 lines, 4 matches

At the source: doi.org/10.5281/zenodo.18217440

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
1 script, each with its path and the digest of its content;
4 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

Publisher: n/a → IOS Press

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 6 keywords, 12 MeSH terms, 5 funders, 65 references.

Cite

This paper

Tsutsui, S., Stepanchuk, A. A., Stys, J. P., Black, S. A. G., Templeton, G. W., Greiner, R., & Stys, P. K. (2026). Fluorescence spectroscopy and machine learning methods for detection of Alzheimer's disease from circulating white blood cells. Journal of Alzheimer's disease : JAD, 112(2), 933-946. https://doi.org/10.1177/13872877261453512

BibTeX

@article{tsutsui2026fluorescence, author = {Tsutsui, Shigeki and Stepanchuk, Anastasiia A and Stys, Julian P and Black, Stefanie A G and Templeton, George W and Greiner, Russell and Stys, Peter K}, title = {{Fluorescence spectroscopy and machine learning methods for detection of Alzheimer's disease from circulating white blood cells}}, journal = {Journal of Alzheimer's disease : JAD}, year = {2026}, month = jun, volume = {112}, number = {2}, pages = {933--946}, publisher = {IOS Press}, issn = {1387-2877}, doi = {10.1177/13872877261453512}, url = {https://doi.org/10.1177/13872877261453512}, pmid = {42231859}, pmcid = {PMC13334060} }

RIS

TY - JOUR AU - Tsutsui, Shigeki AU - Stepanchuk, Anastasiia A AU - Stys, Julian P AU - Black, Stefanie A G AU - Templeton, George W AU - Greiner, Russell AU - Stys, Peter K TI - Fluorescence spectroscopy and machine learning methods for detection of Alzheimer's disease from circulating white blood cells T2 - Journal of Alzheimer's disease : JAD J2 - J Alzheimers Dis PY - 2026 DA - 2026/06/03 VL - 112 IS - 2 SP - 933 EP - 946 SN - 1387-2877 PB - IOS Press DO - 10.1177/13872877261453512 UR - https://doi.org/10.1177/13872877261453512 LA - en ER -

CSL-JSON

{ "id": "10.1177/13872877261453512", "type": "article-journal", "title": "Fluorescence spectroscopy and machine learning methods for detection of Alzheimer's disease from circulating white blood cells", "container-title": "Journal of Alzheimer's disease : JAD", "author": [ { "family": "Tsutsui", "given": "Shigeki" }, { "family": "Stepanchuk", "given": "Anastasiia A" }, { "family": "Stys", "given": "Julian P" }, { "family": "Black", "given": "Stefanie A G" }, { "family": "Templeton", "given": "George W" }, { "family": "Greiner", "given": "Russell" }, { "family": "Stys", "given": "Peter K" } ], "container-title-short": "J Alzheimers Dis", "volume": "112", "issue": "2", "page": "933-946", "DOI": "10.1177/13872877261453512", "PMID": "42231859", "PMCID": "PMC13334060", "ISSN": "1387-2877", "publisher": "IOS Press", "URL": "https://doi.org/10.1177/13872877261453512", "language": "en", "issued": { "date-parts": [ [ 2026, 6, 3 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41592-026-03057-2 [code]: CREsted: modeling genomic and synthetic cell-type-specific enhancers across tissues and species.
Journal: Nature methods
In common: PyWavelets, Keras, UMAP, 8 other tools, methods / tools, mouse
[2] doi:10.1038/s41467-026-72057-9 [code]: Sex-specific behavioral feedback modulates sensorimotor processing and drives flexible social behavior.
Journal: Nature communications
In common: PyWavelets, Keras, UMAP, 8 other tools
[3] doi:10.1038/s42003-026-10957-8 [code]: Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: Keras, UMAP, TensorFlow, 7 other tools, mouse
[4] doi:10.1093/nar/gkag706 [code]: scDifformer: diffusion-based post-training for virtual cell modeling across large-scale single-cell data.
Journal: Nucleic acids research
In common: Keras, UMAP, TensorFlow, 7 other tools
[5] doi:10.1038/s41398-026-04081-8 [code]: Functional system-specific brain aging across the Alzheimer's disease continuum.
Journal: Translational psychiatry
In common: Keras, TensorFlow, statsmodels, 6 other tools, Alzheimer's / dementia, 1 reference
[6] doi:10.1093/jnen/nlaf152 [code]: Clinical and pathologic correlations of machine learning quantification of Aβ deposits across 3 brain regions of decedents with Alzheimer disease.
Journal: Journal of neuropathology and experimental neurology
In common: Keras, TensorFlow, scikit-learn, 4 other tools, Alzheimer's / dementia, 2 references
[7] doi:10.64898/2026.05.06.26352540 [code]: Generating synthetic tau-PET scans in Alzheimer’s disease from MRI, blood biomarkers and demographics with deep learning
Journal: medRxiv (preprint)
In common: Keras, TensorFlow, seaborn, 5 other tools, Alzheimer's / dementia, 1 reference
[8] doi:10.1039/d6ra03343a [code]: A benchmark dataset and interpretable deep learning framework for drug-induced developmental neurotoxicity prediction.
Journal: RSC advances
In common: Keras, UMAP, TensorFlow, 6 other tools, methods / tools
[9] doi:10.1038/s41593-026-02376-z [code]: A framework for comparative analysis of human and mouse cortical neuron dendrites in corresponding brain regions.
Journal: Nature neuroscience
In common: Keras, UMAP, statsmodels, 6 other tools, mouse
[10] doi:10.1016/j.xcrm.2026.102766 [code]: A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Keras, UMAP, TensorFlow, 6 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.