11776819502

Committed 11 Nov 2024 10:34AM UTC coverage: 0.797% (+0.005%) from 0.792%

Build # 11776819502

Build Type

push

github

Committed by

Pablo1990

Commit Message

intercalations seem to work?????

Run Details

0 of 995 branches covered (0.0%)

Branch coverage included in aggregate %.

0 of 1 new or added line in 1 file covered. (0.0%)

355 existing lines in 6 files now uncovered.

51 of 5408 relevant lines covered (0.94%)

0.01 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/src/pyVertexModel/analysis/analyse_simulation.py

import os
import pickle

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit

from src.pyVertexModel.algorithm.vertexModel import VertexModel, logger
from src.pyVertexModel.util.utils import load_state, load_variables, save_variables


def analyse_simulation(folder):
    """
    Analyse the simulation results
    :param folder:
    :return:
    """

    # Check if the pkl file exists
    if not os.path.exists(os.path.join(folder, 'features_per_time.pkl')):
        # return None, None, None, None
        vModel = VertexModel(create_output_folder=False)

        features_per_time = []
        features_per_time_all_cells = []

        # Go through all the files in the folder
        all_files = os.listdir(folder)

        # Sort files by date
        all_files = sorted(all_files, key=lambda x: os.path.getmtime(os.path.join(folder, x)))

        # if all_files has less than 20 files, return None
        if len(all_files) < 20:
            return None, None, None, None

        for file_id, file in enumerate(all_files):
            if file.endswith('.pkl') and not file.__contains__('data_step_before_remodelling') and not file.__contains__('recoil'):
                # Load the state of the model
                load_state(vModel, os.path.join(folder, file))

                # Analyse the simulation
                all_cells, avg_cells = vModel.analyse_vertex_model()
                features_per_time_all_cells.append(all_cells)
                features_per_time.append(avg_cells)

                # # Create a temporary directory to store the images
                # temp_dir = os.path.join(folder, 'images')
                # if not os.path.exists(temp_dir):
                #     os.mkdir(temp_dir)
                # vModel.screenshot(temp_dir)

                # temp_dir = os.path.join(folder, 'images_wound_edge')
                # if not os.path.exists(temp_dir):
                #     os.mkdir(temp_dir)
                # _, debris_cells = vModel.geo.compute_wound_centre()
                # list_of_cell_distances_top = vModel.geo.compute_cell_distance_to_wound(debris_cells, location_filter=0)
                # alive_cells = [cell.ID for cell in vModel.geo.Cells if cell.AliveStatus == 1]
                # wound_edge_cells = []
                # for cell_num, cell_id in enumerate(alive_cells):
                #     if list_of_cell_distances_top[cell_num] == 1:
                #         wound_edge_cells.append(cell_id)
                # vModel.screenshot(temp_dir, wound_edge_cells)

        if not features_per_time:
            return None, None, None, None

        # Export to xlsx
        features_per_time_all_cells_df = pd.DataFrame(np.concatenate(features_per_time_all_cells),
                                                      columns=features_per_time_all_cells[0].columns)
        features_per_time_all_cells_df.sort_values(by='time', inplace=True)
        features_per_time_all_cells_df.to_excel(os.path.join(folder, 'features_per_time_all_cells.xlsx'))

        features_per_time_df = pd.DataFrame(features_per_time)
        features_per_time_df.sort_values(by='time', inplace=True)
        features_per_time_df.to_excel(os.path.join(folder, 'features_per_time.xlsx'))

        # Obtain pre-wound features
        try:
            pre_wound_features = features_per_time_df['time'][features_per_time_df['time'] < vModel.set.TInitAblation]
            pre_wound_features = features_per_time_df[features_per_time_df['time'] ==
                                                      pre_wound_features.iloc[-1]]
        except Exception as e:
            pre_wound_features = features_per_time_df.iloc[0]

        # Obtain post-wound features
        post_wound_features = features_per_time_df[features_per_time_df['time'] >= vModel.set.TInitAblation]

        if not post_wound_features.empty:
            # Reset time to ablation time.
            post_wound_features.loc[:, 'time'] = post_wound_features['time'] - vModel.set.TInitAblation

            # Compare post-wound features with pre-wound features in percentage
            for feature in post_wound_features.columns:
                if np.any(np.isnan(pre_wound_features[feature])) or np.any(np.isnan(post_wound_features[feature])):
                    continue

                if feature == 'time':
                    continue
                post_wound_features.loc[:, feature] = (post_wound_features[feature] / np.array(
                    pre_wound_features[feature])) * 100

            # Export to xlsx
            post_wound_features.to_excel(os.path.join(folder, 'post_wound_features.xlsx'))

            important_features = calculate_important_features(post_wound_features)
        else:
            important_features = {
                'max_recoiling_top': np.nan,
                'max_recoiling_time_top': np.nan,
                'min_height_change': np.nan,
                'min_height_change_time': np.nan,
            }

        # Export to xlsx
        df = pd.DataFrame([important_features])
        df.to_excel(os.path.join(folder, 'important_features.xlsx'))

        # Save dataframes to a single pkl
        with open(os.path.join(folder, 'features_per_time.pkl'), 'wb') as f:
            pickle.dump(features_per_time_df, f)
            pickle.dump(important_features, f)
            pickle.dump(post_wound_features, f)
            pickle.dump(features_per_time_all_cells_df, f)

    else:
        # Load dataframes from pkl
        with open(os.path.join(folder, 'features_per_time.pkl'), 'rb') as f:
            features_per_time_df = pickle.load(f)
            important_features = pickle.load(f)
            post_wound_features = pickle.load(f)
            features_per_time_all_cells_df = pickle.load(f)

        important_features = calculate_important_features(post_wound_features)

    # Plot wound area top evolution over time and save it to a file
    plot_feature(folder, post_wound_features, name='wound_area_top')
    plot_feature(folder, post_wound_features, name='wound_height')
    plot_feature(folder, post_wound_features, name='num_cells_wound_edge_top')

    return features_per_time_df, post_wound_features, important_features, features_per_time_all_cells_df


