11294040206

Committed 11 Oct 2024 02:19PM UTC coverage: 8.803% (-0.001%) from 8.804%

Build # 11294040206

Build Type

push

github

Committed by

Pablo1990

Commit Message

bugfix, but needs thinking

Run Details

0 of 1005 branches covered (0.0%)

Branch coverage included in aggregate %.

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

59 existing lines in 1 file now uncovered.

565 of 5413 relevant lines covered (10.44%)

0.1 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)
        all_files.sort()
        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
        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]]

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

Pablo1990 / pyVertexModel / 11294040206

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