diff --git a/.gitignore b/.gitignore index 4f946e3a..9c52c374 100644 --- a/.gitignore +++ b/.gitignore @@ -84,3 +84,9 @@ dist/ .idea/ src/tyssue/_version.py + +# closure-test rendered videos +tests/behaviors/output/ + +# notebook-generated demo media +notebooks/closure_demo_output/ diff --git a/src/tyssue/behaviors/sheet/actions.py b/src/tyssue/behaviors/sheet/actions.py index 8f6c621a..52150ede 100644 --- a/src/tyssue/behaviors/sheet/actions.py +++ b/src/tyssue/behaviors/sheet/actions.py @@ -39,6 +39,9 @@ def merge_vertices(sheet): collapse_edge(sheet, short[0], allow_two_sided=False) short = sheet.edge_df[sheet.edge_df["length"] < d_min].index.to_numpy() np.random.shuffle(short) + + sheet.network_changed = True + return 0 diff --git a/src/tyssue/core/history.py b/src/tyssue/core/history.py index 7a4c1879..3bbe082d 100644 --- a/src/tyssue/core/history.py +++ b/src/tyssue/core/history.py @@ -34,13 +34,13 @@ class History: """ def __init__( - self, - sheet, - save_every=None, - dt=None, - save_only=None, - extra_cols=None, - save_all=True, + self, + sheet, + save_every=None, + dt=None, + save_only=None, + extra_cols=None, + save_all=True ): """Creates a `SheetHistory` instance. @@ -60,14 +60,14 @@ def __init__( """ if extra_cols is not None: warnings.warn( - "extra_cols and save_all parameters are deprecated." - " Use save_only instead. " - ) + "extra_cols and save_all parameters are deprecated. Use save_only instead. ") + + extra_cols = { + k: list(sheet.datasets[k].columns) for k in sheet.datasets + } if save_only is not None: - extra_cols = defaultdict(list, **save_only) - else: - extra_cols = {k: list(sheet.datasets[k].columns) for k in sheet.datasets} + extra_cols = defaultdict(list, **extra_cols) self.sheet = sheet @@ -80,31 +80,35 @@ def __init__( self.save_every = None self.datasets = {} + self.dicts = {} self.columns = {} vcols = sheet.coords + extra_cols["vert"] vcols = list(set(vcols)) self.vcols = _filter_columns(vcols, sheet.vert_df.columns, "vertex") _vert_h = sheet.vert_df[self.vcols].reset_index(drop=False) - if "time" not in self.vcols: + if not "time" in self.vcols: _vert_h["time"] = 0 self.datasets["vert"] = _vert_h + self.dicts["vert"] = {} self.columns["vert"] = self.vcols fcols = extra_cols["face"] self.fcols = _filter_columns(fcols, sheet.face_df.columns, "face") _face_h = sheet.face_df[self.fcols].reset_index(drop=False) - if "time" not in self.fcols: + if not "time" in self.fcols: _face_h["time"] = 0 self.datasets["face"] = _face_h + self.dicts["face"] = {} self.columns["face"] = self.fcols if sheet.cell_df is not None: ccols = extra_cols["cell"] self.ccols = _filter_columns(ccols, sheet.cell_df.columns, "cell") _cell_h = sheet.cell_df[self.ccols].reset_index(drop=False) - if "time" not in self.ccols: + if not "time" in self.ccols: _cell_h["time"] = 0 self.datasets["cell"] = _cell_h + self.dicts["cell"] = {} self.columns["cell"] = self.ccols extra_cols["edge"].append("cell") @@ -112,9 +116,10 @@ def __init__( ecols = list(set(ecols)) self.ecols = _filter_columns(ecols, sheet.edge_df.columns, "edge") _edge_h = sheet.edge_df[self.ecols].reset_index(drop=False) - if "time" not in self.ecols: + if not "time" in self.ecols: _edge_h["time"] = 0 self.datasets["edge"] = _edge_h + self.dicts["edge"] = {} self.columns["edge"] = self.ecols def __len__(self): @@ -168,46 +173,52 @@ def record(self, time_stamp=None): self.time += 1 if (self.save_every is None) or ( - self.index % (int(self.save_every / self.dt)) == 0 + self.index % (int(self.save_every / self.dt)) == 0 ): for element in self.datasets: hist = self.datasets[element] cols = self.columns[element] df = self.sheet.datasets[element][cols].reset_index(drop=False) - if "time" not in cols: - times = pd.Series(np.ones((df.shape[0],)) * self.time, name="time") - df = pd.concat([df, times], ignore_index=False, axis=1, sort=False) - else: - df["time"] = self.time - - if self.time in hist["time"]: - # erase previously recorded time point - hist = hist[hist["time"] != self.time] + # if "time" not in cols: + # times = pd.Series(np.ones((df.shape[0],)) * self.time, name="time") + # df = pd.concat([df, times], ignore_index=False, axis=1, sort=False) + # else: + df["time"] = self.time - hist = pd.concat([hist, df], ignore_index=True, axis=0, sort=False) + # if self.time in hist["time"]: + # # erase previously recorded time point + # hist = hist[hist["time"] != self.time] - self.datasets[element] = hist + self.dicts[element].update({f"{self.time}": df}) self.index += 1 + def update_datasets(self): + """Concatenate all datasets in self.datasets into self.datasets as pd.DataFrame objects + """ + for element in self.sheet.datasets: + self.datasets[element] = pd.concat(self.dicts[element].values(), ignore_index=True) + def retrieve(self, time): """Return datasets at time `time`. If a specific dataset was not recorded at time time, the closest record before that time is used. """ - if time > self.datasets["vert"]["time"].values[-1]: + times = [float(_time) for _time in self.dicts["vert"].keys()] + + if time > max(times): warnings.warn( """ The time argument you requested is bigger than the maximum recorded time, are you sure you passed the time stamp as parameter, and not an index ? """ ) + t = times[np.argmin([np.abs(t1 - time) for t1 in times])] sheet_datasets = {} for element in self.datasets: - hist = self.datasets[element] + df = self.dicts[element][f"{t}"] cols = self.columns[element] - df = _retrieve(hist, time) df = df.set_index(element)[cols] sheet_datasets[element] = df @@ -229,7 +240,7 @@ def slice(self, start=0, stop=None, size=None, endpoint=True): """ if size is not None: if stop is not None: - time_stamps = self.time_stamps[start : stop + int(endpoint)] + time_stamps = self.time_stamps[start: stop + int(endpoint)] else: time_stamps = self.time_stamps indices = np.round( @@ -237,7 +248,7 @@ def slice(self, start=0, stop=None, size=None, endpoint=True): ).astype(int) times = time_stamps.take(indices.clip(max=time_stamps.size - 1)) elif stop is not None: - times = self.time_stamps[start : stop + int(endpoint)] + times = self.time_stamps[start: stop + int(endpoint)] else: times = self.time_stamps return times @@ -270,13 +281,13 @@ class HistoryHdf5(History): """ def __init__( - self, - sheet=None, - save_every=None, - dt=None, - save_only=None, - hf5file="", - overwrite=False, + self, + sheet=None, + save_every=None, + dt=None, + save_only=None, + hf5file="", + overwrite=False, ): """Creates a `HistoryHdf5` instance. @@ -373,11 +384,10 @@ def time_stamps(self, element="vert"): with pd.HDFStore(self.hf5file, "r") as file: self._time_stamps = file.select("vert", columns=["time"])[ "time" - ].unique() + ].unique() return self._time_stamps - def record(self, time_stamp=None, sheet=None): """Appends a copy of the sheet datasets to the history HDF file. @@ -399,7 +409,6 @@ def record(self, time_stamp=None, sheet=None): # invalidate _time_stamp cache: self._time_stamps = np.empty((0,)) - dtypes_ = {k: df.dtypes for k, df in self.sheet.datasets.items()} for element, df in self.sheet.datasets.items(): @@ -418,7 +427,7 @@ def record(self, time_stamp=None, sheet=None): ) if (self.save_every is None) or ( - self.index % (int(self.save_every / self.dt)) == 0 + self.index % (int(self.save_every / self.dt)) == 0 ): for element, df in self.sheet.datasets.items(): times = pd.Series(np.ones((df.shape[0],)) * self.time, name="time") @@ -429,11 +438,11 @@ def record(self, time_stamp=None, sheet=None): kwargs["min_itemsize"] = {"segment": 8} with pd.HDFStore(self.hf5file, "a") as store: if ( - element in store - and store.select(element, where=f"time == {self.time}")[ - "time" - ].shape[0] - > 0 + element in store + and store.select(element, where=f"time == {self.time}")[ + "time" + ].shape[0] + > 0 ): store.remove(key=element, where=f"time == {self.time}") store.append(key=element, value=df, **kwargs) diff --git a/src/tyssue/core/objects.py b/src/tyssue/core/objects.py index ac07db0a..5214128a 100644 --- a/src/tyssue/core/objects.py +++ b/src/tyssue/core/objects.py @@ -663,7 +663,11 @@ def face_polygons(self, coords=None): for c in coords: self.edge_df["s" + c] = self.upcast_srce(self.vert_df[c]) - polys = self.edge_df.groupby("face").apply(lambda df: df[scoords].to_numpy()) + polys = ( + self.edge_df[scoords + ["face"]] + .groupby("face")[scoords] + .apply(lambda df: df.to_numpy()) + ) return polys def validate(self): @@ -677,7 +681,7 @@ def get_valid(self): """Set the 'is_valid' column to true if the faces are all closed polygons, and the cells closed polyhedra. """ - is_valid_face = self.edge_df.groupby("face").apply(_test_valid) + is_valid_face = self.edge_df.groupby("face").apply(_test_valid, include_groups=True) is_valid = self.upcast_face(is_valid_face) if "cell" in self.data_names: is_valid_cell = self.edge_df.groupby("cell").apply(_is_closed_cell) diff --git a/src/tyssue/core/sheet.py b/src/tyssue/core/sheet.py index a7753895..bd578e43 100644 --- a/src/tyssue/core/sheet.py +++ b/src/tyssue/core/sheet.py @@ -468,28 +468,39 @@ def planar_sheet_3d(cls, identifier, nx, ny, distx, disty, noise=None): def get_opposite(edge_df, raise_if_invalid=False): - """ - Returns the indices opposite to the edges in `edge_df` - """ + srce = edge_df["srce"].to_numpy() + trgt = edge_df["trgt"].to_numpy() + edge_idx = edge_df.index.to_numpy() - st_indexed = ( - edge_df[["srce", "trgt"]].reset_index().set_index(["srce", "trgt"], drop=False) - ) - flipped = st_indexed.index.swaplevel(0, 1) - flipped.names = ["srce", "trgt"] - try: - opposite = st_indexed.reindex(flipped)["edge"].values - except ValueError as e: - dup = flipped.duplicated() - warnings.warn( - "Duplicated (`srce`, `trgt`) values in edge_df, maybe sanitize your input" - ) - opposite = st_indexed[~dup].reindex(flipped)["edge"].values + pairs = np.zeros(len(srce), dtype=[('f0', srce.dtype), ('f1', trgt.dtype)]) + pairs['f0'], pairs['f1'] = srce, trgt + + sort_idx = np.argsort(pairs) + sorted_pairs = pairs[sort_idx] + + if np.any(sorted_pairs[1:] == sorted_pairs[:-1]): + warnings.warn("Duplicated (`srce`, `trgt`) values detected.") if raise_if_invalid: - raise e + raise ValueError("Duplicated pairs detected") + + # Find opposites + # Create the "target" pairs we are looking for: (trgt, srce) + flipped_pairs = np.zeros(len(srce), dtype=[('f0', srce.dtype), ('f1', trgt.dtype)]) + flipped_pairs['f0'], flipped_pairs['f1'] = trgt, srce + + # Find where flipped_pairs would fit into the sorted list of original pairs + match_indices = np.searchsorted(sorted_pairs, flipped_pairs) + + + valid_mask = (match_indices < len(pairs)) + actual_match_mask = valid_mask.copy() + actual_match_mask[valid_mask] = (sorted_pairs[match_indices[valid_mask]] == flipped_pairs[valid_mask]) + + # Map back from the sorted index to the original edge index + opposite = np.full(len(srce), -1, dtype=edge_idx.dtype) + opposite[actual_match_mask] = edge_idx[sort_idx[match_indices[actual_match_mask]]] - opposite[np.isnan(opposite)] = -1 - return opposite.astype(int) + return opposite def get_outer_sheet(eptm): diff --git a/src/tyssue/draw/__init__.py b/src/tyssue/draw/__init__.py index 69d31784..fae6a59e 100644 --- a/src/tyssue/draw/__init__.py +++ b/src/tyssue/draw/__init__.py @@ -1,6 +1,10 @@ +from mpl_toolkits.mplot3d.art3d import Line3DCollection, Poly3DCollection +import matplotlib.collections as mcollections +import numpy as np + from .ipv_draw import browse_history # noqa from .ipv_draw import sheet_view as sheet_view_3d # noqa -from .plt_draw import create_gif, plot_forces, quick_edge_draw # noqa +from .plt_draw import create_gif, create_gif_3d, plot_forces, quick_edge_draw # noqa from .plt_draw import sheet_view as sheet_view_2d # noqa try: @@ -9,7 +13,6 @@ print("vispy won't work") sheet_view_vispy = None - def sheet_view(sheet, coords=["x", "y", "z"], ax=None, mode="2D", **draw_specs_kw): """Main plotting function in 2D or 3D. diff --git a/src/tyssue/draw/plt_draw.py b/src/tyssue/draw/plt_draw.py index d48c3f0f..d46829e8 100644 --- a/src/tyssue/draw/plt_draw.py +++ b/src/tyssue/draw/plt_draw.py @@ -17,15 +17,48 @@ from matplotlib.collections import LineCollection, PatchCollection, PolyCollection from matplotlib.patches import Arc, FancyArrow, PathPatch from matplotlib.path import Path +from mpl_toolkits.mplot3d.art3d import Line3DCollection, Poly3DCollection +import matplotlib.collections as mcollections +import matplotlib.patches as mpatches from ..config.draw import sheet_spec from ..utils.utils import get_sub_eptm, spec_updater COORDS = ["x", "y"] +COORDS3D = ["x", "y", "z"] log = logging.getLogger(__name__) +def deep_update(base, updates): + for key, value in updates.items(): + if key in base and isinstance(base[key], dict) and isinstance(value, dict): + deep_update(base[key], value) + else: + base[key] = value + return base + + +def patch_2d_collections_to_3d(ax): + """Replace any 2D LineCollections on a 3D axes with proper Line3DCollection instances.""" + replacements = [] + for col in ax.collections: + if type(col) is mcollections.LineCollection: + segments = col.get_segments() + segments_3d = [ + seg if seg.shape[1] == 3 else np.hstack([seg, np.zeros((len(seg), 1))]) + for seg in segments + ] + new_col = Line3DCollection(segments_3d) + new_col.set_color(col.get_colors()) + new_col.set_linewidth(col.get_linewidths()) + replacements.append((col, new_col)) + + for old, new in replacements: + old.remove() + ax.add_collection3d(new) + + def browse_history( history, coords=["x", "y"], @@ -75,6 +108,7 @@ def create_gif( interval=None, draw_func=None, margin=5, + dpi=200, **draw_kwds, ): """Creates an animated gif of the recorded history. @@ -118,6 +152,7 @@ def create_gif( for i, (t, sheet) in enumerate(history.browse(start, stop, num_frames)): try: fig, ax = draw_func(sheet, **draw_kwds) + plt.title(f"t = {t:.2f}") except Exception as e: print(f"Droped frame {i}") print(e) @@ -125,11 +160,15 @@ def create_gif( if isinstance(ax, plt.Axes) and margin >= 0: ax.set(xlim=xlim, ylim=ylim) - fig.savefig(graph_dir / f"movie_{i:04d}.png") + fig.savefig( + graph_dir / f"movie_{i:04d}.png", + dpi = dpi, + bbox_inches="tight", + ) plt.close(fig) try: - subprocess.run(["convert", (graph_dir / "movie_*.png").as_posix(), output]) + subprocess.run(["magick", (graph_dir / "movie_*.png").as_posix(), output]) except Exception as e: print( "Converting didn't work, make sure imagemagick is available on your system" @@ -139,8 +178,117 @@ def create_gif( finally: shutil.rmtree(graph_dir) +def create_gif_3d( + history, + output, + num_frames=None, + interval=None, + draw_func=None, + margin=5, + dpi=200, + view_angle=(30, 45), + dynamic_draw_kwds=None, + legend = None, + cull_back_edges=False, + **draw_kwds, +): + """Creates an animated 3D gif of the recorded history. + + You need imagemagick on your system for this function to work. + The draw_func must accept an `ax` keyword argument and plot into + the provided Axes3D instance. + + Parameters + ---------- + history : a :class:`tyssue.History` object + output : path to the output gif file + num_frames : int, the number of frames in the gif + interval : tuple, define begin and end frame of the gif + draw_func : a drawing function + Must take a `sheet` object as first argument and return a + `fig, ax` pair. Must accept an `ax` keyword argument so it + can plot into the pre-created Axes3D. Defaults to sheet_view. + margin : int, graph margins in percent, default 5. + If -1, let the draw function decide. + dpi : int, resolution of each saved frame, default 200 + view_angle : tuple (elev, azim), default (30, 45) + Elevation and azimuth angles for the 3D camera. + dynamic_draw_kwds : list of functions or None + list of functions that are called to update the draw_kds + + Example:: + + dynamic_draw_kwds={ + "face_colors": lambda sheet: sheet.face_df["myogen"].values, + } + + **draw_kwds are passed unchanged to the drawing function + """ + if draw_func is None: + draw_func = sheet_view_3d # default to the 3D view + + draw_kwds.setdefault("view_angle", view_angle) + + if dynamic_draw_kwds is None: + dynamic_draw_kwds = [] + + graph_dir = pathlib.Path(tempfile.mkdtemp()) + + if interval is None: + start, stop = None, None + else: + start, stop = interval[0], interval[1] + + coords = draw_kwds.get("coords", history.sheet.coords[:3]) + x, y, z = coords[0], coords[1], coords[2] + sheet0 = history.retrieve(0) + bounds = sheet0.vert_df[coords].describe().loc[["min", "max"]] + delta = (bounds.loc["max"] - bounds.loc["min"]).max() + margin_val = delta * margin / 100 + xlim = bounds.loc["min", x] - margin_val, bounds.loc["max", x] + margin_val + ylim = bounds.loc["min", y] - margin_val, bounds.loc["max", y] + margin_val + zlim = bounds.loc["min", z] - margin_val, bounds.loc["max", z] + margin_val + + for i, (t, sheet) in enumerate(history.browse(start, stop, num_frames)): + try: + if len(dynamic_draw_kwds) > 0: + for func in dynamic_draw_kwds: + update_kwds = func(sheet) + draw_kwds = deep_update(draw_kwds, update_kwds) + + fig = plt.figure() + ax = fig.add_subplot(111, projection="3d") + ax.view_init(elev=view_angle[0], azim=view_angle[1]) + + fig, ax = draw_func(sheet, ax=ax, legend=legend, cull_back_edges=cull_back_edges, **draw_kwds) + patch_2d_collections_to_3d(ax) + ax.set(xlim=xlim, ylim=ylim, zlim=zlim) + ax.set_title(f"t = {t:.2f}") + + except Exception as e: + print(f"Dropped frame {i}") + print(e) + continue + + fig.savefig( + graph_dir / f"movie_{i:04d}.png", + dpi=dpi, + bbox_inches="tight", + ) + plt.close(fig) -def sheet_view(sheet, coords=COORDS, ax=None, cbar_axis=None, **draw_specs_kw): + try: + subprocess.run(["magick", (graph_dir / "movie_*.png").as_posix(), output]) + except Exception as e: + print( + "Converting didn't work, make sure imagemagick is available on your system" + ) + raise e + + finally: + shutil.rmtree(graph_dir) + +def sheet_view(sheet, coords=COORDS, ax=None, cbar_axis=None, legend=None, **draw_specs_kw): """Base view function, parametrizable through draw_secs The default sheet_spec specification is: @@ -259,8 +407,140 @@ def sheet_view(sheet, coords=COORDS, ax=None, cbar_axis=None, **draw_specs_kw): cb1.set_label("a.u.") else: cb1.set_label(axis_spec.get("color_bar_label")) + + if legend is not None: + handles = [ + mpatches.Patch(color=color, label=label) + for label, color in legend.items() + ] + ax.legend(handles=handles, loc="upper left", bbox_to_anchor=(0, 1)) + return fig, ax +def sheet_view_3d(sheet, coords=COORDS, ax=None, view_angle=(30, 45), cull_back_edges=False, legend=None, draw_order=("face", "vert", "edge"), **draw_specs_kw): + """3D version of sheet_view using Axes3D. + + Parameters + ---------- + sheet : a tyssue Sheet object + coords : list of 3 coordinate names, default COORDS + ax : an Axes3D instance, or None to create a new one + view_angle : tuple (elev, azim), default (30, 45) + draw_order : tuple or list of {"face", "vert", "edge"}, default ("face", "vert", "edge") + Order in which the elements are drawn. Elements drawn later appear on + top. Any element omitted from this sequence is not drawn. + **draw_specs_kw : passed to the draw spec updater + """ + draw_specs = sheet_spec() + spec_updater(draw_specs, draw_specs_kw) + + if ax is None: + fig = plt.figure() + ax = fig.add_subplot(111, projection="3d") + else: + fig = ax.get_figure() + + ax.view_init(elev=view_angle[0], azim=view_angle[1]) + + valid_elements = {"face", "vert", "edge"} + unknown = set(draw_order) - valid_elements + if unknown: + raise ValueError( + f"Unknown element(s) in draw_order: {sorted(unknown)}. " + f"Valid elements are {sorted(valid_elements)}." + ) + + for element in draw_order: + spec = draw_specs[element] + if not spec["visible"]: + continue + if element == "face": + ax = draw_face_3d(sheet, coords, ax, **spec) + elif element == "vert": + ax = draw_vert_3d(sheet, coords, ax, **spec) + elif element == "edge": + ax = draw_edge_3d(sheet, coords, ax, view_angle=view_angle, cull_back_edges=cull_back_edges, **spec) + + if legend is not None: + handles = [ + mpatches.Patch(color=color, label=label) + for label, color in legend.items() + ] + ax.legend(handles=handles, loc="upper left", bbox_to_anchor=(0, 1)) + + ax.autoscale() + _set_axes_proportional_3d(ax) + _auto_tick_fontsize_3d(ax, base_size=8, min_size=4) + return fig, ax + +def draw_faces_highlighted( + sheet, + face_indices, + highlight_color, + coords=COORDS, + ax=None, + alpha=1.0, + background_color=(1.0, 1.0, 1.0, 1.0), + show_edges=True, +): + """ + Draw full tissue, highlighting selected faces in a given color + and rendering all others in white. + + Parameters + ---------- + sheet : Sheet + face_indices : array-like + Indices of faces to highlight + highlight_color : color-like + Matplotlib color (e.g. "#ff0000", "red", RGBA) + coords : tuple + ax : matplotlib axis, optional + alpha : float + background_color : RGBA tuple + show_edges : bool + """ + + from matplotlib.colors import to_rgba + + # Convert highlight color to RGBA + hi_rgba = np.array(to_rgba(highlight_color)) + bg_rgba = np.array(background_color) + + # Build per-face RGBA array + face_colors = np.tile(bg_rgba, (sheet.Nf, 1)) + + face_idx = sheet.face_df.index + mask = face_idx.isin(face_indices) + + face_colors[mask] = hi_rgba + face_colors[:, 3] *= alpha # apply alpha uniformly + + draw_specs = { + "face": { + "visible": True, + "color": face_colors, + }, + "edge": { + "visible": show_edges, + }, + "vert": { + "visible": False, + }, + "axis": { + "autoscale": True, + "color_bar": False, + }, + } + + fig, ax = sheet_view( + sheet, + coords=coords, + ax=ax, + **draw_specs, + ) + + return fig, ax def draw_face(sheet, coords, ax, **draw_spec_kw): """Draws epithelial sheet polygonal faces in matplotlib @@ -295,6 +575,125 @@ def draw_face(sheet, coords, ax, **draw_spec_kw): return ax +def draw_vert(sheet, coords, ax, **draw_spec_kw): + """Draw junction vertices in matplotlib.""" + draw_spec = sheet_spec()["vert"] + draw_spec.update(**draw_spec_kw) + + x, y = coords + if "z_coord" in sheet.vert_df.columns: + pos = sheet.vert_df.sort_values("z_coord")[coords] + else: + pos = sheet.vert_df[coords] + ax.scatter(pos[x], pos[y], **draw_spec_kw) + return ax + + +def draw_edge(sheet, coords, ax, **draw_spec_kw): + """""" + draw_spec = sheet_spec()["edge"] + draw_spec.update(**draw_spec_kw) + arrow_specs, collections_specs = _parse_edge_specs(draw_spec, sheet) + dx, dy = ("d" + c for c in coords) + sx, sy = ("s" + c for c in coords) + tx, ty = ("t" + c for c in coords) + + if draw_spec.get("head_width"): + + app_length = ( + np.hypot(sheet.edge_df[dx], sheet.edge_df[dy]) * sheet.edge_df.length.mean() + ) + patches = [ + FancyArrow(*edge[[sx, sy, dx, dy]], **arrow_specs) + for idx, edge in sheet.edge_df[app_length > 1e-6].iterrows() + ] + ax.add_collection( + PatchCollection(patches, match_original=False, **collections_specs) + ) + else: + segments = sheet.edge_df[[sx, sy, tx, ty]].to_numpy().reshape((-1, 2, 2)) + ax.add_collection(LineCollection(segments, **collections_specs)) + return ax + + +def draw_vert_3d(sheet, coords, ax, **draw_spec_kw): + """Draw junction vertices in 3D matplotlib.""" + draw_spec = sheet_spec()["vert"] + draw_spec.update(**draw_spec_kw) + + x, y, z = coords + if "z_coord" in sheet.vert_df.columns: + pos = sheet.vert_df.sort_values("z_coord")[coords] + else: + pos = sheet.vert_df[coords] + + ax.scatter(pos[x], pos[y], pos[z], **draw_spec_kw) + ax.autoscale() + return ax + + +def draw_edge_3d(sheet, coords, ax, view_angle=(30, 45), cull_back_edges=False, **draw_spec_kw): + draw_spec = sheet_spec()["edge"] + draw_spec.update(**draw_spec_kw) + _, collections_specs = _parse_edge_specs(draw_spec, sheet) + + sx, sy, sz = ("s" + c for c in coords) + tx, ty, tz = ("t" + c for c in coords) + + edge_df = sheet.edge_df + + if cull_back_edges: + mx = (edge_df[sx] + edge_df[tx]) / 2 + my = (edge_df[sy] + edge_df[ty]) / 2 + + # Only use azimuth — culling is purely in xy for a z-axis cylinder + azim = np.deg2rad(view_angle[1]) + view_dir_xy = np.array([np.cos(azim), np.sin(azim)]) + + # Centroid in xy only + cx, cy = mx.mean(), my.mean() + outward_xy = np.stack([mx - cx, my - cy], axis=1) + + dots = outward_xy @ view_dir_xy + edge_df = edge_df[dots > 0] + + segments = ( + edge_df[[sx, sy, sz, tx, ty, tz]] + .to_numpy() + .reshape((-1, 2, 3)) + ) + ax.add_collection3d(Line3DCollection(segments, **collections_specs)) + return ax + +def draw_face_3d(sheet, coords, ax, **draw_spec_kw): + """Draw epithelial sheet polygonal faces as a Poly3DCollection.""" + draw_spec = sheet_spec()["face"] + draw_spec.update(**draw_spec_kw) + collection_specs = parse_face_specs(draw_spec, sheet) + + if "visible" in sheet.face_df.columns: + edges = sheet.edge_df[sheet.upcast_face(sheet.face_df["visible"])].index + if edges.shape[0]: + _sheet = get_sub_eptm(sheet, edges) + sheet = _sheet + color = collection_specs["facecolors"] + if isinstance(color, np.ndarray): + faces = sheet.face_df["face_o"].values.astype(np.uint32) + collection_specs["facecolors"] = color.take(faces, axis=0) + else: + warnings.warn("No face is visible") + + if not sheet.is_ordered: + sheet_ = sheet.copy() + sheet_.reset_index(order=True) + polys = sheet_.face_polygons(coords) + else: + polys = sheet.face_polygons(coords) + + p = Poly3DCollection(polys, closed=True, **collection_specs) + ax.add_collection3d(p) + return ax + def parse_face_specs(face_draw_specs, sheet): collection_specs = {} @@ -340,47 +739,6 @@ def _face_color_from_sequence(face_spec, sheet): ) -def draw_vert(sheet, coords, ax, **draw_spec_kw): - """Draw junction vertices in matplotlib.""" - draw_spec = sheet_spec()["vert"] - draw_spec.update(**draw_spec_kw) - - x, y = coords - if "z_coord" in sheet.vert_df.columns: - pos = sheet.vert_df.sort_values("z_coord")[coords] - else: - pos = sheet.vert_df[coords] - ax.scatter(pos[x], pos[y], **draw_spec_kw) - return ax - - -def draw_edge(sheet, coords, ax, **draw_spec_kw): - """""" - draw_spec = sheet_spec()["edge"] - draw_spec.update(**draw_spec_kw) - arrow_specs, collections_specs = _parse_edge_specs(draw_spec, sheet) - dx, dy = ("d" + c for c in coords) - sx, sy = ("s" + c for c in coords) - tx, ty = ("t" + c for c in coords) - - if draw_spec.get("head_width"): - - app_length = ( - np.hypot(sheet.edge_df[dx], sheet.edge_df[dy]) * sheet.edge_df.length.mean() - ) - patches = [ - FancyArrow(*edge[[sx, sy, dx, dy]], **arrow_specs) - for idx, edge in sheet.edge_df[app_length > 1e-6].iterrows() - ] - ax.add_collection( - PatchCollection(patches, match_original=False, **collections_specs) - ) - else: - segments = sheet.edge_df[[sx, sy, tx, ty]].to_numpy().reshape((-1, 2, 2)) - ax.add_collection(LineCollection(segments, **collections_specs)) - return ax - - def _parse_edge_specs(edge_draw_specs, sheet): arrow_keys = ["head_width", "length_includes_head", "shape"] @@ -472,6 +830,23 @@ def _get_lines(sheet, coords): return lines_x, lines_y +def _set_axes_proportional_3d(ax): + x_range = ax.get_xlim3d()[1] - ax.get_xlim3d()[0] + y_range = ax.get_ylim3d()[1] - ax.get_ylim3d()[0] + z_range = ax.get_zlim3d()[1] - ax.get_zlim3d()[0] + ax.set_box_aspect([x_range, y_range, z_range]) + +def _auto_tick_fontsize_3d(ax, base_size=8, min_size=4): + ranges = np.array([ + ax.get_xlim3d()[1] - ax.get_xlim3d()[0], + ax.get_ylim3d()[1] - ax.get_ylim3d()[0], + ax.get_zlim3d()[1] - ax.get_zlim3d()[0], + ]) + max_range = ranges.max() + size = max(min_size, round(base_size * min(ranges) / max_range)) + for ax_obj in [ax.xaxis, ax.yaxis, ax.zaxis]: + ax_obj.set_tick_params(labelsize=size) + def plot_forces( sheet, geom, model, coords, scaling, ax=None, approx_grad=None, **draw_specs_kw ): diff --git a/src/tyssue/dynamics/__init__.py b/src/tyssue/dynamics/__init__.py index cd9b2a98..6aaa3d09 100644 --- a/src/tyssue/dynamics/__init__.py +++ b/src/tyssue/dynamics/__init__.py @@ -1,5 +1,5 @@ """dynamics""" from .bulk_model import BulkModel, LaminaModel # noqa -from .factory import model_factory # noqa +from .factory import model_factory, model_factory_vessel, model_factory_cylinder # noqa from .planar_vertex_model import PlanarModel # noqa from .sheet_vertex_model import SheetModel # noqa diff --git a/src/tyssue/dynamics/effectors.py b/src/tyssue/dynamics/effectors.py index b0e2bfef..8aacf1ea 100644 --- a/src/tyssue/dynamics/effectors.py +++ b/src/tyssue/dynamics/effectors.py @@ -1,21 +1,30 @@ """ Generic forces and energies """ -import numpy as np import pandas as pd +import numpy as np from ..utils import to_nd from . import units -from .bulk_gradients import lumen_volume_grad, volume_grad + from .planar_gradients import area_grad as area_grad2d from .planar_gradients import lumen_area_grad -from .sheet_gradients import area_grad, height_grad +from .sheet_gradients import height_grad, area_grad +from .bulk_gradients import volume_grad, lumen_volume_grad def elastic_force(element_df, var, elasticity, prefered): - params = {"x": var, "K": elasticity, "x0": prefered} - force = element_df.eval("{K} * ({x} - {x0})".format(**params)) - return force + """ + K can be: + - column name (str) + - callable: f(df) -> array + """ + if callable(elasticity): + K_val = elasticity(element_df) + else: + K_val = element_df[elasticity] + + return K_val * (element_df[var] - element_df[prefered]) def _elastic_force(element_df, x, elasticity, prefered): @@ -24,8 +33,20 @@ def _elastic_force(element_df, x, elasticity, prefered): def elastic_energy(element_df, var, elasticity, prefered): - params = {"x": var, "K": elasticity, "x0": prefered} - energy = element_df.eval("0.5 * {K} * ({x} - {x0}) ** 2".format(**params)) + """ + elasticity can be: + - str: column name + - callable: f(df) -> array-like + """ + x = element_df[var] + x0 = element_df[prefered] + + if callable(elasticity): + K = elasticity(element_df) + else: + K = element_df[elasticity] + + energy = 0.5 * K * (x - x0) ** 2 return energy @@ -35,7 +56,7 @@ def _elastic_energy(element_df, x, elasticity, prefered): class AbstractEffector: - """The effector class is used by model factories + """ The effector class is used by model factories to construct a model. @@ -71,50 +92,10 @@ def get_nrj_norm(specs): # Works on an `Epithelium` object's {cls.element} elements. # """ -class Repulsion(AbstractEffector): - """ - Repulsion to avoid intersection between two cells. - Effector for 2D lateral model. - """ - dimensions = units.line_elasticity - magnitude = "cell_repulsion" - label = "Cell Repulsion" - element = "vert" - specs = { - "vert": {"force_repulsion": 1.0, - "v_repulsion": 0.0} - } - - @staticmethod - def energy(eptm): - grid = eptm.vert_df.loc[0, "grid"][0] - x = np.argmin(np.abs([grid[0][:, 0] - x for x in eptm.vert_df["x"]]), axis=1) - y = np.argmin(np.abs([grid[1][0, :] - y for y in eptm.vert_df["y"]]), axis=1) - repulse = [eptm.vert_df.loc[v, "v_repulsion"][0][x[v], y[v]] for v in range(eptm.Nv)] - return np.array(eptm.specs['vert']["force_repulsion"]) * repulse - - @staticmethod - def gradient(eptm): - repulse_u = [] - repulse_v = [] - grid = eptm.vert_df.loc[0, "grid"][0] - - x = np.argmin(np.abs([grid[0][:, 0] - x for x in eptm.vert_df["x"]]), axis=1) - y = np.argmin(np.abs([grid[1][0, :] - y for y in eptm.vert_df["y"]]), axis=1) - for v in range(eptm.Nv): - U, V = np.gradient(eptm.vert_df.loc[v, "v_repulsion"][0], 1, 1) - repulse_u.append(U[x[v], y[v]]) - repulse_v.append(V[x[v], y[v]]) - - grad = np.array(eptm.specs['vert']["force_repulsion"]) * pd.DataFrame(np.array([repulse_u, repulse_v]).T) - grad.columns = ["g" + c for c in eptm.coords] - return grad, None class LengthElasticity(AbstractEffector): - """ - Elastic half edge elasticity using the formula - ..math: 1/2*length_elasticity*(length-prefered_length)**2 + """Elastic half edge """ dimensions = units.line_elasticity @@ -138,19 +119,19 @@ class LengthElasticity(AbstractEffector): @staticmethod def get_nrj_norm(specs): return ( - specs["edge"]["length_elasticity"] * specs["edge"]["prefered_length"] ** 2 + specs["edge"]["length_elasticity"] * specs["edge"]["prefered_length"] ** 2 ) @staticmethod def energy(eptm): return elastic_energy( - eptm.edge_df, "length", "length_elasticity * is_active", "prefered_length" + eptm.edge_df, "length", lambda df: df["length_elasticity"] * df["is_alive"], "prefered_length" ) @staticmethod def gradient(eptm): kl_l0 = elastic_force( - eptm.edge_df, "length", "length_elasticity * is_active", "prefered_length" + eptm.edge_df, "length", lambda df: df["length_elasticity"] * df["is_alive"], "prefered_length" ) grad = eptm.edge_df[eptm.ucoords] * to_nd(kl_l0, eptm.dim) grad.columns = ["g" + u for u in eptm.coords] @@ -158,11 +139,9 @@ def gradient(eptm): class PerimeterElasticity(AbstractEffector): + """From Mapeng Bi et al. https://doi.org/10.1038/nphys3471 """ - Face perimeter elasticity using the formula - ..math: 1/2*perimeter_elasticity*(perimeter-prefered_perimeter)**2 - From Mapeng Bi et al. https://doi.org/10.1038/nphys3471 - """ + dimensions = units.line_elasticity magnitude = "perimeter_elasticity" label = "Perimeter Elasticity" @@ -180,34 +159,33 @@ class PerimeterElasticity(AbstractEffector): @staticmethod def energy(eptm): - return elastic_energy( - eptm.face_df, - "perimeter", - "perimeter_elasticity * is_alive", - "prefered_perimeter", - ) + df = eptm.face_df + diff = df["perimeter"] - df["prefered_perimeter"] + return 0.5 * df["is_alive"] * df["perimeter_elasticity"] * diff * diff @staticmethod def gradient(eptm): - gamma_ = elastic_force( - eptm.face_df, - "perimeter", - "perimeter_elasticity * is_alive", - "prefered_perimeter", - ) + # Compute gamma directly + df = eptm.face_df + gamma_ = df["perimeter_elasticity"] * df["is_alive"] * (df["perimeter"] - df["prefered_perimeter"]) + + # Upcast gamma gamma = eptm.upcast_face(gamma_) - grad_srce = -eptm.edge_df[eptm.ucoords] * to_nd(gamma, len(eptm.coords)) - grad_srce.columns = ["g" + u for u in eptm.coords] + # Convert gamma to node-level array + gamma_nd = to_nd(gamma, len(eptm.coords)) + + # Compute gradient at edges + grad_srce = -eptm.edge_df[eptm.ucoords].to_numpy() * gamma_nd + grad_srce = pd.DataFrame(grad_srce, columns=["g" + u for u in eptm.coords]) + + # grad_trgt is just the negative grad_trgt = -grad_srce return grad_srce, grad_trgt class FaceAreaElasticity(AbstractEffector): - """ - Face area elasticity using the formula - ..math: 1/2*area_elasticity*(area-prefered_area)**2 - """ + dimensionless = False dimensions = units.area_elasticity magnitude = "area_elasticity" @@ -232,13 +210,13 @@ def get_nrj_norm(specs): @staticmethod def energy(eptm): return elastic_energy( - eptm.face_df, "area", "area_elasticity * is_alive", "prefered_area" + eptm.face_df, "area", lambda df: df["area_elasticity"] * df["is_alive"], "prefered_area" ) @staticmethod def gradient(eptm): ka_a0_ = elastic_force( - eptm.face_df, "area", "area_elasticity * is_alive", "prefered_area" + eptm.face_df, "area", lambda df: df["area_elasticity"] * df["is_alive"], "prefered_area" ) ka_a0 = to_nd(eptm.upcast_face(ka_a0_), len(eptm.coords)) @@ -257,12 +235,7 @@ def gradient(eptm): class FaceVolumeElasticity(AbstractEffector): - """ - Face volume elasticity using the formula - ..math: 1/2*volume_elasticity*(volume-prefered_volume)**2 - Effector for 2.5D model, where a volume of a cell is taking into account where only apical surface is modeled - """ dimensions = units.vol_elasticity magnitude = "vol_elasticity" label = "Volume elasticity" @@ -282,13 +255,13 @@ def get_nrj_norm(specs): @staticmethod def energy(eptm): return elastic_energy( - eptm.face_df, "vol", "vol_elasticity * is_alive", "prefered_vol" + eptm.face_df, "vol", lambda df: df["vol_elasticity"] * df["is_alive"], "prefered_vol" ) @staticmethod def gradient(eptm): kv_v0_ = elastic_force( - eptm.face_df, "vol", "vol_elasticity * is_alive", "prefered_vol" + eptm.face_df, "vol", lambda df: df["vol_elasticity"] * df["is_alive"], "prefered_vol" ) kv_v0 = to_nd(eptm.upcast_face(kv_v0_), 3) @@ -309,10 +282,7 @@ def gradient(eptm): class CellAreaElasticity(AbstractEffector): - """ - Cell area elasticity using the formula - ..math: 1/2*area_elasticity*(area-prefered_area)**2 - """ + dimensions = units.area_elasticity magnitude = "area_elasticity" label = "Area elasticity" @@ -338,7 +308,7 @@ def energy(eptm): @staticmethod def gradient(eptm): ka_a0_ = elastic_force( - eptm.cell_df, "area", "area_elasticity * is_alive", "prefered_area" + eptm.cell_df, "area", lambda df: df["area_elasticity"] * df["is_alive"], "prefered_area" ) ka_a0 = to_nd(eptm.upcast_cell(ka_a0_), 3) @@ -354,10 +324,7 @@ def gradient(eptm): class CellVolumeElasticity(AbstractEffector): - """ - Cell volume elasticity using the formula - ..math: 1/2*volumne_elasticity*(volume-prefered_volume)**2 - """ + dimensions = units.vol_elasticity magnitude = "vol_elasticity" label = "Volume elasticity" @@ -379,7 +346,7 @@ def energy(eptm): @staticmethod def gradient(eptm): kv_v0_ = elastic_force( - eptm.cell_df, "vol", "vol_elasticity * is_alive", "prefered_vol" + eptm.cell_df, "vol", lambda df: df["vol_elasticity"] * df["is_alive"], "prefered_vol" ) kv_v0 = to_nd(eptm.upcast_cell(kv_v0_), 3) @@ -395,9 +362,7 @@ def gradient(eptm): class LumenVolumeElasticity(AbstractEffector): """ - Global volume elasticity of the object. using the formula - ..math: 1/2*lumen_elasticity*(lumen-prefered_lumen)**2 - + Global volume elasticity of the object. For example the volume of the yolk in the Drosophila embryo """ @@ -418,12 +383,13 @@ class LumenVolumeElasticity(AbstractEffector): @staticmethod def get_nrj_norm(specs): return ( - specs["settings"]["lumen_vol_elasticity"] - * specs["settings"]["lumen_prefered_vol"] ** 2 + specs["settings"]["lumen_vol_elasticity"] + * specs["settings"]["lumen_prefered_vol"] ** 2 ) @staticmethod def energy(eptm): + return _elastic_energy( eptm.settings, "lumen_vol", "lumen_vol_elasticity", "lumen_prefered_vol" ) @@ -445,10 +411,7 @@ def gradient(eptm): class LineTension(AbstractEffector): - """ - Half edge line tension using the formula - ..math: line_tension*length/2 - """ + dimensions = units.line_tension magnitude = "line_tension" label = "Line tension" @@ -459,25 +422,26 @@ class LineTension(AbstractEffector): @staticmethod def energy(eptm): - return eptm.edge_df.eval( - "line_tension" "* is_active" "* length / 2" - ) # accounts for half edges + df = eptm.edge_df + return 0.5 * df["line_tension"] * df["is_active"] * df["length"] # accounts for half edges @staticmethod def gradient(eptm): - grad_srce = -eptm.edge_df[eptm.ucoords] * to_nd( - eptm.edge_df.eval("line_tension * is_active/2"), len(eptm.coords) + edge_df = eptm.edge_df + + coeff = (edge_df["line_tension"] * edge_df["is_active"]) * 0.5 + + grad_srce = -edge_df[eptm.ucoords] * to_nd( + coeff, len(eptm.coords) ) + grad_srce.columns = ["g" + u for u in eptm.coords] grad_trgt = -grad_srce return grad_srce, grad_trgt class FaceContractility(AbstractEffector): - """ - Face contractility using the formula - ..math: 1/2*contractility*perimeter**2 - """ + dimensions = units.line_elasticity magnitude = "contractility" label = "Contractility" @@ -488,24 +452,32 @@ class FaceContractility(AbstractEffector): @staticmethod def energy(eptm): - return eptm.face_df.eval("0.5 * is_alive * contractility * perimeter ** 2") + df = eptm.face_df + return 0.5 * df["is_alive"] * df["contractility"] * df["perimeter"] * df["perimeter"] @staticmethod def gradient(eptm): - gamma_ = eptm.face_df.eval("contractility * perimeter * is_alive") + # Compute gamma directly + df = eptm.face_df + gamma_ = df["contractility"] * df["perimeter"] * df["is_alive"] + + # Upcast gamma to edges gamma = eptm.upcast_face(gamma_) - grad_srce = -eptm.edge_df[eptm.ucoords] * to_nd(gamma, len(eptm.coords)) - grad_srce.columns = ["g" + u for u in eptm.coords] + # Convert gamma to node-level array + gamma_nd = to_nd(gamma, len(eptm.coords)) + + # Compute gradient at edges + grad_srce = -eptm.edge_df[eptm.ucoords].to_numpy() * gamma_nd + grad_srce = pd.DataFrame(grad_srce, columns=["g" + u for u in eptm.coords]) + + # grad_trgt is just the negative grad_trgt = -grad_srce return grad_srce, grad_trgt class SurfaceTension(AbstractEffector): - """ - Face surface tension using the formula - ..math: surface_tension*area - """ + dimensions = units.area_tension magnitude = "surface_tension" @@ -517,10 +489,12 @@ class SurfaceTension(AbstractEffector): @staticmethod def energy(eptm): - return eptm.face_df.eval("surface_tension * area") + df = eptm.face_df + return df["surface_tension"] * df["area"] @staticmethod def gradient(eptm): + G = to_nd(eptm.upcast_face(eptm.face_df["surface_tension"]), len(eptm.coords)) grad_a_srce, grad_a_trgt = area_grad(eptm) @@ -533,9 +507,7 @@ def gradient(eptm): class LineViscosity(AbstractEffector): - """ - Edge line viscosity - """ + dimensions = units.line_viscosity magnitude = "edge_viscosity" @@ -555,10 +527,6 @@ def gradient(eptm): class BorderElasticity(AbstractEffector): - """ - Edge border elasticity using the formula - ..math: border_elasticity*prefered_length**2 - """ dimensions = units.line_elasticity label = "Border edges elasticity" magnitude = "border_elasticity" @@ -578,7 +546,7 @@ class BorderElasticity(AbstractEffector): @staticmethod def get_nrj_norm(specs): return ( - specs["edge"]["border_elasticity"] * specs["edge"]["prefered_length"] ** 2 + specs["edge"]["border_elasticity"] * specs["edge"]["prefered_length"] ** 2 ) @staticmethod @@ -586,16 +554,17 @@ def energy(eptm): return elastic_energy( eptm.edge_df, "length", - "border_elasticity * is_active * is_border / 2", + lambda df: df["border_elasticity"] * df["is_active"] * df["is_border"]/2, "prefered_length", ) @staticmethod def gradient(eptm): + kl_l0 = elastic_force( eptm.edge_df, var="length", - elasticity="border_elasticity * is_active * is_border", + elasticity= lambda df: df["border_elasticity"] * df["is_active"] * df["is_border"], prefered="prefered_length", ) grad = eptm.edge_df[eptm.ucoords] * to_nd(kl_l0, eptm.dim) @@ -604,11 +573,6 @@ def gradient(eptm): class LumenAreaElasticity(AbstractEffector): - """ - Lumen area elasticity for 2D simulation using the formula - ..math: \frac{K_Y}{2}(A_{\mathrm{lumen}} - A_{0,\mathrm{lumen}})^2 - - """ dimensions = units.area_elasticity label = "Lumen volume constraint" @@ -642,8 +606,7 @@ def gradient(eptm): class RadialTension(AbstractEffector): """ - Apply a tension perpendicular to a face divide equally on each vertex - ..math: height*radialTension + Apply a tension perpendicular to a face. """ dimensions = units.line_tension @@ -654,12 +617,13 @@ class RadialTension(AbstractEffector): @staticmethod def energy(eptm): - return eptm.face_df.eval("height * radial_tension") + df = eptm.face_df + return df["height"] * df["radial_tension"] @staticmethod def gradient(eptm): upcast_tension = eptm.upcast_face( - eptm.face_df.eval("radial_tension / num_sides") + eptm.face_df["radial_tension"] / eptm.face_df["num_sides"] ) upcast_height = eptm.upcast_srce(height_grad(eptm)) @@ -670,8 +634,7 @@ def gradient(eptm): class BarrierElasticity(AbstractEffector): """ - Barrier use to maintain the tissue integrity, for 2.5D geometry - ..math: \frac{1}{2} K_barrier \detha_\rho^2 + Barrier use to maintain the tissue integrity. """ dimensions = units.line_elasticity @@ -684,65 +647,216 @@ class BarrierElasticity(AbstractEffector): @staticmethod def energy(eptm): - return eptm.vert_df.eval("delta_rho**2 * barrier_elasticity/2") + df = eptm.vert_df + return 0.5 * df["barrier_elasticity"] * df["delta_rho"] * df["delta_rho"] @staticmethod def gradient(eptm): - grad = height_grad(eptm) * to_nd( - eptm.vert_df.eval("barrier_elasticity * delta_rho"), 3 - ) + # Compute the vertex-level factor + df = eptm.vert_df + factor = df["barrier_elasticity"] * df["delta_rho"] + + # Convert to node-level array + factor_nd = to_nd(factor, 3) + + # Compute gradient + grad = height_grad(eptm) * factor_nd grad.columns = ["g" + c for c in eptm.coords] + return grad, None -class MidlineBoundary(AbstractEffector): - """ - Elastic boundary at the x-axis, to be used with MidlineBoundaryGeometry. - Intended to use with a high midline_boundary_stiffness to prevent vertices - from crossing the midline. +class ChiralTorque(AbstractEffector): + + dimensions = units.line_elasticity + magnitude = "torque_coef" + label = "Apply Chiral Torque to Cells" + element = "face" + specs = { + "face": {"torque_coef": 0.0, "is_alive": 1} + } + + @staticmethod + def energy(eptm): + return np.zeros(eptm.Nv) + + @staticmethod + def gradient(eptm): + torque = eptm.face_df['torque_coef'] + torque = to_nd(eptm.upcast_face(torque), len(eptm.coords)) + grad_srce = np.multiply(eptm.edge_df[["r" + z for z in eptm.coords]].values, torque) + normal = eptm.edge_df[["n" + u for u in eptm.coords]].values + grad_srce = np.cross(grad_srce, normal) + srce_active = eptm.upcast_srce(eptm.vert_df['is_active']) + grad_srce = grad_srce * \ + to_nd(srce_active, len(eptm.coords)) #* srce_bound_coords + grad_srce = pd.DataFrame(grad_srce) + grad_srce.columns = ["g" + u for u in eptm.coords] + + return grad_srce, None + + +class ActiveMigration(AbstractEffector): + """Active cell migration force along a specified direction. + + This is a non-conservative force that drives cells to migrate + along a specified vector direction with constant magnitude. """ + dimensions = units.force # or appropriate force units + magnitude = "migration_strength" + label = "Active Migration" + element = "face" + specs = { + "face": { + "is_alive": 1, + "migration_strength": 0.1, # magnitude of migration force + "migration_dir_x": 1.0, # x-component of migration direction + "migration_dir_y": 0.0, # y-component of migration direction + "migration_dir_z": 0.0, # z-component (if 3D) + } + } + + spatial_ref = "migration_strength", units.force + + @staticmethod + def energy(eptm): + """Non-conservative force - return zero energy.""" + return np.zeros(eptm.Nv) + + @staticmethod + def gradient(eptm): + """Compute migration forces on vertices. + + Each cell exerts a constant force in its migration direction, + distributed among its vertices. + """ + df = eptm.face_df + + # Get migration force vector for each face + # Force magnitude scaled by migration_strength and is_alive + force_magnitude = df["migration_strength"] * df["is_alive"] + + # Build force vector (normalize direction first) + migration_dirs = df[[f"m{c}" for c in eptm.coords]].to_numpy() + + # Normalize direction vectors (per face) + norms = np.linalg.norm(migration_dirs, axis=1, keepdims=True) + norms = np.where(norms > 0, norms, 1.0) # avoid division by zero + migration_dirs_normalized = migration_dirs / norms + + # Scale by force magnitude + force_vectors = migration_dirs_normalized * force_magnitude.to_numpy()[:, np.newaxis] + + # Upcast from face to edge level + force_per_coord = { + coord: eptm.upcast_face(force_vectors[:, i]) + for i, coord in enumerate(eptm.coords) + } + + # Convert to node-level arrays + force_nd = np.column_stack([ + to_nd(force_per_coord[coord], len(eptm.coords)) + for coord in eptm.coords + ]) + + # Distribute force equally to source vertices + # (could also distribute based on edge length or other schemes) + grad_srce = pd.DataFrame(force_nd, columns=["g" + u for u in eptm.coords]) + + return -grad_srce, None + + +class SurfaceElasticity(AbstractEffector): + dimensions = units.line_elasticity - magnitude = "midline_boundary" - label = "Midline boundary" + magnitude = "surface_elasticity" + label = "Apply Surface Elasticity to vertices such that a flat surface is prefered" element = "vert" - specs = {"vert": {"boundary_K": 280, "is_active": 1, "delta_boundary": 0}} + specs = { + "vert": {"torque_coef": 0.0, "is_alive": 1} + } @staticmethod def energy(eptm): - return eptm.vert_df.eval( - "0.5 * delta_boundary**2 * {}".format( - str(eptm.settings["midline_boundary_stiffness"]) - ) + return elastic_energy( + eptm.vert_df, "dev_length", lambda df: df["surface_elasticity"] * df["is_active"], "prefered_deviation" ) @staticmethod def gradient(eptm): - # just a bunch of zeros - kl_l0 = elastic_force(eptm.vert_df, "delta_boundary", "0", "0") - grad = eptm.vert_df[eptm.coords] * to_nd(kl_l0, eptm.dim) + ka_a0_ = elastic_force( + eptm.vert_df, "dev_length", lambda df: df["surface_elasticity"] * df["is_active"], "prefered_deviation" + ) + + ka_a0 = to_nd(ka_a0_, len(eptm.coords)) + + grad = eptm.vert_df[["d" + x for x in eptm.coords]].to_numpy() + + grad = pd.DataFrame(grad * ka_a0) + grad.columns = ["g" + u for u in eptm.coords] - grad["gy"] = 0 - if "z" in eptm.coords: - grad["gz"] = 0 - return grad, grad + return grad, None +class VesselSurfaceElasticity(AbstractEffector): + """ + Applies an elastic force to maintain vertex distance from xy origin + """ + + dimensions = units.line_elasticity + magnitude = "surface_elasticity" + label = "Apply Surface Elasticity to vertices such that a flat surface is prefered" + element = "vert" + specs = { + "vert": {"torque_coef": 0.0, "is_alive": 1} + } + + @staticmethod + def energy(eptm): + return elastic_energy( + eptm.vert_df, "distance_origin", lambda df: df["vessel_elasticity"] * df["is_alive"], "prefered_radius" + ) + + @staticmethod + def gradient(eptm): + coords = ["x", "y", "z"] + if "axis" in eptm.settings.keys(): + axis = eptm.settings["axis"] + else: + axis = "z" + coords.remove(axis) + ka_a0_ = elastic_force( + eptm.vert_df, "distance_origin", lambda df: df["vessel_elasticity"] * df["is_alive"], "prefered_radius" + ) + + ka_a0 = to_nd(ka_a0_, len(eptm.coords)) + + grad = eptm.vert_df[["o" + x for x in coords]].copy() + grad["o" + f"{axis}"] = 0 + grad = grad.to_numpy() + + grad = pd.DataFrame(grad * ka_a0) + + grad.columns = ["g" + u for u in eptm.coords] + + return grad, None def _exponants(dimensions, ref_dimensions, spatial_unit=None, temporal_unit=None): + spatial_exponant = time_exponant = 0 rel_dimensionality = (dimensions / ref_dimensions).dimensionality if spatial_unit is not None: spatial_exponant = ( - rel_dimensionality.get(units.length, 0) - / spatial_unit.dimensionality[units.length] + rel_dimensionality.get(units.length, 0) + / spatial_unit.dimensionality[units.length] ) if temporal_unit is not None: time_exponant = ( - rel_dimensionality.get(units.time, 0) - / temporal_unit.dimensionality[units.time] + rel_dimensionality.get(units.time, 0) + / temporal_unit.dimensionality[units.time] ) return spatial_exponant, time_exponant @@ -758,9 +872,9 @@ def scaler(nondim_specs, dim_specs, effector, ref_effector): ref_magnitude = ref_effector.magnitude ref_element = ref_effector.element factor = ( - dim_specs[ref_element][ref_magnitude] - * dim_specs[ref_element].get(spatial_val, 1) ** s_expo - * dim_specs[ref_element].get(temporal_val, 1) ** t_expo + dim_specs[ref_element][ref_magnitude] + * dim_specs[ref_element].get(spatial_val, 1) ** s_expo + * dim_specs[ref_element].get(temporal_val, 1) ** t_expo ) return factor diff --git a/src/tyssue/dynamics/factory.py b/src/tyssue/dynamics/factory.py index bbb60b3d..0a753f32 100644 --- a/src/tyssue/dynamics/factory.py +++ b/src/tyssue/dynamics/factory.py @@ -1,6 +1,9 @@ import warnings from copy import deepcopy +import pandas as pd +import numpy as np + from .effectors import dimensionalize as dimensionalize from .effectors import normalize as normalize @@ -122,3 +125,373 @@ def compute_gradient(eptm, components=False): return grad_i / norm_factor return NewModel + +def model_factory_vessel(effectors, ref_effector=None): + """Produces a Model class with the provided effectors. + + Parameters + ---------- + effectors : list of :class:`.effectors.AbstractEffectors` classes. + ref_effector : optional, default None + if passed, will be used for normalization, + by default, the last effector in the list is used + + Returns + ------- + NewModel : a Model derived class with compute_enregy and compute_gradient + methods + + """ + if ref_effector is None: + ref_effector = effectors[-1] + + class NewModel: + + labels = [] + specs = { + "cell": {}, + "face": {}, + "edge": {}, + "vert": {}, + "settings": {"nrj_norm_factor": 1.0}, + } + + _effectors = effectors + + for f in effectors: + labels.append(f.label) + try: + for k in specs: + specs[k].update(f.specs.get(k, {})) + except ValueError: + warnings.warn( + """ +Since 0.7, you need to provide a default value for each of the +specs parameters, e.g. + specs = { + "face": { + "perimeter": 1.0, + "perimeter_elasticity": 0.1, + "prefered_perimeter": 3.81, + } + } + +Setting all default values to 1.0 for now +""" + ) + for k in specs: + specs[k].update({key: 1.0 for key in f.specs.get(k, {})}) + + @staticmethod + def dimensionalize(nondim_specs): + dim_specs = deepcopy(nondim_specs) + for effector in effectors: + if effector == ref_effector: + continue + dimensionalize(nondim_specs, dim_specs, effector, ref_effector) + + ref_nrj = ref_effector.get_nrj_norm(dim_specs) + dim_specs["settings"]["nrj_norm_factor"] = ref_nrj + return dim_specs + + @classmethod + def dimentionalize(cls, nondim_specs): + warnings.warn( + """This badly worded method is deprecated, + use dimensionalize instead""" + ) + return cls.dimensionalize(nondim_specs) + + @staticmethod + def normalize(dim_specs): + nondim_specs = deepcopy(dim_specs) + for effector in effectors: + normalize(dim_specs, nondim_specs, effector, ref_effector) + + @staticmethod + def compute_energy(eptm, full_output=False): + energies = [f.energy(eptm) for f in effectors] + norm_factor = eptm.specs["settings"].get("nrj_norm_factor", 1) + if full_output: + return [E / norm_factor for E in energies] + + return sum(E.sum() for E in energies) / norm_factor + + @staticmethod + def compute_gradient(eptm, components=False): + norm_factor = eptm.specs["settings"].get("nrj_norm_factor", 1) + grads = [f.gradient(eptm) for f in effectors] + if components: + return grads + + grad_s, grad_t, grad_v = None, None, None + + srce_grads = [g[0] for g in grads if g[0].shape[0] == eptm.Ne] + if srce_grads: + grad_s = eptm.sum_srce(sum(srce_grads)) + trgt_grads = [ + g[1] for g in grads if (g[1] is not None) and (g[1].shape[0] == eptm.Ne) + ] + if trgt_grads: + grad_t = eptm.sum_trgt(sum(trgt_grads)) + vert_grads = [g[0] for g in grads if g[0].shape[0] == eptm.Nv] + if vert_grads: + grad_v = sum(vert_grads) + + grad_i = sum([g for g in (grad_s, grad_t, grad_v) if g is not None]) + + grad_i.loc[(eptm.vert_df["boundary"]==1) & (eptm.vert_df["z"] < 1), "gz"] = 0 + + # g = grad_i[["gx", "gy", "gz"]].values + # n = eptm.vert_df[["x", "y", "z"]].values + # n_norm = np.linalg.norm(n, axis=1, keepdims=True) + # n_hat = n/n_norm + # dot_product = (g * n_hat).sum(axis=1, keepdims=True) + # g_proj_plane = g - dot_product * n_hat + # grad_i = pd.DataFrame(g_proj_plane, columns=['gx', 'gy', 'gz']) + + return grad_i / norm_factor + + return NewModel + +def model_factory_cylinder(effectors, ref_effector=None): + """Produces a Model class with the provided effectors. + + Parameters + ---------- + effectors : list of :class:`.effectors.AbstractEffectors` classes. + ref_effector : optional, default None + if passed, will be used for normalization, + by default, the last effector in the list is used + + Returns + ------- + NewModel : a Model derived class with compute_enregy and compute_gradient + methods + + """ + if ref_effector is None: + ref_effector = effectors[-1] + + class NewModel: + + labels = [] + specs = { + "cell": {}, + "face": {}, + "edge": {}, + "vert": {}, + "settings": {"nrj_norm_factor": 1.0}, + } + + _effectors = effectors + + for f in effectors: + labels.append(f.label) + try: + for k in specs: + specs[k].update(f.specs.get(k, {})) + except ValueError: + warnings.warn( + """ +Since 0.7, you need to provide a default value for each of the +specs parameters, e.g. + specs = { + "face": { + "perimeter": 1.0, + "perimeter_elasticity": 0.1, + "prefered_perimeter": 3.81, + } + } + +Setting all default values to 1.0 for now +""" + ) + for k in specs: + specs[k].update({key: 1.0 for key in f.specs.get(k, {})}) + + @staticmethod + def dimensionalize(nondim_specs): + dim_specs = deepcopy(nondim_specs) + for effector in effectors: + if effector == ref_effector: + continue + dimensionalize(nondim_specs, dim_specs, effector, ref_effector) + + ref_nrj = ref_effector.get_nrj_norm(dim_specs) + dim_specs["settings"]["nrj_norm_factor"] = ref_nrj + return dim_specs + + @classmethod + def dimentionalize(cls, nondim_specs): + warnings.warn( + """This badly worded method is deprecated, + use dimensionalize instead""" + ) + return cls.dimensionalize(nondim_specs) + + @staticmethod + def normalize(dim_specs): + nondim_specs = deepcopy(dim_specs) + for effector in effectors: + normalize(dim_specs, nondim_specs, effector, ref_effector) + + @staticmethod + def compute_energy(eptm, full_output=False): + energies = [f.energy(eptm) for f in effectors] + norm_factor = eptm.specs["settings"].get("nrj_norm_factor", 1) + if full_output: + return [E / norm_factor for E in energies] + + return sum(E.sum() for E in energies) / norm_factor + + @staticmethod + def compute_gradient(eptm, components=False): + norm_factor = eptm.specs["settings"].get("nrj_norm_factor", 1) + grads = [f.gradient(eptm) for f in effectors] + if components: + return grads + + grad_s, grad_t, grad_v = None, None, None + + srce_grads = [g[0] for g in grads if g[0].shape[0] == eptm.Ne] + if srce_grads: + grad_s = eptm.sum_srce(sum(srce_grads)) + trgt_grads = [ + g[1] for g in grads if (g[1] is not None) and (g[1].shape[0] == eptm.Ne) + ] + if trgt_grads: + grad_t = eptm.sum_trgt(sum(trgt_grads)) + vert_grads = [g[0] for g in grads if g[0].shape[0] == eptm.Nv] + if vert_grads: + grad_v = sum(vert_grads) + + grad_i = sum([g for g in (grad_s, grad_t, grad_v) if g is not None]) + + grad_i.loc[(eptm.vert_df["boundary"] == 1) & (eptm.vert_df["z"] < 0), ["gx", "gy", "gz"]] = 0 + grad_i.loc[(eptm.vert_df["boundary"] == 1) & (eptm.vert_df["z"] > 0), ["gx", "gy"]] = 0 + + return grad_i / norm_factor + + return NewModel + +def model_factory_bound(effectors, ref_effector=None): + """Produces a Model class with the provided effectors. + + Parameters + ---------- + effectors : list of :class:`.effectors.AbstractEffectors` classes. + ref_effector : optional, default None + if passed, will be used for normalization, + by default, the last effector in the list is used + + Returns + ------- + NewModel : a Model derived class with compute_enregy and compute_gradient + methods + + """ + if ref_effector is None: + ref_effector = effectors[-1] + + class NewModel: + + labels = [] + specs = { + "cell": {}, + "face": {}, + "edge": {}, + "vert": {}, + "settings": {"nrj_norm_factor": 1.0}, + } + + _effectors = effectors + + for f in effectors: + labels.append(f.label) + try: + for k in specs: + specs[k].update(f.specs.get(k, {})) + except ValueError: + warnings.warn( + """ +Since 0.7, you need to provide a default value for each of the +specs parameters, e.g. + specs = { + "face": { + "perimeter": 1.0, + "perimeter_elasticity": 0.1, + "prefered_perimeter": 3.81, + } + } + +Setting all default values to 1.0 for now +""" + ) + for k in specs: + specs[k].update({key: 1.0 for key in f.specs.get(k, {})}) + + @staticmethod + def dimensionalize(nondim_specs): + dim_specs = deepcopy(nondim_specs) + for effector in effectors: + if effector == ref_effector: + continue + dimensionalize(nondim_specs, dim_specs, effector, ref_effector) + + ref_nrj = ref_effector.get_nrj_norm(dim_specs) + dim_specs["settings"]["nrj_norm_factor"] = ref_nrj + return dim_specs + + @classmethod + def dimentionalize(cls, nondim_specs): + warnings.warn( + """This badly worded method is deprecated, + use dimensionalize instead""" + ) + return cls.dimensionalize(nondim_specs) + + @staticmethod + def normalize(dim_specs): + nondim_specs = deepcopy(dim_specs) + for effector in effectors: + normalize(dim_specs, nondim_specs, effector, ref_effector) + + @staticmethod + def compute_energy(eptm, full_output=False): + energies = [f.energy(eptm) for f in effectors] + norm_factor = eptm.specs["settings"].get("nrj_norm_factor", 1) + if full_output: + return [E / norm_factor for E in energies] + + return sum(E.sum() for E in energies) / norm_factor + + @staticmethod + def compute_gradient(eptm, components=False): + norm_factor = eptm.specs["settings"].get("nrj_norm_factor", 1) + grads = [f.gradient(eptm) for f in effectors] + if components: + return grads + + grad_s, grad_t, grad_v = None, None, None + + srce_grads = [g[0] for g in grads if g[0].shape[0] == eptm.Ne] + if srce_grads: + grad_s = eptm.sum_srce(sum(srce_grads)) + trgt_grads = [ + g[1] for g in grads if (g[1] is not None) and (g[1].shape[0] == eptm.Ne) + ] + if trgt_grads: + grad_t = eptm.sum_trgt(sum(trgt_grads)) + vert_grads = [g[0] for g in grads if g[0].shape[0] == eptm.Nv] + if vert_grads: + grad_v = sum(vert_grads) + + grad_i = sum([g for g in (grad_s, grad_t, grad_v) if g is not None]) + + mask = eptm.vert_df["boundary"].values + grad_i.loc[mask == 1, ["gx", "gy", "gz"]] = 0 + + return grad_i / norm_factor + + return NewModel \ No newline at end of file diff --git a/src/tyssue/dynamics/sheet_gradients.py b/src/tyssue/dynamics/sheet_gradients.py index 58115ab1..1bb924d3 100644 --- a/src/tyssue/dynamics/sheet_gradients.py +++ b/src/tyssue/dynamics/sheet_gradients.py @@ -33,10 +33,22 @@ def area_grad(sheet): coords = sheet.coords ncoords = sheet.ncoords - inv_area = sheet.edge_df.eval("1 / (4 * sub_area)") - # Some segmentations create null areas - inv_area.replace(np.inf, 0, inplace=True) - inv_area.replace(-np.inf, 0, inplace=True) + # inv_area = sheet.edge_df.eval("1 / (4 * sub_area)") + # # Some segmentations create null areas + # inv_area.replace(np.inf, 0, inplace=True) + # inv_area.replace(-np.inf, 0, inplace=True) + + sub_area = sheet.edge_df["sub_area"].to_numpy() + + inv_area = pd.Series( + np.divide( + 1.0, + 4.0 * sub_area, + out=np.zeros_like(sub_area, dtype=float), + where=sub_area != 0, + ), + index=sheet.edge_df.index, + ) face_pos = sheet.edge_df[["f" + c for c in coords]].values srce_pos = sheet.edge_df[["s" + c for c in coords]].values diff --git a/src/tyssue/generation/shapes.py b/src/tyssue/generation/shapes.py index 98dea48b..c050beaa 100644 --- a/src/tyssue/generation/shapes.py +++ b/src/tyssue/generation/shapes.py @@ -21,7 +21,7 @@ from ..topology import type1_transition from .from_voronoi import from_3d_voronoi -from .._mesh_generation import make_spherical +# from .._mesh_generation import make_spherical # except ImportError: # print( @@ -380,43 +380,43 @@ def ellipsoid_sheet(a, b, c, n_zs, **kwargs): return eptm -def spherical_sheet(radius, Nf, Lloyd_relax=False, **kwargs): - """Returns a spherical sheet with the given radius and (approximately) - the given number of cells - """ - - centers = np.array(make_spherical(Nf)) - eptm = sheet_from_cell_centers(centers, **kwargs) - - rhos = (eptm.vert_df[eptm.coords] ** 2).sum(axis=1).mean() - ClosedSheetGeometry.scale(eptm, radius / rhos, eptm.coords) - - ClosedSheetGeometry.update_all(eptm) - if Lloyd_relax: - eptm = Lloyd_relaxation( - eptm, ClosedSheetGeometry, steps=100, update_method=update_on_sphere - ) - - return eptm - - -def spherical_monolayer(R_in, R_out, Nc, apical="out", Lloyd_relax=False): - """Returns a spherical monolayer with the given inner and - outer radii, and approximately the gieven number of cells. - - The `apical` argument can be 'in' out 'out' to specify wether - the apical face of the cells faces inward or outward, reespectively. - """ - sheet = spherical_sheet(R_in, Nc, Lloyd_relax=Lloyd_relax) - delta_R = R_out - R_in - mono = Monolayer("mono", extrude(sheet.datasets, method="normals", scale=-delta_R)) - if apical == "out": - swap_apico_basal(mono) - else: - mono.settings["lumen_side"] = "apical" - - ClosedMonolayerGeometry.update_all(mono) - return mono +# def spherical_sheet(radius, Nf, Lloyd_relax=False, **kwargs): +# """Returns a spherical sheet with the given radius and (approximately) +# the given number of cells +# """ +# +# centers = np.array(make_spherical(Nf)) +# eptm = sheet_from_cell_centers(centers, **kwargs) +# +# rhos = (eptm.vert_df[eptm.coords] ** 2).sum(axis=1).mean() +# ClosedSheetGeometry.scale(eptm, radius / rhos, eptm.coords) +# +# ClosedSheetGeometry.update_all(eptm) +# if Lloyd_relax: +# eptm = Lloyd_relaxation( +# eptm, ClosedSheetGeometry, steps=100, update_method=update_on_sphere +# ) +# +# return eptm + + +# def spherical_monolayer(R_in, R_out, Nc, apical="out", Lloyd_relax=False): +# """Returns a spherical monolayer with the given inner and +# outer radii, and approximately the gieven number of cells. +# +# The `apical` argument can be 'in' out 'out' to specify wether +# the apical face of the cells faces inward or outward, reespectively. +# """ +# sheet = spherical_sheet(R_in, Nc, Lloyd_relax=Lloyd_relax) +# delta_R = R_out - R_in +# mono = Monolayer("mono", extrude(sheet.datasets, method="normals", scale=-delta_R)) +# if apical == "out": +# swap_apico_basal(mono) +# else: +# mono.settings["lumen_side"] = "apical" +# +# ClosedMonolayerGeometry.update_all(mono) +# return mono def sheet_from_cell_centers(points, noise=0, interp_s=1e-4): diff --git a/src/tyssue/geometry/base_geometry.py b/src/tyssue/geometry/base_geometry.py index a6ef791d..6b957cd0 100644 --- a/src/tyssue/geometry/base_geometry.py +++ b/src/tyssue/geometry/base_geometry.py @@ -4,7 +4,6 @@ class BaseGeometry: - """ """ @staticmethod def update_all(sheet): diff --git a/src/tyssue/geometry/cylinder_geometry.py b/src/tyssue/geometry/cylinder_geometry.py new file mode 100644 index 00000000..4c3119a5 --- /dev/null +++ b/src/tyssue/geometry/cylinder_geometry.py @@ -0,0 +1,215 @@ +from cmath import sqrt +import numpy as np +import pandas as pd +import math + +from .sheet_geometry import SheetGeometry + +class CylinderGeometryInit(SheetGeometry): + + @staticmethod + def update_boundary_index(sheet): + + sheet.vert_df['boundary'] = 0 + sheet.edge_df['boundary'] = 0 + + sheet.get_opposite() + + sheet.edge_df.loc[sheet.edge_df['opposite'] == -1, 'boundary'] = 1 + boundary_verts = sheet.edge_df.loc[sheet.edge_df['opposite'] == -1, 'trgt'].to_numpy() + + sheet.vert_df.loc[boundary_verts, "boundary"] = 1 + + @staticmethod + def update_tangents(sheet): + + vert_coords = sheet.vert_df[sheet.coords] + vert_coords.loc[:, "z"] = 0 + vert_coords = vert_coords.values + normal = np.column_stack((np.zeros(sheet.Nv), np.zeros(sheet.Nv), np.ones(sheet.Nv))) + + tangent = np.cross(vert_coords, normal) + tangent = pd.DataFrame(tangent) + + tangent.columns = ["t" + u for u in sheet.coords] + + length = pd.DataFrame(tangent.eval("sqrt(tx**2 + ty**2 +tz**2)"), columns=['length']) + tangent["length"] = length["length"] + + tangent = tangent[['tx', 'ty', 'tz']].div(length.length, axis=0) + + for u in sheet.coords: + sheet.vert_df["t" + u] = tangent["t" + u] + + @staticmethod + def update_face_tangents(sheet): + + face_coords = sheet.face_df[sheet.coords] + face_coords["z"] = 0 + face_coords = sheet.face_df[sheet.coords].values + normal = np.column_stack((np.zeros(sheet.Nf), np.zeros(sheet.Nf), np.ones(sheet.Nf))) + + tangent = np.cross(face_coords, normal) + tangent = pd.DataFrame(tangent) + + tangent.columns = ["t" + u for u in sheet.coords] + + length = pd.DataFrame(tangent.eval("sqrt(tx**2 + ty**2 +tz**2)"), columns=['length']) + tangent["length"] = length["length"] + + tangent = tangent[['tx', 'ty', 'tz']].div(length.length, axis=0) + + for u in sheet.coords: + sheet.face_df["t" + u] = tangent["t" + u] + + @staticmethod + def update_face_distance(sheet): + sheet.face_df['distance_z_axis'] = sheet.face_df.eval( + "sqrt(x** 2 + y** 2)" + ) + + @staticmethod + def update_vert_distance(sheet): + sheet.vert_df['distance_z_axis'] = sheet.vert_df.eval( + "sqrt(x** 2 + y** 2)" + ) + + @staticmethod + def update_vert_deviation(sheet): + if "dev_length" not in sheet.vert_df.columns: + sheet.vert_df["dev_length"] = np.nan + if "dx" not in sheet.vert_df.columns: + sheet.vert_df["dx"] = np.nan + if "dy" not in sheet.vert_df.columns: + sheet.vert_df["dy"] = np.nan + if "dz" not in sheet.vert_df.columns: + sheet.vert_df["dz"] = np.nan + + edge_np = sheet.edge_df.to_numpy() + edge_dict = dict(zip(sheet.edge_df.columns, + list(range(0, len(sheet.edge_df.columns))))) + vert_np = sheet.vert_df.to_numpy() + vert_dict = dict(zip(sheet.vert_df.columns, + list(range(0, len(sheet.vert_df.columns))))) + + grad1 = np.nan + grad2 = np.nan + gradt = np.nan + lenth = [] + + for i in sheet.vert_df.index: + mask = (i == edge_np[:, edge_dict["srce"]]) + neighbor_verts = edge_np[mask, edge_dict["trgt"]].tolist() + neighbor_verts = vert_np[neighbor_verts][:, [vert_dict["x"], vert_dict["y"], vert_dict["z"]]] + vert_coords = vert_np[i, [vert_dict["x"], vert_dict["y"], vert_dict["z"]]] + + if len(neighbor_verts) >= 3: + center = (neighbor_verts[0] + neighbor_verts[1] + neighbor_verts[2]) / 3 + + grad = np.array(vert_coords) - np.array(center) + + length = np.linalg.norm(grad) + + grad = grad / length + + else: + grad = np.array([0, 0, 0]) + + if i == 0: + grad1 = grad + + if i == 1: + grad2 = grad + gradt = np.vstack((grad1, grad2)) + + if i >= 2: + gradt = np.vstack((gradt, grad)) + + lenth.append(length) + + gradt = pd.DataFrame(gradt, columns=['dx', 'dy', 'dz']) + sheet.vert_df[['dx', 'dy', 'dz']] = gradt + sheet.vert_df['dev_length'] = lenth + + @staticmethod + def update_vert_deviation2(sheet): + if "dev_length" not in sheet.vert_df.columns: + sheet.vert_df["dev_length"] = np.nan + if "dx" not in sheet.vert_df.columns: + sheet.vert_df["dx"] = np.nan + if "dy" not in sheet.vert_df.columns: + sheet.vert_df["dy"] = np.nan + if "dz" not in sheet.vert_df.columns: + sheet.vert_df["dz"] = np.nan + + for i in sheet.vert_df.index: + vert = sheet.vert_df.loc[i, ["x", "y", "z"]].to_numpy() + neighbors = sheet.edge_df.loc[sheet.edge_df["srce"] == 6, ["tx", "ty", "tz"]].to_numpy() + center = neighbors.sum(axis=0) / len(neighbors) + grad = vert - center + length = np.linalg.norm(grad) + sheet.vert_df[['dx', 'dy', 'dz']] = grad + sheet.vert_df['dev_length'] = length + + @staticmethod + def update_lumen_vol(sheet): + lumen_pos_faces = sheet.edge_df[["f" + c for c in sheet.coords]].to_numpy() + lumen_sub_vol = ( + np.sum((lumen_pos_faces) * sheet.edge_df[sheet.ncoords].to_numpy(), axis=1) + / 6 + ) + lumen_volume_gross = sum(lumen_sub_vol) + + top_verts = sheet.vert_df.loc[(sheet.vert_df["boundary"] == 1) & (sheet.vert_df["z"] > 0)] + top_radius = top_verts["distance_z_axis"].values.mean() + top_height = top_verts["z"].values.mean() + top_volume = (1 / 3) * math.pi * top_radius ** 2 * top_height + + bot_verts = sheet.vert_df.loc[(sheet.vert_df["boundary"] == 1) & (sheet.vert_df["z"] < 0)] + bot_radius = bot_verts["distance_z_axis"].values.mean() + bot_height = -bot_verts["z"].values.mean() + bot_volume = (1 / 3) * math.pi * bot_radius ** 2 * bot_height + + sheet.settings["lumen_vol"] = top_volume + bot_volume + lumen_volume_gross + + @staticmethod + def update_vol_cell(sheet): + sheet.settings["vol_cell"] = sheet.settings["lumen_vol"]/len(sheet.face_df) + + @staticmethod + def update_boundary_radius(sheet): + sheet.vert_df[["cx", "cy", "cz"]] = np.nan + sheet.vert_df.loc[(sheet.vert_df["boundary"] == 1) & (sheet.vert_df["z"] <= 0), ["cx", "cy", "cz"]] = ( + sheet.vert_df.loc[(sheet.vert_df["boundary"] == 1) & (sheet.vert_df["z"] <= 0), ["x", "y", "z"]] - + sheet.settings["bot_center"]).to_numpy() + sheet.vert_df.loc[(sheet.vert_df["boundary"] == 1) & (sheet.vert_df["z"] >= 0), ["cx", "cy", "cz"]] = ( + sheet.vert_df.loc[(sheet.vert_df["boundary"] == 1) & (sheet.vert_df["z"] >= 0), ["x", "y", "z"]] - + sheet.settings["top_center"]).to_numpy() + sheet.vert_df["bound_rad"] = sheet.vert_df.eval("(cx**2 + cy**2 + cz**2) ** 0.5") + sheet.vert_df[["cx", "cy", "cz"]] = sheet.vert_df[["cx", "cy", "cz"]].div(sheet.vert_df["bound_rad"], axis=0) + sheet.vert_df.fillna(0, inplace=True) + + @classmethod + def update_all(cls, sheet): + super().update_all(sheet) + cls.update_boundary_index(sheet) + cls.update_tangents(sheet) + # cls.update_face_tangents(sheet) + # cls.update_face_distance(sheet) + cls.update_vert_distance(sheet) + cls.update_vert_deviation(sheet) + cls.update_lumen_vol(sheet) + # cls.update_vol_cell(sheet) + # cls.update_boundary_radius(sheet) + + +class CylinderGeometry(CylinderGeometryInit): + + @classmethod + def update_all(cls, sheet): + super().update_all(sheet) + # cls.update_preflumen_volume(sheet) + + @staticmethod + def update_preflumen_volume(sheet): + sheet.settings["lumen_prefered_vol"] = sheet.settings["vol_cell"] * len(sheet.face_df) \ No newline at end of file diff --git a/src/tyssue/geometry/sheet_geometry.py b/src/tyssue/geometry/sheet_geometry.py index c7201635..d7a218ae 100644 --- a/src/tyssue/geometry/sheet_geometry.py +++ b/src/tyssue/geometry/sheet_geometry.py @@ -1,8 +1,10 @@ +from cmath import sqrt import numpy as np import pandas as pd +import math from .planar_geometry import PlanarGeometry -from .utils import rotation_matrices, rotation_matrix +from .utils import rotation_matrix, rotation_matrices class SheetGeometry(PlanarGeometry): @@ -25,11 +27,12 @@ def update_all(cls, sheet): cls.update_ucoords(sheet) cls.update_length(sheet) cls.update_centroid(sheet) - cls.update_height(sheet) + # cls.update_height(sheet) cls.update_normals(sheet) cls.update_areas(sheet) cls.update_perimeters(sheet) cls.update_vol(sheet) + cls.update_boundary_index(sheet) @staticmethod def update_normals(sheet): @@ -49,7 +52,7 @@ def update_areas(sheet): Updates the normal coordniate of each (srce, trgt, face) face. """ sheet.edge_df["sub_area"] = ( - np.linalg.norm(sheet.edge_df[sheet.ncoords], axis=1) / 2 + np.linalg.norm(sheet.edge_df[sheet.ncoords], axis=1) / 2 ) sheet.face_df["area"] = sheet.sum_face(sheet.edge_df["sub_area"]) @@ -61,7 +64,7 @@ def update_vol(sheet): """ sheet.edge_df["sub_vol"] = ( - sheet.upcast_srce(sheet.vert_df["height"]) * sheet.edge_df["sub_area"] + sheet.upcast_srce(sheet.vert_df["height"]) * sheet.edge_df["sub_area"] ) sheet.face_df["vol"] = sheet.sum_face(sheet.edge_df["sub_vol"]) @@ -123,7 +126,29 @@ def update_height(cls, sheet): edge_height = sheet.upcast_srce(sheet.vert_df[["height", "rho"]]) edge_height.set_index(sheet.edge_df["face"], append=True, inplace=True) - sheet.face_df[["height", "rho"]] = edge_height.groupby(level="face").mean() + sheet.face_df[["height", "rho"]] = edge_height.mean(level="face") + + @staticmethod + def update_boundary_index(sheet): + # Reset boundary flags + sheet.vert_df['boundary'] = 0 + sheet.edge_df['boundary'] = 0 + sheet.face_df['boundary'] = 0 + + # Update opposite edges + sheet.get_opposite() + + # Identify boundary edges + boundary_edges = sheet.edge_df['opposite'] == -1 + sheet.edge_df.loc[boundary_edges, 'boundary'] = 1 + + # Set boundary vertices + boundary_verts = sheet.edge_df.loc[boundary_edges, 'trgt'] + sheet.vert_df.loc[boundary_verts.unique(), 'boundary'] = 1 + + # Set boundary faces + boundary_faces = sheet.edge_df.loc[boundary_edges, 'face'] + sheet.face_df.loc[boundary_faces.dropna().unique().astype(int), 'boundary'] = 1 @classmethod def reset_scafold(cls, sheet): @@ -193,16 +218,17 @@ def face_projected_pos(sheet, face, psi=0): The rotated, relative positions of the face's vertices """ - rel_pos = sheet.edge_df.query(f"face == {face}")[["srce", "rx", "ry", "rz"]] - rel_pos = rel_pos.set_index("srce") - rel_pos.index.name = "vert" - _, _, rotation = np.linalg.svd( - rel_pos.to_numpy().astype(float), full_matrices=False + face_orbit = sheet.edge_df[sheet.edge_df["face"] == face]["srce"] + rel_pos = ( + sheet.vert_df.loc[face_orbit.to_numpy(), sheet.coords].to_numpy() + - sheet.face_df.loc[face, sheet.coords].to_numpy(dtype=float) ) + _, _, rotation = np.linalg.svd( + rel_pos, full_matrices=False) if psi: rotation = np.dot(rotation_matrix(psi, [0, 0, 1]), rotation) rot_pos = pd.DataFrame( - np.dot(rel_pos, rotation.T), index=rel_pos.index, columns=sheet.coords + np.dot(rel_pos, rotation.T), index=face_orbit, columns=sheet.coords ) return rot_pos @@ -213,8 +239,8 @@ def face_rotations(cls, sheet, method="normal", output_as="edge"): vertices are mostly in the u, v plane. If method is 'normal', face is oriented with it's normal along w - if method is 'svd', the u, v, w is determined through - singular value decompostion of the face vertices relative positions. + if method is 'svd', the u, v, w is determined through singular value decompostion + of the face vertices relative positions. svd is slower but more effective at reducing face dimensionality. @@ -297,6 +323,24 @@ def get_phis(cls, sheet, method="normal"): rotated = np.einsum("ikj, ik -> ij", rots, rel_srce_pos) return np.arctan2(rotated[:, 1], rotated[:, 0]) + @staticmethod + def update_boundary_index(sheet): + """Updates the vert_df and edge_df dataframes with a 'boundary' column + that takes values 0 and 1 with 1 denoting that an edge or vertex lies + on the tissue boundary, and 0 when it does not. + + """ + + sheet.vert_df['boundary'] = 0 + sheet.edge_df['boundary'] = 0 + + sheet.get_opposite() + + sheet.edge_df.loc[sheet.edge_df['opposite'] == -1, 'boundary'] = 1 + boundary_verts = sheet.edge_df.loc[sheet.edge_df['opposite'] == -1, 'trgt'].to_numpy() + + sheet.vert_df.loc[boundary_verts, "boundary"] = 1 + class ClosedSheetGeometry(SheetGeometry): """Geometry for a closed 2.5D sheet. @@ -314,54 +358,24 @@ def update_all(cls, sheet): def update_lumen_vol(sheet): lumen_pos_faces = sheet.edge_df[["f" + c for c in sheet.coords]].to_numpy() lumen_sub_vol = ( - np.sum((lumen_pos_faces) * sheet.edge_df[sheet.ncoords].to_numpy(), axis=1) - / 6 + np.sum((lumen_pos_faces) * sheet.edge_df[sheet.ncoords].to_numpy(), axis=1) + / 6 ) sheet.settings["lumen_vol"] = sum(lumen_sub_vol) -class MidlineBoundaryGeometry(ClosedSheetGeometry): - @classmethod - def update_all(cls, eptm): - super().update_all(eptm) - cls.update_delta_boundary(eptm) - - @staticmethod - def update_delta_boundary(eptm): - midline_boudary_stiffness = eptm.settings.get( - "midline_boundary_stiffness", False - ) - # update boundary transgression - # leftright = 1|-1 depending on x position at start - # x / abs(x) = 1|-1 depending on current x position - # (x/abs(x)) - leftright = 0 if both are equal (vert has not crossed midline) - # = -2 if both are unequal and current x is negative - # = 2 if both are unequal and current x is positive - # hence, take (0|2|-2)*0.5x to get the distance from x axis as a positive number - - if midline_boudary_stiffness is not False: - if "leftright" not in eptm.vert_df.columns: - eptm.vert_df["leftright"] = np.sign(eptm.vert_df["x"]) - eptm.vert_df["delta_boundary"] = ( - (np.sign(eptm.vert_df["x"]) - eptm.vert_df["leftright"]) - * eptm.vert_df["x"] - / 2 - ) - - class EllipsoidGeometry(ClosedSheetGeometry): @staticmethod def update_height(eptm): - a, b, c = eptm.settings["abc"] eptm.vert_df["theta"] = np.arcsin((eptm.vert_df.z / c).clip(-1, 1)) eptm.vert_df["vitelline_rho"] = a * np.cos(eptm.vert_df["theta"]) eptm.vert_df["basal_shift"] = ( - eptm.vert_df["vitelline_rho"] - eptm.specs["vert"]["basal_shift"] + eptm.vert_df["vitelline_rho"] - eptm.specs["vert"]["basal_shift"] ) eptm.vert_df["delta_rho"] = ( - np.linalg.norm(eptm.vert_df[["x", "y"]], axis=1) - - eptm.vert_df["vitelline_rho"] + np.linalg.norm(eptm.vert_df[["x", "y"]], axis=1) + - eptm.vert_df["vitelline_rho"] ).clip(lower=0) SheetGeometry.update_height(eptm) @@ -371,20 +385,15 @@ def scale(eptm, scale, coords): SheetGeometry.scale(eptm, scale, coords) eptm.settings["abc"] = [u * scale for u in eptm.settings["abc"]] - class WeightedPerimeterEllipsoidLameGeometry(ClosedSheetGeometry): """ EllipsoidLameGeometry correspond to a super-egg geometry with a calculation of perimeter is based on weight of each junction. Meaning if all junction of a cell have the same weight, perimeter is - calculated as a usual perimeter calculation - .. math:: - p = \\sum l_{ij} + calculated as a usual perimeter calculation (p = l_ij + l_jk + l_km + l_mn + l_ni) Otherwise, weight parameter allowed more or less importance of a junction in the - perimeter calculation - .. math:: - p = \\sum w_{ij} \\, l_{ij} + perimeter calculation (p = w_ij*l_ij + w_jk*l_jk + w_km*l_km + w_mn*l_mn + w_ni*l_ni) In this geometry, a sphere surrounding the tissue, meaning a force is apply only at the extremity of the tissue; `eptm.vert_df['delta_rho']` is computed as the @@ -422,7 +431,6 @@ def normalize_weights(sheet): @staticmethod def update_height(eptm): - eptm.vert_df["rho"] = np.linalg.norm(eptm.vert_df[eptm.coords], axis=1) r = eptm.settings["barrier_radius"] eptm.vert_df["delta_rho"] = (eptm.vert_df["rho"] - r).clip(0) @@ -430,7 +438,6 @@ def update_height(eptm): def face_svd_(faces): - rel_pos = faces[["rx", "ry", "rz"]] _, _, rotation = np.linalg.svd(rel_pos.astype(float), full_matrices=False) - return rotation + return rotation \ No newline at end of file diff --git a/src/tyssue/geometry/vessel_geometry.py b/src/tyssue/geometry/vessel_geometry.py new file mode 100644 index 00000000..120205c3 --- /dev/null +++ b/src/tyssue/geometry/vessel_geometry.py @@ -0,0 +1,58 @@ +from cmath import sqrt +import numpy as np +import pandas as pd +import math + +from .sheet_geometry import SheetGeometry +from .utils import rotation_matrix, rotation_matrices + + +class VesselGeometry(SheetGeometry): + + @staticmethod + def update_tangents(sheet): + # Extract coordinates as NumPy arrays + x = sheet.vert_df['x'].to_numpy() + y = sheet.vert_df['y'].to_numpy() + + # Analytical Cross Product + tx = y + ty = -x + + # Calculate Length (Vectorized) + length = np.hypot(tx, ty) + + # Normalize + with np.errstate(divide='ignore', invalid='ignore'): + inv_length = 1.0 / length + + # 5. Assign directly to DataFrame + # This assumes sheet.coords = ['x', 'y', 'z'], matching the original output naming + sheet.vert_df['tx'] = tx * inv_length + sheet.vert_df['ty'] = ty * inv_length + sheet.vert_df['tz'] = 0.0 # Z-component is always 0 in this projection + + @staticmethod + def update_vert_distance(sheet): + coords = ["x", "y", "z"] + if "axis" in sheet.settings.keys(): + axis = sheet.settings["axis"] + else: + axis = "z" + coords.remove(axis) + distances = np.sqrt(sheet.vert_df[coords[0]].to_numpy()**2 + sheet.vert_df[coords[1]].to_numpy()**2) + sheet.vert_df['distance_origin'] = distances + sheet.vert_df["o"+f"{coords[0]}"] = sheet.vert_df[f"{coords[0]}"]/distances + sheet.vert_df["o"+f"{coords[1]}"] = sheet.vert_df[f"{coords[1]}"]/distances + + @classmethod + def update_all(cls, sheet): + super().update_all(sheet) + cls.update_tangents(sheet) + cls.update_vert_distance(sheet) + +def face_svd_(faces): + + rel_pos = faces[["rx", "ry", "rz"]] + _, _, rotation = np.linalg.svd(rel_pos.astype(float), full_matrices=False) + return rotation diff --git a/src/tyssue/solvers/viscous.py b/src/tyssue/solvers/viscous.py index fcb4d351..9773385e 100644 --- a/src/tyssue/solvers/viscous.py +++ b/src/tyssue/solvers/viscous.py @@ -3,13 +3,20 @@ """ import logging +import numpy as np +import pandas as pd import warnings +import random + +from itertools import count +from scipy.integrate import solve_ivp -import numpy as np +from ..core.history import History from ..behaviors.event_manager import EventManager from ..behaviors.sheet.basic_events import reconnect -from ..core.history import History +from ..topology.sheet_topology import cell_division + log = logging.getLogger(__name__) MAX_ITER = 1000 @@ -27,7 +34,11 @@ def set_pos(eptm, geom, pos): class EulerSolver: - """Explicit Euler solver""" + """Explicit Euler solver + + + + """ def __init__( self, @@ -79,7 +90,7 @@ def __init__( if auto_reconnect: if manager is None: manager = EventManager() - if "reconnect" not in [n[0].__name__ for n in manager.next]: + if not "reconnect" in [n[0].__name__ for n in manager.next]: manager.append(reconnect) self.manager = manager @@ -92,7 +103,8 @@ def current_pos(self): ].values.ravel() def set_pos(self, pos): - """Updates the eptm vertices position""" + """Updates the eptm vertices position + """ return self._set_pos(self.eptm, self.geom, pos) def record(self, t): @@ -132,6 +144,9 @@ def solve(self, tf, dt, on_topo_change=None, topo_change_args=()): self.eptm.topo_changed = False self.record(t) + if t == tf: + self.history.update_datasets() + def ode_func(self, t, pos): """Computes the models' gradient. @@ -139,11 +154,6 @@ def ode_func(self, t, pos): Returns ------- dot_r : 1D np.ndarray of shape (self.eptm.Nv * self.eptm.dim, ) - - .. math:: - - \frac{dr_i}{dt} = -\frac{\nabla U_i}{\eta_i} - """ grad_U = self.model.compute_gradient(self.eptm).loc[self.eptm.active_verts] @@ -152,6 +162,108 @@ def ode_func(self, t, pos): / self.eptm.vert_df.loc[self.eptm.active_verts, "viscosity"].values[:, None] ).ravel() +class EulerSolverDivision(EulerSolver): + + def solve(self, tf, dt, torque, on_topo_change=None, topo_change_args=()): + """Solves the system of differential equations from the current time + to tf with steps of dt with a forward Euler method. + + Parameters + ---------- + tf : float, final time when we stop solving + dt : float, time step + on_topo_change : function, optional, default None + function of `self.eptm` + topo_change_args : tuple, arguments passed to `on_topo_change` + + """ + self.eptm.settings["dt"] = dt + + for t in np.arange(self.prev_t, tf + dt, dt): + + pos = self.current_pos + dot_r = self.ode_func(t, pos) + if self.bounds is not None: + dot_r = np.clip(dot_r, *self.bounds) + pos = pos + dot_r * dt + self.set_pos(pos) + self.prev_t = t + + if t%1==0 or t%0.5==0: + try: + mother = random.choice(self.eptm.face_df.loc[(self.eptm.face_df["ventral"] == 1) & (self.eptm.face_df['area'] >= 0.7)].index) + mother_torque = self.eptm.face_df.loc[mother, "torque_coef"] + daughter = cell_division(self.eptm, mother, self.geom, angle = np.pi/2) + self.eptm.face_df.loc[self.eptm.face_df.index == daughter, "torque_coef"] = mother_torque + self.eptm.face_df.loc[self.eptm.face_df.index == daughter, "ventral"] = 1 + + except: + pass + + if self.manager is not None: + self.manager.execute(self.eptm) + self.geom.update_all(self.eptm) + self.manager.update() + + if self.eptm.topo_changed: + log.info("Topology changed") + if on_topo_change is not None: + on_topo_change(*topo_change_args) + self.eptm.topo_changed = False + self.record(t) + + if t == tf: + self.history.update_datasets() + +class EulerSolverDiv(EulerSolver): + + def solve(self, tf, dt, torque, on_topo_change=None, topo_change_args=()): + """Solves the system of differential equations from the current time + to tf with steps of dt with a forward Euler method. + + Parameters + ---------- + tf : float, final time when we stop solving + dt : float, time step + on_topo_change : function, optional, default None + function of `self.eptm` + topo_change_args : tuple, arguments passed to `on_topo_change` + + """ + self.eptm.settings["dt"] = dt + + for t in np.arange(self.prev_t, tf + dt, dt): + + pos = self.current_pos + dot_r = self.ode_func(t, pos) + if self.bounds is not None: + dot_r = np.clip(dot_r, *self.bounds) + pos = pos + dot_r * dt + self.set_pos(pos) + self.prev_t = t + + if t%1==0 or t%0.5==0: + try: + mother = random.choice(self.eptm.face_df.loc[(self.eptm.face_df["torque_coef"] == torque) & (self.eptm.face_df['area'] >= 0.7)].index) + daughter = cell_division(self.eptm, mother, self.geom, angle = np.pi/2) + self.eptm.face_df.loc[self.eptm.face_df.index == daughter, "torque_coef"] = torque + except: + pass + + if self.manager is not None: + self.manager.execute(self.eptm) + self.geom.update_all(self.eptm) + self.manager.update() + + if self.eptm.topo_changed: + log.info("Topology changed") + if on_topo_change is not None: + on_topo_change(*topo_change_args) + self.eptm.topo_changed = False + self.record(t) + + if t == tf: + self.history.update_datasets() class IVPSolver: def __init__(self, *args, **kwargs): diff --git a/src/tyssue/topology/__init__.py b/src/tyssue/topology/__init__.py index 172bd0a1..1a61d5b3 100644 --- a/src/tyssue/topology/__init__.py +++ b/src/tyssue/topology/__init__.py @@ -18,12 +18,17 @@ merge_vertices, split_vert, ) -from .bulk_topology import ( +from .bulk_topology import ( # noqa: F401 HI_transition, IH_transition, + all_lateral_fusions, find_HIs, find_IHs, + find_fusion_nucleations, + find_fusion_propagations, + find_fusion_splits, find_rearangements, + fuse_lateral_faces, ) from .sheet_topology import remove_face, type1_transition diff --git a/src/tyssue/topology/base_topology.py b/src/tyssue/topology/base_topology.py index 920f527d..13b0e2a7 100644 --- a/src/tyssue/topology/base_topology.py +++ b/src/tyssue/topology/base_topology.py @@ -193,7 +193,7 @@ def remove_face(sheet, face): verts = edges["srce"].unique() new_vert_data = sheet.vert_df.loc[verts[0] : verts[0]].copy() - new_vert_data[sheet.coords] = sheet.vert_df.loc[verts, sheet.coords].mean() + new_vert_data[sheet.coords] = sheet.vert_df.loc[verts, sheet.coords].mean().to_numpy() sheet.vert_df = pd.concat( [sheet.vert_df, pd.DataFrame(new_vert_data)], ignore_index=True ) @@ -223,6 +223,8 @@ def remove_face(sheet, face): sheet.reset_index() sheet.reset_topo() + sheet.network_changed = True + return new_vert diff --git a/src/tyssue/topology/bulk_topology.py b/src/tyssue/topology/bulk_topology.py index 43459b54..3cac58b0 100644 --- a/src/tyssue/topology/bulk_topology.py +++ b/src/tyssue/topology/bulk_topology.py @@ -4,6 +4,9 @@ import numpy as np import pandas as pd +from scipy.spatial import cKDTree +from scipy.spatial.distance import cdist +from scipy.optimize import linear_sum_assignment from ..core.monolayer import Monolayer from ..core.objects import _is_closed_cell, euler_characteristic @@ -551,3 +554,422 @@ def fix_pinch(eptm): eptm.face_df = eptm.face_df.drop(bad_faces) eptm.reset_index() eptm.reset_topo() + + +def _face_centroids_normals(eptm, faces): + """Returns the centroid and unit (Newell) normal of each face in `faces`. + + Computed directly from vertex positions, so no geometry update is required + and the result is independent of the half-edge order within a face. + """ + centroids, normals = {}, {} + grouped = eptm.edge_df[eptm.edge_df["face"].isin(faces)].groupby("face") + for face, edges in grouped: + srce_pos = eptm.vert_df.loc[edges["srce"], eptm.coords].to_numpy() + trgt_pos = eptm.vert_df.loc[edges["trgt"], eptm.coords].to_numpy() + centroids[face] = srce_pos.mean(axis=0) + normal = np.cross(srce_pos, trgt_pos).sum(axis=0) + norm = np.linalg.norm(normal) + normals[face] = normal / norm if norm > 0 else normal + return centroids, normals + + +def _ordered_boundary(eptm, face): + """Vertex indices around `face`, in half-edge order.""" + edges = eptm.edge_df[eptm.edge_df["face"] == face] + nxt = dict(zip(edges["srce"].astype(int), edges["trgt"].astype(int))) + start = int(edges["srce"].iloc[0]) + order, v = [start], nxt[start] + while v != start and len(order) <= len(nxt): + order.append(v) + v = nxt[v] + return order + + +def _lateral_apico_basal(eptm, face): + """Returns (apical_verts, basal_verts) of a lateral `face`, from the apical + / basal labels carried by its vertices.""" + verts = _ordered_boundary(eptm, face) + seg = eptm.vert_df.loc[verts, "segment"] + apical = [v for v in verts if seg[v] == "apical"] + basal = [v for v in verts if seg[v] == "basal"] + return apical, basal + + +def _cell_adjacency(eptm): + """Set of frozenset({cell_a, cell_b}) for cells sharing a face.""" + eptm.get_opposite_faces() + face_cell = eptm.edge_df.groupby("face")["cell"].first() + adjacency = set() + for face, opp in eptm.face_df["opposite"].items(): + if opp != -1 and face in face_cell.index and opp in face_cell.index: + adjacency.add(frozenset((int(face_cell[face]), int(face_cell[opp])))) + return adjacency + + +def _free_lateral_faces(eptm, segment="lateral"): + """Indices of the free (border) faces, restricted to `segment` if set.""" + eptm.get_opposite_faces() + border = eptm.face_df[eptm.face_df["opposite"] == -1] + if segment is not None and "segment" in eptm.face_df.columns: + border = border[border["segment"] == segment] + return list(border.index) + + +def _split_lateral_face(eptm, face): + """Cuts a free lateral quad in two across its apico-basal axis. + + A new apical vertex is inserted on the face's apical edge and a new basal + vertex on its basal edge (each edge is split in the apical / basal cap face + too, via :func:`~tyssue.topology.base_topology.add_vert`), and the two are + joined by a new edge, dividing the wall into two coplanar lateral quads of + the same cell. This reconciles a valence mismatch in a closing seam: a single + wide wall apposed to two narrower walls of the other front (a "2-vs-1 + pocket") is split so each piece can fuse with its own partner. + + Returns ``(new_apical_vert, new_basal_vert, new_face)`` (indices valid until + the following ``reset_index``). + """ + def _segment_edge(face, kind): + seg = eptm.vert_df["segment"] + edges = eptm.edge_df[eptm.edge_df["face"] == face] + match = (seg.loc[edges["srce"]].values == kind) & ( + seg.loc[edges["trgt"]].values == kind + ) + return edges.index[match][0] + + na, _, _ = add_vert(eptm, _segment_edge(face, "apical")) + nb, _, _ = add_vert(eptm, _segment_edge(face, "basal")) + cell = int(eptm.edge_df.loc[eptm.edge_df["face"] == face, "cell"].iloc[0]) + + # cyclic boundary is [a1, na, a2, b1, nb, b2]; na and nb are antipodal, so + # the half-loop from na to nb is one of the two new quads (na, a2, b1, nb). + order = _ordered_boundary(eptm, face) + n = len(order) + ina = order.index(na) + half = [order[(ina + k) % n] for k in range(n // 2 + 1)] + assert half[-1] == nb, "unexpected lateral-face boundary in split" + + new_face_row = eptm.face_df.loc[[face]].copy() + new_face_row.index = [eptm.face_df.index.max() + 1] + eptm.face_df = pd.concat([eptm.face_df, new_face_row]) + new_face = eptm.face_df.index[-1] + + def _half_edge(srce, trgt): + edges = eptm.edge_df + return edges.index[ + (edges["face"] == face) & (edges["srce"] == srce) & (edges["trgt"] == trgt) + ][0] + + for srce, trgt in zip(half[:-1], half[1:]): + eptm.edge_df.loc[_half_edge(srce, trgt), "face"] = new_face + + template = eptm.edge_df[eptm.edge_df["face"] == face].iloc[[0]] + for srce, trgt, fc in ((nb, na, new_face), (na, nb, face)): + new_edge = template.copy() + new_edge.index = [eptm.edge_df.index.max() + 1] + eptm.edge_df = pd.concat([eptm.edge_df, new_edge]) + idx = eptm.edge_df.index[-1] + eptm.edge_df.loc[idx, ["srce", "trgt", "face", "cell"]] = [srce, trgt, fc, cell] + + eptm.reset_index() + eptm.reset_topo() + logger.info("split lateral face %d into %d and %d", face, face, new_face) + return na, nb, new_face + + +def fuse_lateral_faces(eptm, face_g, face_f, validate=True): + """Welds two apposed free lateral faces into one internal interface. + + The vertices of `face_g` (cell A) and `face_f` (cell B) are matched + **within their apico-basal segment** -- apical vertices only to apical, basal + only to basal -- and welded pairwise, so the two faces end up with an + identical vertex set and :meth:`Epithelium.get_opposite_faces` registers them + as opposites. Vertices the two faces already share (a seam corner / edge from + a neighbouring fusion) are left in place; only the remaining ones are welded. + + The weld requires equal numbers of unshared apical (and basal) vertices on + the two faces -- the quad-quad case. A genuine valence/registration mismatch + (e.g. differing rim densities) returns -1 rather than fusing. + + Returns 0 on success, -1 if it could not (validly) fuse. + """ + if validate: + vert_bak = eptm.vert_df.copy() + edge_bak = eptm.edge_df.copy() + face_bak = eptm.face_df.copy() + + g_ap, g_ba = _lateral_apico_basal(eptm, face_g) + f_ap, f_ba = _lateral_apico_basal(eptm, face_f) + + shared = set(g_ap + g_ba) & set(f_ap + f_ba) + g_ap = [v for v in g_ap if v not in shared] + g_ba = [v for v in g_ba if v not in shared] + f_ap = [v for v in f_ap if v not in shared] + f_ba = [v for v in f_ba if v not in shared] + + if not (g_ap or g_ba): + return 0 # already coincident + if len(g_ap) != len(f_ap) or len(g_ba) != len(f_ba): + return -1 # valence / registration mismatch -- deferred + + pairs = [] + for gv, fv in (g_ap, f_ap), (g_ba, f_ba): + if not gv: + continue + cost = cdist( + eptm.vert_df.loc[gv, eptm.coords].to_numpy(), + eptm.vert_df.loc[fv, eptm.coords].to_numpy(), + ) + rows, cols = linear_sum_assignment(cost) + pairs += [(gv[r], fv[c]) for r, c in zip(rows, cols)] + + for va, vb in pairs: + v_keep, v_drop = sorted((int(va), int(vb))) + eptm.vert_df.loc[v_keep, eptm.coords] = ( + eptm.vert_df.loc[[v_keep, v_drop], eptm.coords].mean(axis=0).to_numpy() + ) + eptm.edge_df.replace({"srce": v_drop, "trgt": v_drop}, v_keep, inplace=True) + eptm.vert_df.drop(v_drop, axis=0, inplace=True) + + degenerate = eptm.edge_df.query("srce == trgt") + if degenerate.shape[0]: + eptm.edge_df.drop(degenerate.index, axis=0, inplace=True) + + eptm.reset_index() + eptm.reset_topo() + fix_pinch(eptm) + eptm.get_opposite_faces() + eptm.reset_index() + eptm.reset_topo() + + if validate and not eptm.validate(): + eptm.vert_df = vert_bak + eptm.edge_df = edge_bak + eptm.face_df = face_bak + eptm.reset_index() + eptm.reset_topo() + eptm.get_opposite_faces() + logger.info("rolled back invalid lateral fusion (%d, %d)", face_g, face_f) + return -1 + + logger.info("fused lateral faces %d and %d", face_g, face_f) + return 0 + + +def find_fusion_nucleations(eptm, d_max=None, theta_face=45.0, margin=0.5, + segment="lateral"): + """Finds apposed free lateral faces that are approaching but not yet joined. + + A pair ``(g, f)`` is returned when both are free `segment` faces of two + **different, not-yet-adjacent** cells that do **not** already share an edge, + their outward normals are anti-parallel to within `theta_face`, and a vertex + of one lies within `d_max` of the other's plane (and laterally over it, to + within `margin`). This nucleates a seam where two fronts first meet, at any + relative angle -- no face coincidence required. + """ + faces = _free_lateral_faces(eptm, segment) + if len(faces) < 2: + return [] + + if d_max is None: + d_max = eptm.settings.get("fusion_distance", None) + if d_max is None: + bedges = eptm.edge_df[eptm.edge_df["face"].isin(faces)] + if "length" in bedges.columns: + d_max = float(bedges["length"].mean()) + else: + d_max = 0.5 + + centroids, normals = _face_centroids_normals(eptm, faces) + bedges = eptm.edge_df[eptm.edge_df["face"].isin(faces)] + face_cell = bedges.groupby("face")["cell"].first().to_dict() + face_verts = bedges.groupby("face")["srce"].apply(lambda s: set(s)).to_dict() + radius = { + face: np.linalg.norm( + eptm.vert_df.loc[list(face_verts[face]), eptm.coords].to_numpy() + - centroids[face], + axis=1, + ).max() + for face in faces + } + adjacency = _cell_adjacency(eptm) + + cent_arr = np.array([centroids[f] for f in faces]) + tree = cKDTree(cent_arr) + cos_face = np.cos(np.radians(theta_face)) + max_radius = max(radius.values()) + + candidates = [] + for i, j in tree.query_pairs(d_max + 2 * max_radius): + g, f = faces[i], faces[j] + cg, cf = face_cell[g], face_cell[f] + if cg == cf or frozenset((cg, cf)) in adjacency: + continue + if len(face_verts[g] & face_verts[f]) >= 2: + continue # share an edge -> handled by propagation + if np.dot(normals[g], normals[f]) > -cos_face: + continue # not anti-parallel enough + gpos = eptm.vert_df.loc[list(face_verts[g]), eptm.coords].to_numpy() + rel = gpos - centroids[f] + perp = rel @ normals[f] + lateral = np.linalg.norm(rel - np.outer(perp, normals[f]), axis=1) + hit = (np.abs(perp) < d_max) & (lateral < radius[f] * (1 + margin)) + if hit.any(): + candidates.append((g, f, np.abs(perp[hit]).min())) + + candidates.sort(key=lambda c: c[2]) + out, used = [], set() + for g, f, _ in candidates: + if g in used or f in used: + continue + out.append((int(g), int(f))) + used.update((g, f)) + return out + + +def find_fusion_propagations(eptm, theta_fold=20.0, segment="lateral"): + """Finds free lateral faces that have folded together along a shared edge. + + A pair ``(g, f)`` is returned when both are free `segment` faces of two + **different, not-yet-adjacent** cells that **share an edge** (the seam edge + welded by a neighbouring fusion) and whose *fold angle* about that edge is + below `theta_fold` -- i.e. the two flaps have closed onto each other. This + walks the seam outward from a nucleation as the rim closes. Same-front + neighbours are excluded because their cells are already face-adjacent. + """ + faces = _free_lateral_faces(eptm, segment) + if len(faces) < 2: + return [] + + centroids, _ = _face_centroids_normals(eptm, faces) + bedges = eptm.edge_df[eptm.edge_df["face"].isin(faces)] + face_cell = bedges.groupby("face")["cell"].first().to_dict() + adjacency = _cell_adjacency(eptm) + + edge_faces = {} + for face, srce, trgt in bedges[["face", "srce", "trgt"]].to_numpy(): + key = (int(min(srce, trgt)), int(max(srce, trgt))) + edge_faces.setdefault(key, set()).add(int(face)) + + cos_fold = np.cos(np.radians(theta_fold)) + pos = eptm.vert_df[eptm.coords] + candidates = [] + for (a, b), efaces in edge_faces.items(): + if len(efaces) < 2: + continue + e_mid = (pos.loc[a].to_numpy() + pos.loc[b].to_numpy()) / 2 + for g, f in itertools.combinations(sorted(efaces), 2): + cg, cf = face_cell[g], face_cell[f] + if cg == cf or frozenset((cg, cf)) in adjacency: + continue + r_g = centroids[g] - e_mid + r_f = centroids[f] - e_mid + ng, nf = np.linalg.norm(r_g), np.linalg.norm(r_f) + if ng == 0 or nf == 0: + continue + cos_phi = np.dot(r_g, r_f) / (ng * nf) + if cos_phi > cos_fold: # angle below theta_fold -> folded together + candidates.append((g, f, np.arccos(np.clip(cos_phi, -1, 1)))) + + candidates.sort(key=lambda c: c[2]) + out, used = [], set() + for g, f, _ in candidates: + if g in used or f in used: + continue + out.append((int(g), int(f))) + used.update((g, f)) + return out + + +def find_fusion_splits(eptm, theta_fold=20.0, segment="lateral"): + """Finds wide walls straddling a valence mismatch in a closing seam. + + A free `segment` face is returned when it has folded (within `theta_fold`) + against **two or more distinct** free `segment` faces belonging to different, + not-yet-adjacent cells -- i.e. a single wall apposed to two walls of the + other front (a "2-vs-1 pocket", left behind when two fusions trap an unequal + number of cells between them). Such a face must be split + (:func:`_split_lateral_face`) before the pocket can zip shut. Returns the + faces to split, widest mismatch first. + """ + faces = _free_lateral_faces(eptm, segment) + if len(faces) < 3: + return [] + + centroids, _ = _face_centroids_normals(eptm, faces) + bedges = eptm.edge_df[eptm.edge_df["face"].isin(faces)] + face_cell = bedges.groupby("face")["cell"].first().to_dict() + adjacency = _cell_adjacency(eptm) + + edge_faces = {} + for face, srce, trgt in bedges[["face", "srce", "trgt"]].to_numpy(): + key = (int(min(srce, trgt)), int(max(srce, trgt))) + edge_faces.setdefault(key, set()).add(int(face)) + + cos_fold = np.cos(np.radians(theta_fold)) + pos = eptm.vert_df[eptm.coords] + # for each free face, collect the distinct opposing cells it is folded against + opposing = {face: set() for face in faces} + for (a, b), efaces in edge_faces.items(): + if len(efaces) < 2: + continue + e_mid = (pos.loc[a].to_numpy() + pos.loc[b].to_numpy()) / 2 + for g, f in itertools.combinations(sorted(efaces), 2): + cg, cf = face_cell[g], face_cell[f] + if cg == cf or frozenset((cg, cf)) in adjacency: + continue + r_g = centroids[g] - e_mid + r_f = centroids[f] - e_mid + ng, nf = np.linalg.norm(r_g), np.linalg.norm(r_f) + if ng == 0 or nf == 0: + continue + if np.dot(r_g, r_f) / (ng * nf) > cos_fold: # folded together + opposing[g].add(cf) + opposing[f].add(cg) + + splits = [(face, len(cells)) for face, cells in opposing.items() if len(cells) >= 2] + splits.sort(key=lambda c: c[1], reverse=True) + return [int(face) for face, _ in splits] + + +def all_lateral_fusions(eptm, d_max=None, theta_face=45.0, theta_fold=20.0, + margin=0.5, segment="lateral", validate=True): + """Performs every available lateral-face fusion on `eptm`. + + Each round it extends existing seams (:func:`find_fusion_propagations`) and + nucleates new ones (:func:`find_fusion_nucleations`), applies one valid + :func:`fuse_lateral_faces`, then re-detects (indices change). When a wall + straddles a valence mismatch (a "2-vs-1 pocket", :func:`find_fusion_splits`) + it is divided with :func:`_split_lateral_face` so the pocket can finish + zipping. Mismatched or invalid fusions are skipped. + + Returns the number of fusions performed (splits are not counted). + """ + count = 0 + max_rounds = 4 * eptm.face_df.shape[0] + 10 + for _ in range(max_rounds): + splits = find_fusion_splits(eptm, theta_fold=theta_fold, segment=segment) + flagged = set(splits) + candidates = find_fusion_propagations( + eptm, theta_fold=theta_fold, segment=segment + ) + find_fusion_nucleations( + eptm, d_max=d_max, theta_face=theta_face, margin=margin, segment=segment + ) + # don't fuse a wall that needs splitting first -- it would mis-weld the + # far corner of the wide wall onto a single narrow partner. + candidates = [ + (g, f) for g, f in candidates if g not in flagged and f not in flagged + ] + progressed = False + for g, f in candidates: + if fuse_lateral_faces(eptm, g, f, validate=validate) == 0: + count += 1 + progressed = True + break # indices changed, re-detect + if not progressed and splits: + _split_lateral_face(eptm, splits[0]) + progressed = True # a piece can now fuse on the next round + if not progressed: + break + return count diff --git a/src/tyssue/topology/sheet_topology.py b/src/tyssue/topology/sheet_topology.py index 372533bc..f992fe02 100644 --- a/src/tyssue/topology/sheet_topology.py +++ b/src/tyssue/topology/sheet_topology.py @@ -54,6 +54,8 @@ def split_vert( sheet.reset_index() sheet.reset_topo() + sheet.network_changed = True + return new_edges @@ -170,6 +172,9 @@ def cell_division(sheet, mother, geom, angle=None): if sheet.settings.get("boundaries") is not None and mother_on_periodic_boundary: sheet.specs["settings"]["boundaries"] = saved_boundary geom.update_all(sheet) + + sheet.network_changed = True + return daughter @@ -240,7 +245,6 @@ def face_division(sheet, mother, vert_a, vert_b): sheet.reset_topo() return daughter - def drop_face(eptm, face, geom, **kwargs): """ Removes the face indexed by "face" and all associated edges to allow holes @@ -248,9 +252,10 @@ def drop_face(eptm, face, geom, **kwargs): edge = eptm.edge_df.loc[(eptm.edge_df['face'] == face)].index eptm.remove(edge, **kwargs) - eptm.sanitize(trim_borders = True) + eptm.sanitize(trim_borders = False) geom.update_all(eptm) + eptm.network_changed = True def resolve_t1s(sheet, geom, model, solver, max_iter=60):