def plot_feature(folder, post_wound_features, name='wound_area_top'):
    """
    Plot a feature and save it to a file
    :param folder:
    :param post_wound_features:
    :param name:
    :return:
    """
    plt.figure()
    plt.plot(post_wound_features['time'], post_wound_features[name])
    plt.xlabel('Time (h)')
    plt.ylabel(name)
    # Change axis limits
    if np.max(post_wound_features['time']) > 60:
        plt.xlim([0, np.max(post_wound_features['time'])])
    else:
        plt.xlim([0, 60])
    plt.ylim([0, 250])
    plt.savefig(os.path.join(folder, name + '.png'))
    plt.close()


def calculate_important_features(post_wound_features):
    """
    Calculate important features from the post-wound features
    :param post_wound_features:
    :return:
    """
    # Obtain important features for post-wound
    if not post_wound_features['wound_area_top'].empty and post_wound_features['time'].iloc[-1] > 4:
        important_features = {
            'max_recoiling_top': np.max(post_wound_features['wound_area_top']),
            'max_recoiling_time_top': np.array(post_wound_features['time'])[
                np.argmax(post_wound_features['wound_area_top'])],
            'min_recoiling_top': np.min(post_wound_features['wound_area_top']),
            'min_recoiling_time_top': np.array(post_wound_features['time'])[
                np.argmin(post_wound_features['wound_area_top'])],
            'min_height_change': np.min(post_wound_features['wound_height']),
            'min_height_change_time': np.array(post_wound_features['time'])[
                np.argmin(post_wound_features['wound_height'])],
            'last_area_top': post_wound_features['wound_area_top'].iloc[-1],
            'last_area_time_top': post_wound_features['time'].iloc[-1],
        }

        # Extrapolate features to a given time
        times_to_extrapolate = {3.0, 6.0, 9.0, 12.0, 15.0, 21.0, 30.0, 36.0, 45.0, 51.0, 60.0}
        columns_to_extrapolate = {'wound_area_top', 'wound_height'}  # post_wound_features.columns
        for feature in columns_to_extrapolate:
            for time in times_to_extrapolate:
                # Extrapolate results to a given time
                important_features[feature + '_extrapolated_' + str(time)] = np.interp(time,
                                                                                       post_wound_features['time'],
                                                                                       post_wound_features[feature])

        # # Get ratio from area the first time to the other times
        # for time in times_to_extrapolate:
        #     if time != 6.0:
        #         important_features['ratio_area_top_' + str(time)] = (
        #                 important_features['wound_area_top_extrapolated_' + str(time)] /
        #                 important_features['wound_area_top_extrapolated_6.0'])

    else:
        important_features = {
            'max_recoiling_top': np.nan,
            'max_recoiling_time_top': np.nan,
            'min_height_change': np.nan,
            'min_height_change_time': np.nan,
        }

    return important_features


