Anushaoram/Dl_project

1import base642import io3import json4 5import numpy as np6import pandas as pd7import torch8from catboost import CatBoostClassifier9from matplotlib import pyplot as plt10from matplotlib import rcParams11from matplotlib.cm import ScalarMappable12from scipy.ndimage import gaussian_filter1d13from torch.utils.data import DataLoader14 15from model.DeepSpiro import DeepSpiro, MyDataset16 17config = {18    "font.family": 'Times New Roman',19    "axes.unicode_minus": False,20    "font.size": 18,21    "axes.labelsize": 30,22    "xtick.labelsize": 25,23    "ytick.labelsize": 25,24}25rcParams.update(config)26SPIRO_RECORD_SERIES_KEY = 'flow'27 28 29def load_config(config_path):30    with open(config_path, 'r') as f:31        config = json.load(f)32    return config33 34 35def smooth(data, sigma=1):36    smoothed_data = gaussian_filter1d(data, sigma=sigma)37    return smoothed_data38 39 40def compute_flow_volume_by_num_points(series, max_num_points, volume_scale=0.001, time_scale=0.01,41                                      max_interp_volume=6.58):42    volume = (series * volume_scale).astype(np.float32)43 44    flow = np.concatenate(([0.0], np.diff(volume) / time_scale))45 46    def right_pad_array(arr, pad_value, max_num_points):47        if len(arr) > max_num_points:48            return arr[:max_num_points]49        else:50            return np.pad(arr, (0, max_num_points - len(arr)), 'constant', constant_values=(pad_value,))51 52    padded_volume = right_pad_array(volume, 0, max_num_points)53    padded_flow = right_pad_array(flow, 0, max_num_points)54 55    monotonic_volume = np.maximum.accumulate(padded_volume)56    volume_interp_intervals = np.linspace(start=0, stop=max_interp_volume, num=max_num_points)57    flow_volume = np.interp(volume_interp_intervals, xp=monotonic_volume, fp=padded_flow, left=0, right=0)58 59    return volume, flow, flow_volume60 61 62def compute_fef(flow: np.ndarray, volume: np.ndarray, volume_max: float):63    flow_size = len(flow)64    assert flow_size == len(volume), 'Flow and Volume lengths do not match.'65    assert flow_size > 1, 'Flow should have more than one values'66    volumes_over_25 = volume >= (0.25 * volume_max)67    volumes_over_50 = volume >= (0.50 * volume_max)68    volumes_over_75 = volume >= (0.75 * volume_max)69    if not any(volumes_over_75):70        raise ValueError(f'Cannot find FEF75 in volume curve: {volume}')71 72    idx_25 = np.argmax(volumes_over_25)73    idx_50 = np.argmax(volumes_over_50)74    idx_75 = np.argmax(volumes_over_75)75    assert 0 <= idx_25 <= idx_50 <= idx_75 < flow_size76 77    fef25, fef50, fef75 = flow[[idx_25, idx_50, idx_75]]78    fef25_75 = flow[idx_25: (idx_75 + 1)].mean()79    return fef25, fef50, fef75, fef25_7580 81 82def calculate_index(row):83    flow = row['flow_volume']84    last_index = np.argmin(flow[5:])85    PEF_index = len(flow) - 1 - np.argmax(flow[::-1])86    flow_index = len(flow[:last_index + 1])87    index_25 = int(0.25 * flow_index)88    index_50 = int(0.50 * flow_index)89    index_75 = int(0.75 * flow_index)90 91    return PEF_index, index_25, index_50, index_75, last_index92 93 94def calculate_acceleration(row):95    if 'flow_volume' in row and isinstance(row['flow_volume'], np.ndarray):96        flow_volume = row['flow_volume']97        try:98            index_pef = int(row['index_pef'])99            start_index_25 = int(row['index_fef25'])100            end_index_50 = int(row['index_fef50'])101            end_index_75 = int(row['index_fef75'])102            last_index = int(row['last_index'])103            if not (index_pef < start_index_25 < end_index_50 < end_index_75):104                return np.nan, np.nan, np.nan, np.nan, np.nan, np.nan105        except ValueError:106            return np.nan, np.nan, np.nan, np.nan, np.nan, np.nan107 108        def calc_area_under_curve(flow_volume, start_index, end_index):109            slope = (flow_volume[end_index] - flow_volume[start_index]) / (end_index - start_index)110            intercept = flow_volume[start_index] - slope * start_index111            baseline = slope * np.arange(start_index, end_index + 1) + intercept112            differences = -(flow_volume[start_index:end_index + 1] - baseline)113            area_below = np.sum(differences[differences < 0]) * 0.01114            area_above = np.sum(differences[differences > 0]) * 0.01115            adjusted_area = area_below + area_above116            return adjusted_area117 118        area_pef_25 = calc_area_under_curve(flow_volume, index_pef, start_index_25)119        area_25_50 = calc_area_under_curve(flow_volume, start_index_25, end_index_50)120        area_50_75 = calc_area_under_curve(flow_volume, end_index_50, end_index_75)121        area_75 = calc_area_under_curve(flow_volume, end_index_75, last_index - 1)122        area_pef_75 = calc_area_under_curve(flow_volume, start_index_25, last_index - 1)123        area_p = area_pef_25 + area_25_50 - area_50_75 - area_75124        return area_pef_25, area_25_50, area_50_75, area_75, area_pef_75, area_p125 126 127def process_data(row):128    handle = row.copy()129    series = [int(v) for v in row[SPIRO_RECORD_SERIES_KEY].split(',')]130    series = np.array(series)131    series = smooth(series)132    volume, flow, flow_volume = compute_flow_volume_by_num_points(series, len(series))133    fef25, fef50, fef75, fef25_75 = compute_fef(flow, volume, volume.max())134    handle['blow_fef25'] = fef25135    handle['blow_fef50'] = fef50136    handle['blow_fef75'] = fef75137    handle['blow_fef25_75'] = fef25_75138    handle['flow_volume'] = flow_volume139    handle['volume'] = volume140    handle['flow'] = flow141    handle['series'] = series142    handle['PEF'] = row['pef'] if row['pef'] != '' else np.nan143    handle["FEV1"] = row["fev1"] if row["fev1"] != '' else np.nan144    handle["FVC"] = row["fvc"] if row["fvc"] != '' else np.nan145    return handle146 147 148def process_acceleration(row):149    row['index_pef'], row['index_fef25'], row['index_fef50'], row['index_fef75'], row['last_index'] = calculate_index(150        row)151    acceleration_pef_25, acceleration_25_50, acceleration_50_75, acceleration_75, acceleration_pef_75, acceleration_total = calculate_acceleration(152        row)153    acceleration_pef_25 = pd.Series(acceleration_pef_25)154    acceleration_25_50 = pd.Series(acceleration_25_50)155    acceleration_50_75 = pd.Series(acceleration_50_75)156    acceleration_75 = pd.Series(acceleration_75)157    acceleration_pef_75 = pd.Series(acceleration_pef_75)158    acceleration_total = pd.Series(acceleration_total)159    row['PEF_FEF25'] = acceleration_pef_25160    row['FEF25_FEF50'] = acceleration_25_50161    row['FEF50_FEF75'] = acceleration_50_75162    row['FEF75'] = acceleration_75163    row['PEF_FEF75'] = acceleration_pef_75164    row['TOTAL'] = acceleration_total165    return row166 167 168def preprocess_data(input_path, age, sex, smoke):169    if input_path.endswith('.xlsx'):170        df = pd.read_excel(input_path)171        if len(df) == 1:172            row = df.iloc[0]173            row = process_data(row)174            row = process_acceleration(row)175            processed_data = pd.Series(dtype='float64')176            processed_data['flow_volume'] = row['flow_volume']177            processed_data['PEF_FEF25'] = row['PEF_FEF25'].values[0]178            processed_data['FEF25_FEF50'] = row['FEF25_FEF50'].values[0]179            processed_data['FEF50_FEF75'] = row['FEF50_FEF75'].values[0]180            processed_data['FEF75'] = row['FEF75'].values[0]181            processed_data['PEF_FEF75'] = row['PEF_FEF75'].values[0]182            processed_data['TOTAL'] = row['TOTAL'].values[0]183            processed_data['AGE'] = age184            processed_data['SEX'] = sex185            processed_data['smoke'] = smoke186            processed_data['blow_ratio'] = 1 - (row['FEV1'] / row['FVC'])187            processed_data['fef25'] = row['blow_fef25']188            processed_data['fef50'] = row['blow_fef50']189            processed_data['fef75'] = row['blow_fef75']190            processed_data['FEV1'] = row['FEV1']191            processed_data['FVC'] = row['FVC']192            return processed_data193        else:194            AssertionError("Error: Only one row of data is supported.")195    else:196        AssertionError("Error: Unsupported file format.")197 198 199def load_spiro_encoder(device_str, model_path):200    device = torch.device(device_str if torch.cuda.is_available() else "cpu")201    model = DeepSpiro(202        in_channels=1,203        out_channels=32,204        n_len_seg=30,205        n_classes=2,206        device=device,207        verbose=False208    ).to(device)209    model.load_state_dict(torch.load(model_path, map_location=device, weights_only=True))210    return model211 212 213def load_cb_model(model_path):214    model = CatBoostClassifier()215    model.load_model(model_path)216    return model217 218 219def run_spiro_encoder(model, data, device):220    dataset = MyDataset([data['flow_volume']], 30)221    data_loader = DataLoader(dataset, batch_size=1, shuffle=False, num_workers=0)222    model.eval()223    predictions = []224    attention_weights = []225    all_input_x = []226    with torch.no_grad():227        for data, mask in data_loader:228            data = data.to(device)229            mask = mask.to(device)230            output = model(data, mask)231            probabilities = torch.softmax(output, dim=1)232            predictions.append(probabilities.cpu().numpy())233 234            temporal_attention_weights = model.temporal_attention.attention_weights.detach().cpu().numpy()235            temporal_attention_weights = np.squeeze(temporal_attention_weights, axis=-1)236            input_x1 = data.detach().cpu().numpy()237            input_x1 = input_x1.reshape((1, -1, 1))238            attention_weights_padded = np.zeros((1, data.shape[1], 1))239            attention_weights_padded[:, :temporal_attention_weights.shape[1], :] = temporal_attention_weights[:, :,240                                                                                   None]241            attention_weights_expanded = np.repeat(attention_weights_padded, 30, axis=1)242            attention_weights.append(attention_weights_expanded)243            all_input_x.append(input_x1)244    return predictions, attention_weights, all_input_x245 246 247def run_spiro_explainer(model, data, threshold, spiro_encoder_original_result, attention_weights, all_input_x,248                        is_show=True):249    spiro_encoder_result = spiro_encoder_original_result[0][0][1]250    buf = plt_attention(251        all_input_x, attention_weights, data['fef25'], data['fef50'],252        data['fef75'], data['FEV1'], data['FVC'], is_show=is_show253    )254    image_base64 = base64.b64encode(buf.read()).decode('utf-8')255 256    data['copd_detection'] = spiro_encoder_result257 258    X_pred = [[data['AGE'], data['SEX'], data['smoke'], data['blow_ratio'], data['copd_detection']]]259    probabilities = model.predict_proba(X_pred)260    detection = (probabilities[0][1] >= threshold).astype(int)261    return detection, image_base64262 263 264def run_spiro_predictor(model, data):265    X_pred = [266        [data['PEF_FEF25'], data['FEF25_FEF50'], data['FEF50_FEF75'], data['FEF75'], data['PEF_FEF75'], data['TOTAL'],267         data['copd_detection']]]268    probabilities = model.predict_proba(X_pred)269    return probabilities270 271 272def find_closest_x(y_target, y_data):273    y_data = np.array(y_data)274    index = np.abs(y_data - y_target).argmin()275    return index276 277 278def plt_attention(input_x1, attention, fef25, fef50, fef75, fev1_value, fvc_value, is_show=True):279    input_x = input_x1[0][0, :, 0]280    attention_weights = attention[0][0, :, 0]281    length_actual = len(input_x)282    input_x = input_x[:int(length_actual)]283    attention_weights = attention_weights[:int(length_actual)]284    y_data = input_x285    pef_max = np.max(y_data)286    x_pef_max = np.argmax(y_data)287    x_pef25 = (x_pef_max + find_closest_x(fef25, y_data[x_pef_max:])) / 100.0288    x_pef50 = (x_pef_max + find_closest_x(fef50, y_data[x_pef_max:])) / 100.0289    x_pef75 = (x_pef_max + find_closest_x(fef75, y_data[x_pef_max:])) / 100.0290    x_pef_max = x_pef_max / 100.0291 292    fig, ax = plt.subplots(figsize=(12, 8))293    ax.set_xlim(0, 15)294    ax.set_ylim(-0.05, 12)295    ax.spines['bottom'].set_position(('data', -0.05))296    ax.spines['left'].set_position(('data', 0))297    ax.spines['top'].set_visible(False)298    ax.spines['right'].set_visible(False)299    cmap = plt.get_cmap('Reds')300    attention_weights = (attention_weights - attention_weights.min()) / (301            attention_weights.max() - attention_weights.min())302    colors = cmap(attention_weights)303    for j in range(len(input_x) - 1):304        ax.plot([j / 100.0, (j + 1) / 100.0], [input_x[j], input_x[j + 1]], color=colors[j])305    max_y = ax.get_ylim()[1]306    for pef_value, pef_x, label, color in zip([fef25, fef50, fef75, pef_max],307                                              [x_pef25, x_pef50, x_pef75, x_pef_max],308                                              ['FEF25', 'FEF50', 'FEF75', 'PEF'],309                                              ['red', 'blue', 'green', 'purple']):310        ax.vlines(x=pef_x, ymin=0, ymax=pef_value, colors=color, linestyles='--', label=label)311        ax.text(pef_x + 0.1, pef_value, f'{label}', verticalalignment='bottom', horizontalalignment='left',312                color=color)313    ax.annotate(f'FEV1: {fev1_value:.2f}L', xy=(fev1_value + 0.15, 0.0),314                xytext=(fev1_value, -max_y * 0.15),315                textcoords='data',316                va='top',317                ha='center',318                fontsize=6,319                arrowprops=dict(facecolor='orange', shrink=0.05))320    ax.annotate(f'FVC: {fvc_value:.2f}L', xy=(fvc_value - 0.15, 0.0),321                xytext=(fvc_value, -max_y * 0.15),322                textcoords='data',323                va='top',324                ha='center',325                fontsize=6,326                arrowprops=dict(facecolor='grey', shrink=0.05))327    ax.legend()328    ax.set_xlabel('Volume(L)')329    ax.set_ylabel('Flow(L/s)')330    sm = ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=attention_weights.min(), vmax=attention_weights.max()))331    sm.set_array([])332 333    fig.colorbar(sm, ax=ax, orientation='vertical', fraction=0.046, pad=0.04)334    buf = io.BytesIO()335    if is_show:336        plt.show()337    plt.savefig(buf, format='png')338    buf.seek(0)339    return buf340