def analyse_edge_recoil(file_name_v_model, type_of_ablation='recoil_edge_info_apical', n_ablations=2, location_filter=0, t_end=0.5):
    """
    Analyse how much an edge recoil if we ablate an edge of a cell
    :param type_of_ablation:
    :param t_end: Time to iterate after the ablation
    :param file_name_v_model: file nae of the Vertex model
    :param n_ablations: Number of ablations to perform
    :param location_filter: Location filter
    :return:
    """

    v_model = VertexModel(create_output_folder=False)
    load_state(v_model, file_name_v_model)

    # Cells to ablate
    # cell_to_ablate = np.random.choice(possible_cells_to_ablate, 1)
    cell_to_ablate = [v_model.geo.Cells[0]]

    #Pick the neighbouring cell to ablate
    neighbours = cell_to_ablate[0].compute_neighbours(location_filter)

    # Random order of neighbours
    np.random.seed(0)
    np.random.shuffle(neighbours)

    list_of_dicts_to_save = []
    for num_ablation in range(n_ablations):
        load_state(v_model, file_name_v_model)
        try:
            vars = load_variables(file_name_v_model.replace('before_ablation.pkl', type_of_ablation + '.pkl'))
            list_of_dicts_to_save_loaded = vars['recoiling_info_df_apical']

            cell_to_ablate_ID = list_of_dicts_to_save_loaded['cell_to_ablate'][num_ablation]
            neighbour_to_ablate_ID = list_of_dicts_to_save_loaded['neighbour_to_ablate'][num_ablation]
            edge_length_init = list_of_dicts_to_save_loaded['edge_length_init'][num_ablation]
            edge_length_final = list_of_dicts_to_save_loaded['edge_length_final'][num_ablation]
            if 'edge_length_final_normalized' in list_of_dicts_to_save_loaded:
                edge_length_final_normalized = list_of_dicts_to_save_loaded['edge_length_final_normalized'][
                    num_ablation]
            else:
                edge_length_final_normalized = (edge_length_final - edge_length_init) / edge_length_init

            initial_recoil = list_of_dicts_to_save_loaded['initial_recoil_in_s'][num_ablation]
            K = list_of_dicts_to_save_loaded['K'][num_ablation]
            scutoid_face = list_of_dicts_to_save_loaded['scutoid_face'][num_ablation]
            distance_to_centre = list_of_dicts_to_save_loaded['distance_to_centre'][num_ablation]
            if 'time_steps' in list_of_dicts_to_save_loaded:
                time_steps = list_of_dicts_to_save_loaded['time_steps'][num_ablation]
            else:
                time_steps = np.arange(0, len(edge_length_final)) * 6

            if edge_length_final[0] == 0:
                # Remove the first element
                edge_length_final = edge_length_final[1:]
                time_steps = time_steps[1:]
        except Exception as e:
            logger.info('Performing the analysis...' + str(e))
            # Change name of folder and create it
            if type_of_ablation == 'recoil_info_apical':
                v_model.set.OutputFolder = v_model.set.OutputFolder + '_ablation_' + str(num_ablation)
            else:
                v_model.set.OutputFolder = v_model.set.OutputFolder + '_ablation_edge_' + str(num_ablation)

            if not os.path.exists(v_model.set.OutputFolder):
                os.mkdir(v_model.set.OutputFolder)

            neighbour_to_ablate = [neighbours[num_ablation]]

            # Calculate if the cell is neighbour on both sides
            scutoid_face = None
            neighbours_other_side = []
            if location_filter == 0:
                neighbours_other_side = cell_to_ablate[0].compute_neighbours(location_filter=2)
                scutoid_face = np.nan
            elif location_filter == 2:
                neighbours_other_side = cell_to_ablate[0].compute_neighbours(location_filter=0)
                scutoid_face = np.nan

            if scutoid_face is not None:
                if neighbour_to_ablate[0] in neighbours_other_side:
                    scutoid_face = True
                else:
                    scutoid_face = False

            # Get the centre of the tissue
            centre_of_tissue = v_model.geo.compute_centre_of_tissue()
            neighbour_to_ablate_cell = [cell for cell in v_model.geo.Cells if cell.ID == neighbour_to_ablate[0]][0]
            distance_to_centre = np.mean([cell_to_ablate[0].compute_distance_to_centre(centre_of_tissue),
                                          neighbour_to_ablate_cell.compute_distance_to_centre(centre_of_tissue)])

            # Pick the neighbour and put it in the list
            cells_to_ablate = [cell_to_ablate[0].ID, neighbour_to_ablate[0]]

            # Get the edge that share both cells
            edge_length_init = v_model.geo.get_edge_length(cells_to_ablate, location_filter)

            # Ablate the edge
            v_model.set.ablation = True
            v_model.geo.cellsToAblate = cells_to_ablate
            v_model.set.TInitAblation = v_model.t
            if type_of_ablation == 'recoil_info_apical':
                v_model.geo.ablate_cells(v_model.set, v_model.t, combine_cells=False)
                v_model.geo.y_ablated = []
            elif type_of_ablation == 'recoil_edge_info_apical':
                v_model.geo.y_ablated = v_model.geo.ablate_edge(v_model.set, v_model.t, domain=location_filter,
                                                                adjacent_surface=False)

            # Relax the system
            initial_time = v_model.t
            v_model.set.tend = v_model.t + t_end
            if type_of_ablation == 'recoil_info_apical':
                v_model.set.dt = 0.005
            elif type_of_ablation == 'recoil_edge_info_apical':
                v_model.set.dt = 0.005

            v_model.set.Remodelling = False

            v_model.set.dt0 = v_model.set.dt
            if type_of_ablation == 'recoil_edge_info_apical':
                v_model.set.RemodelingFrequency = 0.05
            else:
                v_model.set.RemodelingFrequency = 100
            v_model.set.ablation = False
            v_model.set.export_images = True
            if v_model.set.export_images and not os.path.exists(v_model.set.OutputFolder + '/images'):
                os.mkdir(v_model.set.OutputFolder + '/images')
            edge_length_final_normalized = []
            edge_length_final = []
            recoil_speed = []
            time_steps = []

            # if os.path.exists(v_model.set.OutputFolder):
            #     list_of_files = os.listdir(v_model.set.OutputFolder)
            #     # Get file modification times and sort files by date
            #     files_with_dates = [(file, os.path.getmtime(os.path.join(v_model.set.OutputFolder, file))) for file in
            #                         list_of_files]
            #     files_with_dates.sort(key=lambda x: x[1])
            #     for file in files_with_dates:
            #         load_state(v_model, os.path.join(v_model.set.OutputFolder, file[0]))
            #         compute_edge_length_v_model(cells_to_ablate, edge_length_final, edge_length_final_normalized,
            #                                     edge_length_init, initial_time, location_filter, recoil_speed,
            #                                     time_steps,
            #                                     v_model)

            while v_model.t <= v_model.set.tend and not v_model.didNotConverge:
                gr = v_model.single_iteration()

                compute_edge_length_v_model(cells_to_ablate, edge_length_final, edge_length_final_normalized,
                                            edge_length_init, initial_time, location_filter, recoil_speed, time_steps,
                                            v_model)

                if np.isnan(gr):
                    break

            cell_to_ablate_ID = cell_to_ablate[0].ID
            neighbour_to_ablate_ID = neighbour_to_ablate[0]

        K, initial_recoil, error_bars = fit_ablation_equation(edge_length_final, time_steps)

        # Generate a plot with the edge length final and the fit for each ablation
        plt.figure()
        plt.plot(time_steps, edge_length_final_normalized, 'o')
        # Plot fit line of the Kelvin-Voigt model
        plt.plot(time_steps, recoil_model(np.array(time_steps), initial_recoil, K), 'r')
        plt.xlabel('Time (s)')
        plt.ylabel('Edge length final')
        plt.title('Ablation fit - ' + str(cell_to_ablate_ID) + ' ' + str(neighbour_to_ablate_ID))

        # Save plot
        if type_of_ablation == 'recoil_info_apical':
            plt.savefig(
                os.path.join(file_name_v_model.replace('before_ablation.pkl', 'ablation_fit_' + str(num_ablation) + '.png'))
            )
        elif type_of_ablation == 'recoil_edge_info_apical':
            plt.savefig(
                os.path.join(file_name_v_model.replace('before_ablation.pkl', 'ablation_edge_fit_' + str(num_ablation) + '.png'))
            )
        plt.close()

        # Save the results
        dict_to_save = {
            'cell_to_ablate': cell_to_ablate_ID,
            'neighbour_to_ablate': neighbour_to_ablate_ID,
            'edge_length_init': edge_length_init,
            'edge_length_final': edge_length_final,
            'edge_length_final_normalized': edge_length_final_normalized,
            'initial_recoil_in_s': initial_recoil,
            'K': K,
            'scutoid_face': scutoid_face,
            'location_filter': location_filter,
            'distance_to_centre': distance_to_centre,
            'time_steps': time_steps,
        }
        list_of_dicts_to_save.append(dict_to_save)

    recoiling_info_df_apical = pd.DataFrame(list_of_dicts_to_save)
    recoiling_info_df_apical.to_excel(file_name_v_model.replace('before_ablation.pkl', type_of_ablation+'.xlsx'))
    save_variables({'recoiling_info_df_apical': recoiling_info_df_apical},
                   file_name_v_model.replace('before_ablation.pkl', type_of_ablation+'.pkl'))

    return list_of_dicts_to_save


def recoil_model(x, initial_recoil, K):
    """
    Model of the recoil based on a Kelvin-Voigt model
    :param x:
    :param initial_recoil:
    :param K:
    :return:   Recoil
    """
    return (initial_recoil / K) * (1 - np.exp(-K * x))


def fit_ablation_equation(edge_length_final_normalized, time_steps):
    """
    Fit the ablation equation. Thanks to Veronika Lachina.
    :param edge_length_final_normalized:
    :param time_steps:
    :return:    K, initial_recoil
    """

    # Normalize the edge length
    edge_length_init = edge_length_final_normalized[0]
    edge_length_final_normalized = (edge_length_final_normalized - edge_length_init) / edge_length_init

    # Fit the model to the data
    [params, covariance] = curve_fit(recoil_model, time_steps, edge_length_final_normalized,
                                     p0=[0.00001, 3], bounds=(0, np.inf))

    # Get the error
    error_bars = np.sqrt(np.diag(covariance))

    initial_recoil, K = params
    return K, initial_recoil, error_bars


def compute_edge_length_v_model(cells_to_ablate, edge_length_final, edge_length_final_normalized, edge_length_init,
                                initial_time, location_filter, recoil_speed, time_steps, v_model):
    """
    Compute the edge length of the edge that share the cells_to_ablate
    :param cells_to_ablate:
    :param edge_length_final:
    :param edge_length_final_normalized:
    :param edge_length_init:
    :param initial_time:
    :param location_filter:
    :param recoil_speed:
    :param time_steps:
    :param v_model:
    :return:
    """
    if v_model.t == initial_time:
        return
    # Get the edge length
    edge_length_final.append(v_model.geo.get_edge_length(cells_to_ablate, location_filter))
    edge_length_final_normalized.append((edge_length_final[-1] - edge_length_init) / edge_length_init)
    print('Edge length final: ', edge_length_final[-1])
    # In seconds. 1 t = 1 minute = 60 seconds
    time_steps.append((v_model.t - initial_time) * 60)
    # Calculate the recoil
    recoil_speed.append(edge_length_final_normalized[-1] / time_steps[-1])

1	import os	×
2	import pickle	×
3
4	import numpy as np	×
5	import pandas as pd	×
6	from matplotlib import pyplot as plt	×
7	from scipy.optimize import curve_fit	×
8
9	from src.pyVertexModel.algorithm.vertexModel import VertexModel, logger	×
10	from src.pyVertexModel.util.utils import load_state, load_variables, save_variables	×
11
12
13	def analyse_simulation(folder):	×
14	"""
15	Analyse the simulation results
16	:param folder:
17	:return:
18	"""
19
20	# Check if the pkl file exists
21	if not os.path.exists(os.path.join(folder, 'features_per_time.pkl')):	×
22	# return None, None, None, None
23	vModel = VertexModel(create_output_folder=False)	×
24
25	features_per_time = []	×
26	features_per_time_all_cells = []	×
27
28	# Go through all the files in the folder
29	all_files = os.listdir(folder)	×
30
31	# Sort files by date
32	all_files = sorted(all_files, key=lambda x: os.path.getmtime(os.path.join(folder, x)))	×
33
34	# if all_files has less than 20 files, return None
UNCOV 35	if len(all_files) < 20:	×
36	return None, None, None, None	×
37
UNCOV 38	for file_id, file in enumerate(all_files):	×
39	if file.endswith('.pkl') and not file.__contains__('data_step_before_remodelling') and not file.__contains__('recoil'):	×
40	# Load the state of the model
41	load_state(vModel, os.path.join(folder, file))	×
42
43	# Analyse the simulation
UNCOV 44	all_cells, avg_cells = vModel.analyse_vertex_model()	×
UNCOV 45	features_per_time_all_cells.append(all_cells)	×
UNCOV 46	features_per_time.append(avg_cells)	×
47
48	# # Create a temporary directory to store the images
49	# temp_dir = os.path.join(folder, 'images')
50	# if not os.path.exists(temp_dir):
51	# os.mkdir(temp_dir)
52	# vModel.screenshot(temp_dir)
53
54	# temp_dir = os.path.join(folder, 'images_wound_edge')
55	# if not os.path.exists(temp_dir):
56	# os.mkdir(temp_dir)
57	# _, debris_cells = vModel.geo.compute_wound_centre()
58	# list_of_cell_distances_top = vModel.geo.compute_cell_distance_to_wound(debris_cells, location_filter=0)
59	# alive_cells = [cell.ID for cell in vModel.geo.Cells if cell.AliveStatus == 1]
60	# wound_edge_cells = []
61	# for cell_num, cell_id in enumerate(alive_cells):
62	# if list_of_cell_distances_top[cell_num] == 1:
63	# wound_edge_cells.append(cell_id)
64	# vModel.screenshot(temp_dir, wound_edge_cells)
65
UNCOV 66	if not features_per_time:	×
67	return None, None, None, None	×
68
69	# Export to xlsx
70	features_per_time_all_cells_df = pd.DataFrame(np.concatenate(features_per_time_all_cells),	×
71	columns=features_per_time_all_cells[0].columns)
72	features_per_time_all_cells_df.sort_values(by='time', inplace=True)	×
UNCOV 73	features_per_time_all_cells_df.to_excel(os.path.join(folder, 'features_per_time_all_cells.xlsx'))	×
74
75	features_per_time_df = pd.DataFrame(features_per_time)	×
76	features_per_time_df.sort_values(by='time', inplace=True)	×
77	features_per_time_df.to_excel(os.path.join(folder, 'features_per_time.xlsx'))	×
78
79	# Obtain pre-wound features
80	try:	×
UNCOV 81	pre_wound_features = features_per_time_df['time'][features_per_time_df['time'] < vModel.set.TInitAblation]	×
UNCOV 82	pre_wound_features = features_per_time_df[features_per_time_df['time'] ==	×
83	pre_wound_features.iloc[-1]]
UNCOV 84	except Exception as e:	×
85	pre_wound_features = features_per_time_df.iloc[0]	×
86
87	# Obtain post-wound features
UNCOV 88	post_wound_features = features_per_time_df[features_per_time_df['time'] >= vModel.set.TInitAblation]	×
89
90	if not post_wound_features.empty:	×
91	# Reset time to ablation time.
92	post_wound_features.loc[:, 'time'] = post_wound_features['time'] - vModel.set.TInitAblation	×
93
94	# Compare post-wound features with pre-wound features in percentage
95	for feature in post_wound_features.columns:	×
96	if np.any(np.isnan(pre_wound_features[feature])) or np.any(np.isnan(post_wound_features[feature])):	×
UNCOV 97	continue	×
98
UNCOV 99	if feature == 'time':	×
100	continue	×
UNCOV 101	post_wound_features.loc[:, feature] = (post_wound_features[feature] / np.array(	×
102	pre_wound_features[feature])) * 100
103
104	# Export to xlsx
UNCOV 105	post_wound_features.to_excel(os.path.join(folder, 'post_wound_features.xlsx'))	×
106
UNCOV 107	important_features = calculate_important_features(post_wound_features)	×
108	else:
UNCOV 109	important_features = {	×
110	'max_recoiling_top': np.nan,
111	'max_recoiling_time_top': np.nan,
112	'min_height_change': np.nan,
113	'min_height_change_time': np.nan,
114	}
115
116	# Export to xlsx
117	df = pd.DataFrame([important_features])	×
118	df.to_excel(os.path.join(folder, 'important_features.xlsx'))	×
119
120	# Save dataframes to a single pkl
UNCOV 121	with open(os.path.join(folder, 'features_per_time.pkl'), 'wb') as f:	×
UNCOV 122	pickle.dump(features_per_time_df, f)	×
UNCOV 123	pickle.dump(important_features, f)	×
124	pickle.dump(post_wound_features, f)	×
125	pickle.dump(features_per_time_all_cells_df, f)	×
126
127	else:
128	# Load dataframes from pkl
UNCOV 129	with open(os.path.join(folder, 'features_per_time.pkl'), 'rb') as f:	×
130	features_per_time_df = pickle.load(f)	×
UNCOV 131	important_features = pickle.load(f)	×
UNCOV 132	post_wound_features = pickle.load(f)	×
133	features_per_time_all_cells_df = pickle.load(f)	×
134
135	important_features = calculate_important_features(post_wound_features)	×
136
137	# Plot wound area top evolution over time and save it to a file
UNCOV 138	plot_feature(folder, post_wound_features, name='wound_area_top')	×
UNCOV 139	plot_feature(folder, post_wound_features, name='wound_height')	×
140	plot_feature(folder, post_wound_features, name='num_cells_wound_edge_top')	×
141
UNCOV 142	return features_per_time_df, post_wound_features, important_features, features_per_time_all_cells_df	×
143
144
UNCOV 145	def plot_feature(folder, post_wound_features, name='wound_area_top'):	×
146	"""
147	Plot a feature and save it to a file
148	:param folder:
149	:param post_wound_features:
150	:param name:
151	:return:
152	"""
153	plt.figure()	×
154	plt.plot(post_wound_features['time'], post_wound_features[name])	×
155	plt.xlabel('Time (h)')	×
156	plt.ylabel(name)	×
157	# Change axis limits
UNCOV 158	if np.max(post_wound_features['time']) > 60:	×
159	plt.xlim([0, np.max(post_wound_features['time'])])	×
160	else:
UNCOV 161	plt.xlim([0, 60])	×
UNCOV 162	plt.ylim([0, 250])	×
UNCOV 163	plt.savefig(os.path.join(folder, name + '.png'))	×
UNCOV 164	plt.close()	×
165
166
167	def calculate_important_features(post_wound_features):	×
168	"""
169	Calculate important features from the post-wound features
170	:param post_wound_features:
171	:return:
172	"""
173	# Obtain important features for post-wound
UNCOV 174	if not post_wound_features['wound_area_top'].empty and post_wound_features['time'].iloc[-1] > 4:	×
UNCOV 175	important_features = {	×
176	'max_recoiling_top': np.max(post_wound_features['wound_area_top']),
177	'max_recoiling_time_top': np.array(post_wound_features['time'])[
178	np.argmax(post_wound_features['wound_area_top'])],
179	'min_recoiling_top': np.min(post_wound_features['wound_area_top']),
180	'min_recoiling_time_top': np.array(post_wound_features['time'])[
181	np.argmin(post_wound_features['wound_area_top'])],
182	'min_height_change': np.min(post_wound_features['wound_height']),
183	'min_height_change_time': np.array(post_wound_features['time'])[
184	np.argmin(post_wound_features['wound_height'])],
185	'last_area_top': post_wound_features['wound_area_top'].iloc[-1],
186	'last_area_time_top': post_wound_features['time'].iloc[-1],
187	}
188
189	# Extrapolate features to a given time
UNCOV 190	times_to_extrapolate = {3.0, 6.0, 9.0, 12.0, 15.0, 21.0, 30.0, 36.0, 45.0, 51.0, 60.0}	×
UNCOV 191	columns_to_extrapolate = {'wound_area_top', 'wound_height'} # post_wound_features.columns	×
UNCOV 192	for feature in columns_to_extrapolate:	×
UNCOV 193	for time in times_to_extrapolate:	×
194	# Extrapolate results to a given time
UNCOV 195	important_features[feature + '_extrapolated_' + str(time)] = np.interp(time,	×
196	post_wound_features['time'],
197	post_wound_features[feature])
198
199	# # Get ratio from area the first time to the other times
200	# for time in times_to_extrapolate:
201	# if time != 6.0:
202	# important_features['ratio_area_top_' + str(time)] = (
203	# important_features['wound_area_top_extrapolated_' + str(time)] /
204	# important_features['wound_area_top_extrapolated_6.0'])
205
206	else:
UNCOV 207	important_features = {	×
208	'max_recoiling_top': np.nan,
209	'max_recoiling_time_top': np.nan,
210	'min_height_change': np.nan,
211	'min_height_change_time': np.nan,
212	}
213
UNCOV 214	return important_features	×
215
216
UNCOV 217	def analyse_edge_recoil(file_name_v_model, type_of_ablation='recoil_edge_info_apical', n_ablations=2, location_filter=0, t_end=0.5):	×
218	"""
219	Analyse how much an edge recoil if we ablate an edge of a cell
220	:param type_of_ablation:
221	:param t_end: Time to iterate after the ablation
222	:param file_name_v_model: file nae of the Vertex model
223	:param n_ablations: Number of ablations to perform
224	:param location_filter: Location filter
225	:return:
226	"""
227
228	v_model = VertexModel(create_output_folder=False)	×
UNCOV 229	load_state(v_model, file_name_v_model)	×
230
231	# Cells to ablate
232	# cell_to_ablate = np.random.choice(possible_cells_to_ablate, 1)
UNCOV 233	cell_to_ablate = [v_model.geo.Cells[0]]	×
234
235	#Pick the neighbouring cell to ablate
236	neighbours = cell_to_ablate[0].compute_neighbours(location_filter)	×
237
238	# Random order of neighbours
239	np.random.seed(0)	×
UNCOV 240	np.random.shuffle(neighbours)	×
241
242	list_of_dicts_to_save = []	×
243	for num_ablation in range(n_ablations):	×
244	load_state(v_model, file_name_v_model)	×
245	try:	×
246	vars = load_variables(file_name_v_model.replace('before_ablation.pkl', type_of_ablation + '.pkl'))	×
UNCOV 247	list_of_dicts_to_save_loaded = vars['recoiling_info_df_apical']	×
248
249	cell_to_ablate_ID = list_of_dicts_to_save_loaded['cell_to_ablate'][num_ablation]	×
UNCOV 250	neighbour_to_ablate_ID = list_of_dicts_to_save_loaded['neighbour_to_ablate'][num_ablation]	×
251	edge_length_init = list_of_dicts_to_save_loaded['edge_length_init'][num_ablation]	×
252	edge_length_final = list_of_dicts_to_save_loaded['edge_length_final'][num_ablation]	×
253	if 'edge_length_final_normalized' in list_of_dicts_to_save_loaded:	×
254	edge_length_final_normalized = list_of_dicts_to_save_loaded['edge_length_final_normalized'][	×
255	num_ablation]
256	else:
UNCOV 257	edge_length_final_normalized = (edge_length_final - edge_length_init) / edge_length_init	×
258
UNCOV 259	initial_recoil = list_of_dicts_to_save_loaded['initial_recoil_in_s'][num_ablation]	×
260	K = list_of_dicts_to_save_loaded['K'][num_ablation]	×
UNCOV 261	scutoid_face = list_of_dicts_to_save_loaded['scutoid_face'][num_ablation]	×
262	distance_to_centre = list_of_dicts_to_save_loaded['distance_to_centre'][num_ablation]	×
263	if 'time_steps' in list_of_dicts_to_save_loaded:	×
264	time_steps = list_of_dicts_to_save_loaded['time_steps'][num_ablation]	×
265	else:
UNCOV 266	time_steps = np.arange(0, len(edge_length_final)) * 6	×
267
268	if edge_length_final[0] == 0:	×
269	# Remove the first element
270	edge_length_final = edge_length_final[1:]	×
UNCOV 271	time_steps = time_steps[1:]	×
272	except Exception as e:	×
273	logger.info('Performing the analysis...' + str(e))	×
274	# Change name of folder and create it
275	if type_of_ablation == 'recoil_info_apical':	×
UNCOV 276	v_model.set.OutputFolder = v_model.set.OutputFolder + '_ablation_' + str(num_ablation)	×
277	else:
278	v_model.set.OutputFolder = v_model.set.OutputFolder + '_ablation_edge_' + str(num_ablation)	×
279
280	if not os.path.exists(v_model.set.OutputFolder):	×
281	os.mkdir(v_model.set.OutputFolder)	×
282
283	neighbour_to_ablate = [neighbours[num_ablation]]	×
284
285	# Calculate if the cell is neighbour on both sides
UNCOV 286	scutoid_face = None	×
287	neighbours_other_side = []	×
288	if location_filter == 0:	×
289	neighbours_other_side = cell_to_ablate[0].compute_neighbours(location_filter=2)	×
UNCOV 290	scutoid_face = np.nan	×
291	elif location_filter == 2:	×
UNCOV 292	neighbours_other_side = cell_to_ablate[0].compute_neighbours(location_filter=0)	×
UNCOV 293	scutoid_face = np.nan	×
294
295	if scutoid_face is not None:	×
296	if neighbour_to_ablate[0] in neighbours_other_side:	×
UNCOV 297	scutoid_face = True	×
298	else:
UNCOV 299	scutoid_face = False	×
300
301	# Get the centre of the tissue
UNCOV 302	centre_of_tissue = v_model.geo.compute_centre_of_tissue()	×
303	neighbour_to_ablate_cell = [cell for cell in v_model.geo.Cells if cell.ID == neighbour_to_ablate[0]][0]	×
UNCOV 304	distance_to_centre = np.mean([cell_to_ablate[0].compute_distance_to_centre(centre_of_tissue),	×
305	neighbour_to_ablate_cell.compute_distance_to_centre(centre_of_tissue)])
306
307	# Pick the neighbour and put it in the list
308	cells_to_ablate = [cell_to_ablate[0].ID, neighbour_to_ablate[0]]	×
309
310	# Get the edge that share both cells
311	edge_length_init = v_model.geo.get_edge_length(cells_to_ablate, location_filter)	×
312
313	# Ablate the edge
UNCOV 314	v_model.set.ablation = True	×
UNCOV 315	v_model.geo.cellsToAblate = cells_to_ablate	×
UNCOV 316	v_model.set.TInitAblation = v_model.t	×
317	if type_of_ablation == 'recoil_info_apical':	×
318	v_model.geo.ablate_cells(v_model.set, v_model.t, combine_cells=False)	×
319	v_model.geo.y_ablated = []	×
320	elif type_of_ablation == 'recoil_edge_info_apical':	×
321	v_model.geo.y_ablated = v_model.geo.ablate_edge(v_model.set, v_model.t, domain=location_filter,	×
322	adjacent_surface=False)
323
324	# Relax the system
UNCOV 325	initial_time = v_model.t	×
326	v_model.set.tend = v_model.t + t_end	×
327	if type_of_ablation == 'recoil_info_apical':	×
328	v_model.set.dt = 0.005	×
UNCOV 329	elif type_of_ablation == 'recoil_edge_info_apical':	×
330	v_model.set.dt = 0.005	×
331
332	v_model.set.Remodelling = False	×
333
334	v_model.set.dt0 = v_model.set.dt	×
335	if type_of_ablation == 'recoil_edge_info_apical':	×
336	v_model.set.RemodelingFrequency = 0.05	×
337	else:
338	v_model.set.RemodelingFrequency = 100	×
UNCOV 339	v_model.set.ablation = False	×
UNCOV 340	v_model.set.export_images = True	×
UNCOV 341	if v_model.set.export_images and not os.path.exists(v_model.set.OutputFolder + '/images'):	×
UNCOV 342	os.mkdir(v_model.set.OutputFolder + '/images')	×
UNCOV 343	edge_length_final_normalized = []	×
UNCOV 344	edge_length_final = []	×
UNCOV 345	recoil_speed = []	×
UNCOV 346	time_steps = []	×
347
348	# if os.path.exists(v_model.set.OutputFolder):
349	# list_of_files = os.listdir(v_model.set.OutputFolder)
350	# # Get file modification times and sort files by date
351	# files_with_dates = [(file, os.path.getmtime(os.path.join(v_model.set.OutputFolder, file))) for file in
352	# list_of_files]
353	# files_with_dates.sort(key=lambda x: x[1])
354	# for file in files_with_dates:
355	# load_state(v_model, os.path.join(v_model.set.OutputFolder, file[0]))
356	# compute_edge_length_v_model(cells_to_ablate, edge_length_final, edge_length_final_normalized,
357	# edge_length_init, initial_time, location_filter, recoil_speed,
358	# time_steps,
359	# v_model)
360
361	while v_model.t <= v_model.set.tend and not v_model.didNotConverge:	×
UNCOV 362	gr = v_model.single_iteration()	×
363
364	compute_edge_length_v_model(cells_to_ablate, edge_length_final, edge_length_final_normalized,	×
365	edge_length_init, initial_time, location_filter, recoil_speed, time_steps,
366	v_model)
367
UNCOV 368	if np.isnan(gr):	×
369	break	×
370
UNCOV 371	cell_to_ablate_ID = cell_to_ablate[0].ID	×
372	neighbour_to_ablate_ID = neighbour_to_ablate[0]	×
373
374	K, initial_recoil, error_bars = fit_ablation_equation(edge_length_final, time_steps)	×
375
376	# Generate a plot with the edge length final and the fit for each ablation
UNCOV 377	plt.figure()	×
378	plt.plot(time_steps, edge_length_final_normalized, 'o')	×
379	# Plot fit line of the Kelvin-Voigt model
UNCOV 380	plt.plot(time_steps, recoil_model(np.array(time_steps), initial_recoil, K), 'r')	×
UNCOV 381	plt.xlabel('Time (s)')	×
382	plt.ylabel('Edge length final')	×
383	plt.title('Ablation fit - ' + str(cell_to_ablate_ID) + ' ' + str(neighbour_to_ablate_ID))	×
384
385	# Save plot
386	if type_of_ablation == 'recoil_info_apical':	×
UNCOV 387	plt.savefig(	×
388	os.path.join(file_name_v_model.replace('before_ablation.pkl', 'ablation_fit_' + str(num_ablation) + '.png'))
389	)
UNCOV 390	elif type_of_ablation == 'recoil_edge_info_apical':	×
UNCOV 391	plt.savefig(	×
392	os.path.join(file_name_v_model.replace('before_ablation.pkl', 'ablation_edge_fit_' + str(num_ablation) + '.png'))
393	)
UNCOV 394	plt.close()	×
395
396	# Save the results
UNCOV 397	dict_to_save = {	×
398	'cell_to_ablate': cell_to_ablate_ID,
399	'neighbour_to_ablate': neighbour_to_ablate_ID,
400	'edge_length_init': edge_length_init,
401	'edge_length_final': edge_length_final,
402	'edge_length_final_normalized': edge_length_final_normalized,
403	'initial_recoil_in_s': initial_recoil,
404	'K': K,
405	'scutoid_face': scutoid_face,
406	'location_filter': location_filter,
407	'distance_to_centre': distance_to_centre,
408	'time_steps': time_steps,
409	}
UNCOV 410	list_of_dicts_to_save.append(dict_to_save)	×
411
412	recoiling_info_df_apical = pd.DataFrame(list_of_dicts_to_save)	×
UNCOV 413	recoiling_info_df_apical.to_excel(file_name_v_model.replace('before_ablation.pkl', type_of_ablation+'.xlsx'))	×
UNCOV 414	save_variables({'recoiling_info_df_apical': recoiling_info_df_apical},	×
415	file_name_v_model.replace('before_ablation.pkl', type_of_ablation+'.pkl'))
416
UNCOV 417	return list_of_dicts_to_save	×
418
419
420	def recoil_model(x, initial_recoil, K):	×
421	"""
422	Model of the recoil based on a Kelvin-Voigt model
423	:param x:
424	:param initial_recoil:
425	:param K:
426	:return: Recoil
427	"""
UNCOV 428	return (initial_recoil / K) * (1 - np.exp(-K * x))	×
429
430
UNCOV 431	def fit_ablation_equation(edge_length_final_normalized, time_steps):	×
432	"""
433	Fit the ablation equation. Thanks to Veronika Lachina.
434	:param edge_length_final_normalized:
435	:param time_steps:
436	:return: K, initial_recoil
437	"""
438
439	# Normalize the edge length
440	edge_length_init = edge_length_final_normalized[0]	×
UNCOV 441	edge_length_final_normalized = (edge_length_final_normalized - edge_length_init) / edge_length_init	×
442
443	# Fit the model to the data
UNCOV 444	[params, covariance] = curve_fit(recoil_model, time_steps, edge_length_final_normalized,	×
445	p0=[0.00001, 3], bounds=(0, np.inf))
446
447	# Get the error
UNCOV 448	error_bars = np.sqrt(np.diag(covariance))	×
449
UNCOV 450	initial_recoil, K = params	×
UNCOV 451	return K, initial_recoil, error_bars	×
452
453
UNCOV 454	def compute_edge_length_v_model(cells_to_ablate, edge_length_final, edge_length_final_normalized, edge_length_init,	×
455	initial_time, location_filter, recoil_speed, time_steps, v_model):
456	"""
457	Compute the edge length of the edge that share the cells_to_ablate
458	:param cells_to_ablate:
459	:param edge_length_final:
460	:param edge_length_final_normalized:
461	:param edge_length_init:
462	:param initial_time:
463	:param location_filter:
464	:param recoil_speed:
465	:param time_steps:
466	:param v_model:
467	:return:
468	"""
UNCOV 469	if v_model.t == initial_time:	×
470	return	×
471	# Get the edge length
UNCOV 472	edge_length_final.append(v_model.geo.get_edge_length(cells_to_ablate, location_filter))	×
UNCOV 473	edge_length_final_normalized.append((edge_length_final[-1] - edge_length_init) / edge_length_init)	×
UNCOV 474	print('Edge length final: ', edge_length_final[-1])	×
475	# In seconds. 1 t = 1 minute = 60 seconds
UNCOV 476	time_steps.append((v_model.t - initial_time) * 60)	×
477	# Calculate the recoil
UNCOV 478	recoil_speed.append(edge_length_final_normalized[-1] / time_steps[-1])	×

Pablo1990 / pyVertexModel / 11776819502

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous