diff --git a/.env b/.env index f440861..9102a85 100644 --- a/.env +++ b/.env @@ -2,8 +2,7 @@ PYTHONPATH=. # JAX_ENABLE_X64=1 # XLA_PYTHON_CLIENT_PREALLOCATE=0 -# CUDA_VISIBLE_DEVICES=0 -# CUDA_VISIBLE_DEVICES=-1 +# CUDA_VISIBLE_DEVICES=1 #-1 # JAX_PLATFORM_NAME=cpu # JAX_DISABLE_JIT=1 # JAX_DEBUG_NANS=1 diff --git a/.vscode/launch.json b/.vscode/launch.json index 67965f0..ff25232 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -50,6 +50,13 @@ "program": "${workspaceFolder}/examples/draft_examples.py", "justMyCode": false }, + { + "name": "Draft presentation", + "type": "python", + "request": "launch", + "program": "${workspaceFolder}/examples/draft_presentation.py", + "justMyCode": false + }, { "name": "Train graph", "type": "python", diff --git a/.vscode/settings.json b/.vscode/settings.json index 986dc34..573308f 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -1,7 +1,7 @@ { "python.formatting.blackArgs": [ - "--line-length=100", - "--target-version=py37" + // "--line-length=100", + // "--target-version=py37" ], // Pylint "python.linting.pylintArgs": [ diff --git a/README.md b/README.md index 2342b90..6911d2a 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ Contact mechanics models the behavior of physical bodies that come into contact with each other. It examines phenomena such as collisions, normal compliance, and friction. Most contact problems cannot be solved analytically and require a numerical procedure, such as the classical finite element method (FEM). -Conmech3d is an implementation of FEM for soft-body mechanical contact problems. The project is almost entirely self-contained and mainly aimed at research and didactic applications. Conmech3d is written in Python and uses [JAX](https://github.com/google/jax/tree/main), a library for high-performance numerical computing. Besides basic Python libraries, such as [Numpy](https://github.com/numpy/numpy) and [Scipy](https://scipy.org/), it also employs [pygmsh](https://github.com/meshpro/pygmsh) for mesh construction and [Numba](https://github.com/numba/numba) along with [Cython](https://github.com/cython/cython) to increase the speed of initial setup. Various options for visualization of simulation results are included, such as [Blender](https://github.com/blender/blender), [Three.js](https://github.com/mrdoob/three.js/) and [Matplotlib](https://github.com/matplotlib/matplotlib) +Conmech3d is an implementation of FEM for soft-body mechanical contact problems. The project is almost entirely self-contained and mainly aimed at research and didactic applications. Conmech3d is written in Python and uses [JAX](https://github.com/google/jax/tree/main), a library for high-performance numerical computing. Besides basic Python libraries, such as [Numpy](https://github.com/numpy/numpy) and [Scipy](https://scipy.org/), it also employs [pygmsh](https://github.com/meshpro/pygmsh) for mesh construction and [Numba](https://github.com/numba/numba) along with [Cython](https://github.com/cython/cython) to increase the speed of initial setup. Various options for visualization of simulation results are included, such as [Blender](https://github.com/blender/blender), [Three.js](https://github.com/mrdoob/three.js/) and [Matplotlib](https://github.com/matplotlib/matplotlib). Experimental implementations of model reduction techniques that include tetrahedral skinning used in computer graphics and a new approach using Graph Neural Network are included in this repository. diff --git a/conmech/dynamics/dynamics.py b/conmech/dynamics/dynamics.py index e21abea..9760b87 100644 --- a/conmech/dynamics/dynamics.py +++ b/conmech/dynamics/dynamics.py @@ -171,9 +171,11 @@ def get_rotation(self, displacement): displacement, self.matrices.dx_big_jax ) if not state.success: - raise Exception("Error calculating rotation") + raise ArithmeticError("Error calculating rotation") # print(state.iteration, state.norm) - return complete_base(base_seed=np.array(final_rotation, dtype=np.float64)) + return np.array( + complete_base(base_seed=np.array(final_rotation, dtype=np.float64)) + ) # def iterate_self(self, acceleration, temperature=None): # super().iterate_self(acceleration, temperature) diff --git a/conmech/helpers/cmh.py b/conmech/helpers/cmh.py index 16e8112..2f8bc92 100644 --- a/conmech/helpers/cmh.py +++ b/conmech/helpers/cmh.py @@ -98,11 +98,11 @@ def find_files_by_name(directory, name): def get_base_for_comarison(): - print("USING BASE FOR COMPARISON") all_paths = glob( "output/**/scenarios/*skinning_backwards*.scenes_comparer", recursive=True ) assert len(all_paths) == 1 + print("USING BASE FOR COMPARISON") return all_paths[0] diff --git a/conmech/helpers/config.py b/conmech/helpers/config.py index 69567f2..c6bd220 100644 --- a/conmech/helpers/config.py +++ b/conmech/helpers/config.py @@ -40,10 +40,10 @@ class SimulationConfig: use_nonconvex_friction_law: bool use_constant_contact_integral: bool use_lhs_preconditioner: bool - use_pca: bool with_self_collisions: bool mesh_layer_proportion: int = None - mode: str = "normal" # "normal" "skinning" "net" + mode: str = "normal" # "normal" "skinning" "net" "pca" + @dataclass @@ -51,7 +51,7 @@ class Config: shell: bool = False timestamp_skip: int = 10000 run_timestamp: float = int(time.time() * timestamp_skip) - current_time: str = datetime.now().strftime("%m.%d-%H.%M.%S") + current_time: str = datetime.now().strftime("%y.%m.%d-%H.%M.%S") verbose: bool = True animation_backend: str = "three" # "matplotlib blender three" diff --git a/conmech/helpers/interpolation_helpers.py b/conmech/helpers/interpolation_helpers.py index 3c3379a..9a19eec 100644 --- a/conmech/helpers/interpolation_helpers.py +++ b/conmech/helpers/interpolation_helpers.py @@ -126,14 +126,14 @@ def interpolate_3d_corner_vectors( nodes: np.ndarray, base: np.ndarray, corner_vectors: np.ndarray ): # orthonormal matrix; inverse equals transposition - upward_nodes = lnh.get_in_base(nodes, base.T) + upward_nodes = lnh.get_in_base2(nodes, base) scaled_nodes = scale_nodes_to_cube(upward_nodes) upward_vectors_interpolation = interpolate_scaled_nodes_numba( scaled_nodes=scaled_nodes, corner_vectors=corner_vectors, ) - vectors_interpolation = lnh.get_in_base(upward_vectors_interpolation, base) + vectors_interpolation = lnh.get_in_base2(upward_vectors_interpolation, base.T) # assert np.abs(np.mean(vectors_interpolation)) < 0.1 return vectors_interpolation diff --git a/conmech/helpers/lnh.py b/conmech/helpers/lnh.py index 6a58149..4ce9924 100644 --- a/conmech/helpers/lnh.py +++ b/conmech/helpers/lnh.py @@ -2,13 +2,10 @@ linear algebra helpers """ +import jax.numpy as jnp import numpy as np -from conmech.helpers import nph - - -def move_vector(vectors, index): - return np.roll(vectors, -index, axis=0) +from conmech.helpers import jxh, nph def complete_base(base_seed): @@ -19,17 +16,21 @@ def complete_base(base_seed): return base +def __move_vector(vectors, index): + return jnp.roll(vectors, -index, axis=0) + + def __orthonormalize_priority_gram_schmidt(base_seed, index): - prioritized_base_seed = move_vector(vectors=base_seed, index=index) + prioritized_base_seed = __move_vector(vectors=base_seed, index=index) prioritized_base = __orthonormalize_gram_schmidt(prioritized_base_seed) - base = move_vector(vectors=prioritized_base, index=index) + base = __move_vector(vectors=prioritized_base, index=index) return base def __orthonormalize_gram_schmidt(base_seed): - normalized_base_seed = nph.normalize_euclidean_numba(base_seed) + normalized_base_seed = jxh.normalize_euclidean(base_seed) unnormalized_base = __orthogonalize_gram_schmidt(normalized_base_seed) - base = nph.normalize_euclidean_numba(unnormalized_base) + base = jxh.normalize_euclidean(unnormalized_base) return base @@ -37,14 +38,14 @@ def __orthogonalize_gram_schmidt(vectors): # Gramm-Schmidt orthogonalization b0 = vectors[0] if len(vectors) == 1: - return np.array((b0)) + return jnp.array((b0)) b1 = vectors[1] - (vectors[1] @ b0) * b0 if len(vectors) == 2: - return np.array((b0, b1)) + return jnp.array((b0, b1)) - b2 = np.cross(b0, b1) - return np.array((b0, b1, b2)) + b2 = jnp.cross(b0, b1) + return jnp.array((b0, b1, b2)) def generate_base(dimension): @@ -81,5 +82,5 @@ def correct_base(base): return True -def get_in_base(vectors, base): - return vectors @ base.T +def get_in_base2(vectors, base): + return vectors @ base diff --git a/conmech/helpers/pca.py b/conmech/helpers/pca.py index ffa4a96..453f609 100644 --- a/conmech/helpers/pca.py +++ b/conmech/helpers/pca.py @@ -4,9 +4,13 @@ import jax import jax.numpy as jnp +import numpy as np from tqdm import tqdm -from conmech.helpers import cmh, nph +from conmech.helpers import cmh, lnh, nph +from conmech.properties.mesh_properties import MeshProperties +from conmech.scene.scene import Scene +from conmech.simulations.simulation_runner import create_scene def get_all_indices(data_path): @@ -37,11 +41,14 @@ def get_scenes(): input_path = "/home/michal/Desktop/conmech3d/output" scene_files = cmh.find_files_by_extension(input_path, "scenes") # scenes_data path_id = "/scenarios/" - scene_files = [f for f in scene_files if path_id in f] + scene_files = [f for f in scene_files if path_id in f and "SAVED" not in f] + assert len(scene_files) == 1 # all_arrays_path = max(scene_files, key=os.path.getctime) scenes = [] for all_arrays_path in scene_files: + if "SAVED" in all_arrays_path: + continue all_arrays_name = os.path.basename(all_arrays_path).split("DATA")[0] print(f"FILE: {all_arrays_name}") @@ -57,34 +64,6 @@ def get_scenes(): return scenes -def get_projection(data, latent_dim=200): - projection_mean = 0 * data.mean(axis=0) # columnwise mean = 0 - svd = jax.numpy.linalg.svd(data - projection_mean, full_matrices=False) - # (svd[0] @ jnp.diag(svd[1]) @ svd[2]) - projection_matrix = svd[2][:latent_dim].T - return {"matrix": projection_matrix, "mean": projection_mean.reshape(-1, 1)} - - -def project_to_latent(projection, data_stack): - data_stack_zeroed = data_stack - projection["mean"] - latent = projection["matrix"].T @ data_stack_zeroed - return latent - - -def project_from_latent(projection, latent): - data_stack_zeroed = projection["matrix"] @ latent - data_stack = data_stack_zeroed + projection["mean"] - return data_stack - - -def p_to_vector(projection, vector): - return project_to_latent(projection, vector.reshape(-1, 1)).reshape(-1) - - -def p_from_vector(projection, vector): - return project_from_latent(projection, vector.reshape(-1, 1)).reshape(-1) - - def save_pca(projection, file_path="./output/PCA"): with open(file_path, "wb") as file: pickle.dump(projection, file) @@ -96,45 +75,104 @@ def load_pca(file_path="./output/PCA"): return projection +def get_displacement_new(scene): + velocity = scene.velocity_old + scene.time_step * scene.exact_acceleration + displacement = scene.displacement_old + scene.time_step * velocity + return displacement + + def get_data_scenes(scenes): data_list = [] - count = len(scenes) - for scene in scenes: - u = jnp.array(scene.get_last_displacement_step()) # scene.displacement_old) - u_stack = nph.stack_column(u) + for scene in tqdm(scenes): + # print(scene.moved_base) + u = scene.get_lifted_displacement() + u_stack = nph.stack(u) data_list.append(u_stack) - data = jnp.array(data_list).reshape(count, -1) + data = jnp.array(data_list) return data, u_stack, u -def get_data_dataset(dataloader): +def get_data_dataset(dataloader, scene): data_list = [] - count = 1000 - for _ in tqdm(range(count)): - sample = next(iter(dataloader)) + count = 3000 + print(f"LIMIT TO {count}") + for i, sample in enumerate(tqdm(dataloader)): # check randomness target = sample[0][1] - u = jnp.array(target.reduced_acceleration) - u_stack = nph.stack_column(u) - data_list.append(u_stack) + original_displacement = jnp.array(target["new_displacement"]) - data = jnp.array(data_list).reshape(count, -1) - return data, u_stack, u + original_rotation = scene.get_rotation(original_displacement) + random_rotation = jnp.linalg.qr(np.random.rand(3, 3))[0] + new_rotation = original_rotation.T @ random_rotation + + moved_nodes = scene.initial_nodes + original_displacement + displacement_mean = np.mean(moved_nodes, axis=0) + rotated_moved_nodes = lnh.get_in_base2( + (moved_nodes - displacement_mean), new_rotation + ) + displacement = rotated_moved_nodes - scene.initial_nodes + displacement += displacement_mean # 20 * np.random.rand(3) ### + + displacement_stack = nph.stack(displacement) + data_list.append(displacement_stack) + if i > count: + break + + data = jnp.array( + data_list + ) # Sort by displacement and get max, plot hist # np.linalg.norm(data_list[190]) + return data, displacement_stack, displacement + + +def get_projection(data, latent_dim): + projection_mean = 0 # data.mean(axis=0) + + svd = jax.numpy.linalg.svd(data - projection_mean, full_matrices=False) + # (svd[0] @ jnp.diag(svd[1]) @ svd[2]) + + projection_matrix = svd[2][:latent_dim] + # projection_matrix = jax.experimental.sparse.eye(data.shape[1]) + + return {"matrix": projection_matrix, "mean": projection_mean} + + +def project_to_latent(projection, data): + data_zeroed = data - projection["mean"] + latent = projection["matrix"] @ data_zeroed + return latent + + +def project_from_latent(projection, latent): + data_stack_zeroed = projection["matrix"].T @ latent + data_stack = data_stack_zeroed + projection["mean"] + return data_stack + + +def p_to_vector(projection, data): + return project_to_latent(projection, nph.stack(data)) + + +def p_from_vector(projection, latent): + return nph.unstack(project_from_latent(projection, latent), dim=3) -def run(dataloader): - _ = dataloader - scenes = get_scenes() - data, sample_u_stack, sample_u = get_data_scenes(scenes) - # data, sample_u_stack, sample_u = get_data_dataset(dataloader) +def run(dataloader, latent_dim, scenario): + if dataloader is None: + scenes = get_scenes() + data, sample_u_stack, sample_u = get_data_scenes(scenes) + else: + scene = create_scene(scenario) + data, sample_u_stack, sample_u = get_data_dataset( + dataloader=dataloader, scene=scene + ) - original_projection = get_projection(data) + original_projection = get_projection(data, latent_dim) save_pca(original_projection) - projection = load_pca() - latent = project_to_latent(projection, sample_u_stack) - u_reprojected_stack = project_from_latent(projection, latent) - u_reprojected = nph.unstack(u_reprojected_stack, dim=3) - print("Error max: ", jnp.abs(u_reprojected - sample_u).max()) + # projection = load_pca() + # latent = project_to_latent(projection, sample_u_stack) + # u_reprojected_stack = project_from_latent(projection, latent) + # u_reprojected = nph.unstack(u_reprojected_stack, dim=3) + # print("Error max: ", jnp.abs(u_reprojected - sample_u).max()) return 0 diff --git a/conmech/mesh/mesh_builders.py b/conmech/mesh/mesh_builders.py index 09076d8..20555db 100644 --- a/conmech/mesh/mesh_builders.py +++ b/conmech/mesh/mesh_builders.py @@ -44,7 +44,7 @@ def translate_nodes(nodes: np.ndarray, mesh_prop: MeshProperties): if mesh_prop.mean_at_origin: nodes -= np.mean(nodes, axis=0) if mesh_prop.initial_base is not None: - nodes = lnh.get_in_base(nodes, mesh_prop.initial_base) + nodes = lnh.get_in_base2(nodes, mesh_prop.initial_base.T) if mesh_prop.initial_position is not None: nodes += mesh_prop.initial_position return nodes diff --git a/conmech/plotting/plotter_2d.py b/conmech/plotting/plotter_2d.py index 3a22784..6545060 100644 --- a/conmech/plotting/plotter_2d.py +++ b/conmech/plotting/plotter_2d.py @@ -7,13 +7,18 @@ from conmech.helpers.config import Config from conmech.plotting import plotter_common -from conmech.plotting.plotter_common import PlotAnimationConfig, make_animation +from conmech.plotting.plotter_common import ( + PlotAnimationConfig, + get_recentered_nodes, + make_animation, +) from conmech.scene.scene import Scene from conmech.scene.scene_temperature import SceneTemperature +from conmech.state.body_position import mesh_normalization_decorator def get_fig(): - return plt.figure(figsize=(4, 2)) + return plt.figure(figsize=(4, 3)) def get_axs(fig): @@ -24,9 +29,8 @@ def get_axs(fig): def set_perspective(scale, axes): axes.set_aspect("equal", "box") - padding = 6 - axes.set_xlim(-padding * scale, 18 * scale) - axes.set_ylim(-padding * scale, padding * scale) + axes.set_xlim(-3 * scale, 14 * scale) + axes.set_ylim(-6 * scale, 6 * scale) plotter_common.set_ax(axes) @@ -54,10 +58,11 @@ def plot_animation( ) +@mesh_normalization_decorator def plot_frame( + scene: Scene, fig, axs, - scene: Scene, current_time: float, draw_detailed: bool = True, t_scale: Optional[np.ndarray] = None, @@ -80,31 +85,31 @@ def plot_frame( position = np.array([-3.7, 4.2]) * scale draw_all_sparse(scene, position, axes=axes) - position = np.array([-6.2, -4.2]) * scale + position = np.array([-0.75, -4.2]) * scale shift = 2.5 * scale - position[0] += shift - draw_initial(scene, position, axes=axes) + # position[0] += shift + # draw_initial(scene, position, axes=axes) position[0] += shift draw_forces(scene, position, axes=axes) - position[0] += shift - draw_obstacle_resistance_normalized(scene, position, axes=axes) - position[0] += shift + # position[0] += shift + # draw_obstacle_resistance_normalized(scene, position, axes=axes) + # position[0] += shift # draw_boundary_surfaces_normals(scene, position, axes) # position[0] += shift # draw_boundary_normals(scene, position, axes) - # position[0] += shift + position[0] += shift draw_boundary_resistance_normal(scene, position, axes=axes) position[0] += shift draw_boundary_resistance_tangential(scene, position, axes=axes) - position[0] += shift - draw_boundary_v_tangential(scene, position, axes=axes) - position[0] += shift + # position[0] += shift + # draw_boundary_v_tangential(scene, position, axes=axes) - draw_input_u(scene, position, axes=axes) - position[0] += shift - draw_input_v(scene, position, axes=axes) + # position[0] += shift + # draw_input_u(scene, position, axes=axes) + # position[0] += shift + # draw_input_v(scene, position, axes=axes) position[0] += shift draw_a(scene, position, axes=axes) @@ -326,15 +331,15 @@ def draw_nodes(nodes, position, color, axes): def draw_forces(scene: Scene, position, axes): - return draw_data("F", scene.normalized_inner_forces, scene, position, axes) + return draw_data("F", scene.inner_forces, scene, position, axes) def draw_input_u(scene: Scene, position, axes): - return draw_data("U", scene.normalized_displacement_old, scene, position, axes) + return draw_data("U", scene.displacement_old, scene, position, axes) def draw_input_v(scene: Scene, position, axes): - return draw_data("V", scene.normalized_velocity_old, scene, position, axes) + return draw_data("V", scene.velocity_old, scene, position, axes) def draw_a(scene, position, axes): @@ -349,11 +354,11 @@ def draw_a(scene, position, axes): def draw_data(annotation, data, scene: Scene, position, axes): draw_moved_body(annotation, scene, position, axes) - plot_arrows(scene.normalized_nodes + position, data, axes) + plot_arrows(get_recentered_nodes(scene, position), data, axes) def draw_moved_body(annotation, scene: Scene, position, axes): - draw_triplot(scene.normalized_nodes + position, scene, "tab:blue", axes) + draw_triplot(get_recentered_nodes(scene, position), scene, "tab:blue", axes) add_annotation(annotation, scene, position, axes) diff --git a/conmech/plotting/plotter_3d.py b/conmech/plotting/plotter_3d.py index 28e7d2f..e4b26ba 100644 --- a/conmech/plotting/plotter_3d.py +++ b/conmech/plotting/plotter_3d.py @@ -10,6 +10,7 @@ from conmech.plotting.plotter_common import PlotAnimationConfig, make_animation from conmech.scene.scene import Scene from conmech.scene.scene_temperature import SceneTemperature +from conmech.state.body_position import mesh_normalization_decorator def get_fig(): @@ -67,10 +68,11 @@ def get_axs(fig): return [ax0, ax1, ax2, ax3] +@mesh_normalization_decorator def plot_frame( + scene: Scene, fig, axs, - scene: Scene, current_time, t_scale: Optional[List] = None, ): diff --git a/conmech/plotting/plotter_common.py b/conmech/plotting/plotter_common.py index 6345e49..e2f595f 100644 --- a/conmech/plotting/plotter_common.py +++ b/conmech/plotting/plotter_common.py @@ -32,6 +32,10 @@ def mappable(self): return plt.cm.ScalarMappable(norm=norm, cmap=self.cmap) +def get_recentered_nodes(scene, position): + return scene.moved_nodes - np.mean(scene.moved_nodes, axis=0) + position + + def get_t_scale( scenario: Scenario, index_skip: int, @@ -138,9 +142,9 @@ def animate(step: int, args: AnimationArgs): ) plot_frame( + scene=scene, axs=axs, fig=args.fig, - scene=scene, current_time=step * args.time_skip, t_scale=t_scale, ) diff --git a/conmech/plotting/plotter_functions.py b/conmech/plotting/plotter_functions.py index 16665e8..2df5a4c 100644 --- a/conmech/plotting/plotter_functions.py +++ b/conmech/plotting/plotter_functions.py @@ -88,9 +88,6 @@ def save_results_three(file_path, json_dict): with open(file_path, "w", encoding="utf-8") as file: json.dump(json_dict, file) - list_path = "output/three_list.json" - cmh.clear_file(list_path) - folder_list = [f[0] for f in os.walk("output") if "0.json" in f[2]] # [1:] folder_list.sort(reverse=True) simulations_list, step_list = [], [] @@ -100,12 +97,14 @@ def save_results_three(file_path, json_dict): simulations_list.append(simulation) # os.path.basename(simulation)) step_list.append(max(steps)) + list_path = "output/three_list.json" + # cmh.clear_file(list_path) with open(list_path, "w", encoding="utf-8") as file: json.dump({"simulations": simulations_list, "steps": step_list}, file) def plot_using_blender(output: bool = True): - path = "~/Desktop/Blender/blender-3.2.0-linux-x64/blender" # TODO: Move to configuration + path = "~/Desktop/Blender/blender-4.1.1-linux-x64/blender" # TODO: Move to configuration args = " --background --python ~/Desktop/conmech3d/blender/load.py --render" print("Plotting using Blender...") stdout = sys.stdout if output else subprocess.DEVNULL @@ -150,9 +149,9 @@ def plot_setting( fig = plotter_2d.get_fig() axs = plotter_2d.get_axs(fig) plotter_2d.plot_frame( + scene=scene, fig=fig, axs=axs, - scene=scene, current_time=current_time, draw_detailed=draw_detailed, ) @@ -160,5 +159,5 @@ def plot_setting( else: fig = plotter_3d.get_fig() axs = plotter_3d.get_axs(fig) - plotter_3d.plot_frame(fig=fig, axs=axs, scene=scene, current_time=current_time) + plotter_3d.plot_frame(scene=scene, fig=fig, axs=axs, current_time=current_time) plt_save(path, extension) diff --git a/conmech/scene/energy_functions.py b/conmech/scene/energy_functions.py index a423c65..376ecbe 100644 --- a/conmech/scene/energy_functions.py +++ b/conmech/scene/energy_functions.py @@ -6,7 +6,7 @@ import numpy as np from conmech.dynamics.dynamics import _get_deform_grad -from conmech.helpers import jxh, nph +from conmech.helpers import jxh, lnh, nph from conmech.helpers.config import SimulationConfig @@ -419,53 +419,29 @@ def compute_velocity_energy( self.temperature_cost_function = None - # return - - # def to_displacement(function): - # return lambda displacement, args: function( - # nph.displacement_to_acceleration(displacement, args), - # args, - # ) * (1 / (args.time_step**2)) - - # if not simulation_config.use_pca: - # def to_displacement_by_factor(function): - # return lambda disp_by_factor, args: to_displacement(function)( - # disp_by_factor * (args.time_step**2), args - # ) - - # self.energy_obstacle_free = to_displacement_by_factor(self._energy_obstacle_free) - # self.energy_obstacle_colliding = to_displacement_by_factor( - # self._energy_obstacle_colliding - # ) - - # # self.energy_obstacle_free = self._energy_obstacle_free - # # self.energy_obstacle_colliding = self._energy_obstacle_colliding - - # # self.energy_obstacle_free = lambda vector, args: jnp.float64( - # # self._energy_obstacle_free(jnp.array(vector, dtype=jnp.float32), args) - # # ) - # # self.energy_obstacle_colliding = lambda vector, args: jnp.float64( - # # self._energy_obstacle_colliding(jnp.array(vector, dtype=jnp.float32), args) - # # ) - - # else: - # projection = pca.load_pca() - - # def to_displacement_by_factor_pca(function): - # return lambda disp_by_factor, args: to_displacement(function)( - # (pca.p_from_vector( - # projection, - # pca.p_to_vector(projection, disp_by_factor\ - # - nph.stack_column(args.displacement_old).reshape(-1)), - # ) + nph.stack_column(args.displacement_old).reshape(-1)) - # * (args.time_step**2), - # args, - # ) - - # self.energy_obstacle_free = to_displacement_by_factor_pca(self._energy_obstacle_free) - # self.energy_obstacle_colliding = to_displacement_by_factor_pca( - # self._energy_obstacle_colliding - # ) + if simulation_config.mode == "pca": + from conmech.helpers.pca import load_pca, p_from_vector + + self.projection = load_pca() + + def to_displacement_by_factor_pca(energy_function): + def reformulation(u_latent, args): + u_projected = p_from_vector(self.projection, u_latent) + u_projected_vector = nph.stack(u_projected) + + return energy_function( + nph.displacement_to_acceleration(u_projected_vector, args), + args, + ) + + return reformulation + + self.energy_obstacle_free = to_displacement_by_factor_pca( + self._energy_obstacle_free + ) + self.energy_obstacle_colliding = to_displacement_by_factor_pca( + self._energy_obstacle_colliding + ) @staticmethod def get_manual_modes(): diff --git a/conmech/scene/scene.py b/conmech/scene/scene.py index 74a0990..42a94fb 100644 --- a/conmech/scene/scene.py +++ b/conmech/scene/scene.py @@ -10,7 +10,6 @@ from conmech.helpers import jxh, lnh, nph from conmech.helpers.config import SimulationConfig from conmech.helpers.interpolation_helpers import interpolate_nodes -from conmech.helpers.lnh import get_in_base from conmech.properties.body_properties import TimeDependentBodyProperties from conmech.properties.mesh_properties import MeshProperties from conmech.properties.obstacle_properties import ObstacleProperties @@ -343,20 +342,6 @@ def normalized_exact_acceleration(self): def normalized_lifted_acceleration(self): return self.normalize_rotate(self.lifted_acceleration) - @mesh_normalization_decorator - def force_denormalize(self, acceleration): - return self.denormalize_rotate(acceleration) - - # @property - # @mesh_normalization_decorator - # def norm_exact_new_displacement(self): - # return self.to_normalized_displacement(self.exact_acceleration) - - # @property - # @mesh_normalization_decorator - # def norm_lifted_new_displacement(self): - # return self.to_normalized_displacement(self.lifted_acceleration) - @property @mesh_normalization_decorator def norm_by_reduced_lifted_new_displacement(self): @@ -369,9 +354,9 @@ def _normalize_current_reduced(moved_nodes): base_scene = self.reduced else: base_scene = self # TODO: Add warning - return get_in_base( + return lnh.get_in_base2( (moved_nodes - np.mean(base_scene.moved_nodes, axis=0)), - base_scene.get_rotation(base_scene.displacement_old), + base_scene.get_rotation(base_scene.displacement_old).T, ) displacement_new = self.to_displacement(exact_acceleration) @@ -392,11 +377,10 @@ def from_displacement(self, displacement): @mesh_normalization_decorator def to_normalized_displacement(self, acceleration): displacement_new = self.to_displacement(acceleration) - moved_nodes_new = self.initial_nodes + displacement_new - new_normalized_nodes = get_in_base( + new_normalized_nodes = lnh.get_in_base2( (moved_nodes_new - np.mean(moved_nodes_new, axis=0)), - self.get_rotation(displacement_new), + self.get_rotation(displacement_new).T, ) return new_normalized_nodes - self.normalized_initial_nodes @@ -405,9 +389,9 @@ def to_normalized_displacement(self, acceleration): # displacement_new = self.to_displacement(acceleration) # moved_nodes_new = self.initial_nodes + displacement_new - # new_normalized_nodes = get_in_base( + # new_normalized_nodes = lnh.get_in_base2( # (moved_nodes_new - np.mean(moved_nodes_new, axis=0)), - # self.get_rotation(self.displacement_old), + # self.get_rotation(self.displacement_old).T, # ) # assert np.allclose(new_normalized_nodes, self.normalize_shift_and_rotate(moved_nodes_new)) # return new_normalized_nodes - self.normalized_initial_nodes @@ -416,9 +400,9 @@ def to_normalized_displacement(self, acceleration): def to_normalized_displacement_rotated_displaced(self, acceleration): displacement_new = self.to_displacement(acceleration) moved_nodes_new = self.initial_nodes + displacement_new - new_normalized_nodes = get_in_base( + new_normalized_nodes = lnh.get_in_base2( (moved_nodes_new - np.mean(self.moved_nodes, axis=0)), - self.get_rotation(self.displacement_old), + self.get_rotation(self.displacement_old).T, ) # assert np.allclose( # new_normalized_nodes, @@ -436,16 +420,42 @@ def to_normalized_displacement_rotated_displaced(self, acceleration): # acceleration = self.from_displacement(displacement_new) # return acceleration - def get_last_displacement_step(self): - return self.displacement_old - self.time_step * self.velocity_old + @mesh_normalization_decorator + def force_normalize_pca(self, displacement_new): + # return displacement_new + return displacement_new - np.mean(self.moved_nodes, axis=0) + moved_nodes_new = self.initial_nodes + displacement_new + new_normalized_nodes = self.normalize_rotate( + moved_nodes_new - np.mean(self.moved_nodes, axis=0) + ) + # rotation = self.get_rotation_pca(displacement_new) + # new_normalized_nodes = lnh.get_in_base2( + # (moved_nodes_new - np.mean(self.moved_nodes, axis=0)), rotation.T + # ) + return new_normalized_nodes - self.normalized_initial_nodes + + @mesh_normalization_decorator + def force_denormalize_pca(self, normalized_displacement_new): + # return normalized_displacement_new + return normalized_displacement_new + np.mean(self.moved_nodes, axis=0) + new_normalized_nodes = ( + normalized_displacement_new + self.normalized_initial_nodes + ) + moved_nodes_new = self.denormalize_rotate(new_normalized_nodes) + np.mean( + self.moved_nodes, axis=0 + ) + # moved_nodes_new = lnh.get_in_base2( + # new_normalized_nodes, rotation.T + # ) + np.mean(self.moved_nodes, axis=0) + return moved_nodes_new - self.initial_nodes - def displacement_from_step(self, displacement_step): - return self.displacement_old + displacement_step + def get_lifted_displacement(self): + return self.to_displacement(self.lifted_acceleration) def get_centered_nodes(self, displacement): nodes = self.centered_initial_nodes + displacement - centered_nodes = lnh.get_in_base( - (nodes - nodes.mean(axis=0)), self.get_rotation(displacement) + centered_nodes = lnh.get_in_base2( + (nodes - nodes.mean(axis=0)), self.get_rotation(displacement).T ) return centered_nodes @@ -454,9 +464,7 @@ def get_displacement(self, base, position, base_displacement=None): centered_nodes = self.centered_nodes else: centered_nodes = self.get_centered_nodes(base_displacement) - moved_centered_nodes = ( - lnh.get_in_base(centered_nodes, np.linalg.inv(base)) + position - ) + moved_centered_nodes = lnh.get_in_base2(centered_nodes, base) + position displacement = moved_centered_nodes - self.centered_initial_nodes return displacement diff --git a/conmech/simulations/problem_solver.py b/conmech/simulations/problem_solver.py index d799436..c6c7ba0 100644 --- a/conmech/simulations/problem_solver.py +++ b/conmech/simulations/problem_solver.py @@ -81,7 +81,6 @@ def __init__(self, setup: Problem, body_properties: BodyProperties): use_constant_contact_integral=False, use_lhs_preconditioner=False, with_self_collisions=False, - use_pca=False, ), dynamics_config=DynamicsConfiguration( create_in_subprocess=False, diff --git a/conmech/simulations/simulation_runner.py b/conmech/simulations/simulation_runner.py index d1cae96..7fa64e7 100644 --- a/conmech/simulations/simulation_runner.py +++ b/conmech/simulations/simulation_runner.py @@ -20,6 +20,8 @@ def get_solve_function(simulation_config): if simulation_config.mode == "normal": return Calculator.solve + if simulation_config.mode == "pca": + return Calculator.solve_skinning_backwards_base if simulation_config.mode == "compare_reduced": return Calculator.solve_compare_reduced if simulation_config.mode == "skinning": @@ -59,7 +61,7 @@ def create_scene(scenario): create_in_subprocess = False def get_scene(): - if scenario.simulation_config.mode == "normal": + if scenario.simulation_config.mode in ["normal"]: scene = Scene( mesh_prop=scenario.mesh_prop, body_prop=scenario.body_prop, @@ -68,7 +70,11 @@ def get_scene(): create_in_subprocess=create_in_subprocess, simulation_config=scenario.simulation_config, ) - elif scenario.simulation_config.mode in ["skinning", "skinning_backwards"]: + elif scenario.simulation_config.mode in [ + "skinning", + "skinning_backwards", + "pca", + ]: from deep_conmech.scene.scene_layers import SceneLayers scene = SceneLayers( @@ -180,24 +186,26 @@ def save_scene(scene: Scene, scenes_path: str, save_animation: bool): pkh.append_data(data=blender_data, data_path=blender_data_path, lock=None) # Comparer - comparer_data_path = scenes_path + "_comparer" + if hasattr(scene, "reduced"): + comparer_data_path = scenes_path + "_comparer" + + comparer_data = { + "exact_acceleration": scene.exact_acceleration, + "lifted_acceleration": scene.lifted_acceleration, + "lifted_displacement": scene.get_lifted_displacement(), + "norm_new_displacement": scene.to_normalized_displacement( + scene.lifted_acceleration + ), + "norm_reduced": scene.norm_by_reduced_lifted_new_displacement, + "normalized_by_reduces_nodes": scene.normalized_initial_nodes + + scene.norm_by_reduced_lifted_new_displacement, + "norm_lifted_new_displacement": scene.norm_lifted_new_displacement, + "normalized_nodes": scene.normalized_initial_nodes + + scene.norm_lifted_new_displacement, + "recentered_norm_lifted_new_displacement": scene.recentered_norm_lifted_new_displacement, + } - normalized_nodes = ( - scene.initial_nodes + scene.norm_by_reduced_lifted_new_displacement - ) - comparer_data = { - "displacement_old": scene.displacement_old, - "exact_acceleration": scene.exact_acceleration, - "normalized_nodes": normalized_nodes, - "lifted_acceleration": scene.lifted_acceleration, - "norm_lifted_new_displacement": scene.norm_lifted_new_displacement, - "recentered_norm_lifted_new_displacement": scene.recentered_norm_lifted_new_displacement, - "norm_reduced": scene.get_norm_by_reduced_lifted_new_displacement( - scene.exact_acceleration - ), - } - - pkh.append_data(data=comparer_data, data_path=comparer_data_path, lock=None) + pkh.append_data(data=comparer_data, data_path=comparer_data_path, lock=None) # Matplotlib if save_animation: @@ -254,12 +262,12 @@ def operation_save(scene: Scene): save_scene( scene=scene, scenes_path=scenes_path, save_animation=save_animation ) - if with_reduced: - save_scene( - scene=scene.reduced, - scenes_path=scenes_path_reduced, - save_animation=save_animation, - ) + # if with_reduced: + # save_scene( + # scene=scene.reduced, + # scenes_path=scenes_path_reduced, + # save_animation=save_animation, + # ) if plot_index: plot_scenes_count[0] += 1 step[0] += 1 diff --git a/conmech/solvers/algorithms/lbfgs.py b/conmech/solvers/algorithms/lbfgs.py index dc83493..d014e44 100644 --- a/conmech/solvers/algorithms/lbfgs.py +++ b/conmech/solvers/algorithms/lbfgs.py @@ -1,4 +1,5 @@ """The Limited-Memory Broyden-Fletcher-Goldfarb-Shanno minimization algorithm.""" +# Based on JAX repository # pylint: skip-file import os from functools import partial diff --git a/conmech/solvers/calculator.py b/conmech/solvers/calculator.py index 2d32f76..5cb15d8 100644 --- a/conmech/solvers/calculator.py +++ b/conmech/solvers/calculator.py @@ -54,6 +54,7 @@ def set_compiled_optimization_functions(energy_functions, hes_inv, x0, args): sample_x0=x0, sample_args=args, ) + # energy_functions.energy_obstacle_free(disp_by_factor=jnp.array([2,0.,0.]).reshape(-1,3).repeat(919, axis=0).reshape(-1), args=args, normalize=normalize, denormalize=denormalize) def set_and_get_opti_fun(energy_functions, scene, hes_inv, x0, args): @@ -105,22 +106,37 @@ def minimize_jax( return np.asarray(state.x_k) # , state @staticmethod - def minimize_jax_displacement( - function, initial_vector: np.ndarray, args, hes_inv, verbose: bool = True + def minimize_jax_displacement_pca( + initial_vector: np.ndarray, + args, + hes_inv, + scene, + energy_functions, + verbose: bool = True, ) -> np.ndarray: - range_factor = args.time_step**2 - initial_disp_by_factor = ( - nph.acceleration_to_displacement(initial_vector, args) / range_factor - ) - disp_by_factor = Calculator.minimize_jax( - initial_vector=initial_disp_by_factor, + from conmech.helpers.pca import p_from_vector, p_to_vector + + initial_disp = nph.acceleration_to_displacement(initial_vector, args) + + dim = 3 + initial_u = nph.unstack(initial_disp, dim=dim) + + initial_u_latent = p_to_vector(energy_functions.projection, initial_u) + + u_latent = Calculator.minimize_jax( + initial_vector=initial_u_latent, args=args, hes_inv=hes_inv, - function=function, + scene=scene, + energy_functions=energy_functions, verbose=verbose, ) - return nph.displacement_to_acceleration( - np.asarray(disp_by_factor * range_factor), args + + u_projected = p_from_vector(energy_functions.projection, u_latent) + u_projected_vector = nph.stack(u_projected) + + return np.asarray( + nph.displacement_to_acceleration(np.asarray(u_projected_vector), args) ) @staticmethod @@ -206,6 +222,24 @@ def solve_skinning_backwards( initial_a, initial_t=None, timer=Timer(), + ): + return Calculator.solve_skinning_backwards_base( + scene=scene, + energy_functions=energy_functions, + initial_a=initial_a, + initial_t=initial_t, + timer=timer, + with_base_for_comparison=False, + ) + + @staticmethod + def solve_skinning_backwards_base( + scene: Scene, + energy_functions: EnergyFunctions, + initial_a, + initial_t=None, + timer=Timer(), + with_base_for_comparison=True, ): _ = initial_a, initial_t energy_functions = ( @@ -213,14 +247,21 @@ def solve_skinning_backwards( if hasattr(energy_functions, "__len__") else energy_functions ) + with timer["reduced_solver"]: - exact_acceleration, _ = Calculator.solve( + scene.lifted_acceleration, _ = Calculator.solve( scene=scene, energy_functions=energy_functions, initial_a=scene.exact_acceleration, timer=timer, ) - scene.lifted_acceleration = exact_acceleration + if with_base_for_comparison: + dense_path = cmh.get_base_for_comarison() + exact_acceleration, _ = cmh.get_exact_acceleration( + scene=scene, path=dense_path + ) + else: + exact_acceleration = scene.lifted_acceleration with timer["lift_data"]: scene.reduced.exact_acceleration = scene.lift_acceleration_from_position( exact_acceleration @@ -471,8 +512,12 @@ def solve_acceleration_normalized_optimization_jax( ) with timer["__minimize_jax"]: + if not scene.simulation_config.mode == "pca": + minimize_fun = Calculator.minimize_jax + else: + minimize_fun = Calculator.minimize_jax_displacement_pca normalized_a_vector_np = cmh.profile( - lambda: Calculator.minimize_jax( # _displacement( + lambda: minimize_fun( # solver= energy_functions.get_solver(scene), initial_vector=initial_a_vector, args=args, diff --git a/conmech/state/body_position.py b/conmech/state/body_position.py index e5d9a18..59c5a15 100644 --- a/conmech/state/body_position.py +++ b/conmech/state/body_position.py @@ -125,15 +125,15 @@ def get_surface_per_boundary_node_jax( def mesh_normalization_decorator(func: Callable): - def inner(self, *args, **kwargs): - saved_normalize = self.normalize - self.normalize = True - if hasattr(self, "reduced"): - self.reduced.normalize = True - returned_value = func(self, *args, **kwargs) - self.normalize = saved_normalize - if hasattr(self, "reduced"): - self.reduced.normalize = saved_normalize + def inner(scene: "BodyPosition", *args, **kwargs): + saved_normalize = scene.normalize + scene.normalize = True + if hasattr(scene, "reduced"): + scene.reduced.normalize = True + returned_value = func(scene, *args, **kwargs) + scene.normalize = saved_normalize + if hasattr(scene, "reduced"): + scene.reduced.normalize = saved_normalize return returned_value return inner @@ -350,12 +350,12 @@ def iterate_self(self, acceleration, temperature=None): def normalize_rotate(self, vectors): if not self.normalize: return vectors - return lnh.get_in_base(vectors, self.moved_base) + return lnh.get_in_base2(vectors, self.moved_base.T) def denormalize_rotate(self, vectors): if not self.normalize: return vectors - return lnh.get_in_base(vectors, np.linalg.inv(self.moved_base)) + return lnh.get_in_base2(vectors, self.moved_base) def normalize_shift_and_rotate(self, vectors): return self.normalize_rotate(self._normalize_shift(vectors)) @@ -437,6 +437,6 @@ def input_displacement_old(self): @property def centered_nodes(self): - return lnh.get_in_base( - (self.moved_nodes - np.mean(self.moved_nodes, axis=0)), self.moved_base + return lnh.get_in_base2( + (self.moved_nodes - np.mean(self.moved_nodes, axis=0)), self.moved_base.T ) diff --git a/deep_conmech/data/base_dataset.py b/deep_conmech/data/base_dataset.py index 99de60f..560c2db 100644 --- a/deep_conmech/data/base_dataset.py +++ b/deep_conmech/data/base_dataset.py @@ -534,7 +534,7 @@ def solve_and_prepare_scene( # scene.reduced.exact_acceleration # ) - dense_path = cmh.get_base_for_comarison() + dense_path = None # cmh.get_base_for_comarison() if dense_path is not None: ( scene.exact_acceleration, diff --git a/deep_conmech/run_model.py b/deep_conmech/run_model.py index d83f08a..8b2635b 100644 --- a/deep_conmech/run_model.py +++ b/deep_conmech/run_model.py @@ -2,6 +2,10 @@ # if name == "__main__": # SET_ENV() +from dotenv import load_dotenv + +load_dotenv() + import argparse import os from argparse import ArgumentParser, Namespace @@ -13,7 +17,6 @@ import torch import torch.distributed as dist import torch.multiprocessing -from dotenv import load_dotenv from conmech.helpers import cmh, pca from conmech.helpers.config import Config, SimulationConfig @@ -156,10 +159,8 @@ def plot(config: TrainingConfig): def run_pca(config: TrainingConfig): - dataset = get_train_dataset(config.td.dataset, config=config) - dataset.initialize_data() - dataloader = base_dataset.get_train_dataloader(dataset) - + mesh_density = 32 # 16 #32 + final_time = 2.0 # 8.0 # 2.0 simulation_config = SimulationConfig( use_normalization=False, use_linear_solver=False, @@ -168,43 +169,54 @@ def run_pca(config: TrainingConfig): use_constant_contact_integral=False, use_lhs_preconditioner=False, with_self_collisions=True, - use_pca=False, + mode="pca", ) - final_time = 0.5 all_scenarios = [ scenarios.bunny_fall_3d( - mesh_density=32, + mesh_density=mesh_density, scale=1, final_time=final_time, simulation_config=simulation_config, ), scenarios.bunny_rotate_3d( - mesh_density=32, + mesh_density=mesh_density, scale=1, final_time=final_time, simulation_config=simulation_config, ), ] + dataset = get_train_dataset( + dataset_type=config.td.dataset, config=config, device_count=1 + ) + # datasets = get_all_val_datasets(config=config, rank=0, world_size=1, device_count=1) + # dataset = datasets[1] + dataset.initialize_data() + dataloader = base_dataset.get_train_dataloader(dataset) + + # dataloader = None + + pca.run(dataloader, latent_dim=200, scenario=all_scenarios[0]) + simulation_runner.run_examples( all_scenarios=all_scenarios, file=__file__, plot_animation=True, config=Config(shell=False), + save_all=True, ) - pca.run(dataloader) - def get_train_dataset( dataset_type, config: TrainingConfig, rank: int = 0, world_size: int = 1, - device_count: int = 1, + device_count=None, item_fn=None, ): - device_count = get_device_count(config) + if device_count is None: + device_count = get_device_count(config) if dataset_type == "synthetic": train_dataset = SyntheticDataset( description="train", @@ -320,5 +332,4 @@ def main(args: Namespace): ) # Python 3.9+ args = parser.parse_args() # with jax.disable_jit(): - load_dotenv() main(args) diff --git a/deep_conmech/scene/scene_input.py b/deep_conmech/scene/scene_input.py index 66f68bc..75cb360 100644 --- a/deep_conmech/scene/scene_input.py +++ b/deep_conmech/scene/scene_input.py @@ -283,9 +283,7 @@ def get_target_data(self): self.norm_by_reduced_lifted_new_displacement ) # lower and rotate - target_data.last_displacement_step = thh.to_double( - self.get_last_displacement_step() - ) + target_data.new_displacement = thh.to_double(self.get_lifted_displacement()) ### skinning_acceleration = np.array( diff --git a/deep_conmech/training_config.py b/deep_conmech/training_config.py index 8431733..abd2f1b 100644 --- a/deep_conmech/training_config.py +++ b/deep_conmech/training_config.py @@ -118,7 +118,6 @@ def get_train_config(shell, mode): use_lhs_preconditioner=False, with_self_collisions=True, mesh_layer_proportion=4, # 2 4 - use_pca=False, mode=mode, ) return config diff --git a/examples/draft_compare.py b/examples/draft_compare.py index 962cc6a..f91a166 100644 --- a/examples/draft_compare.py +++ b/examples/draft_compare.py @@ -6,9 +6,8 @@ from tqdm import tqdm from conmech.helpers import cmh -from conmech.helpers.config import Config, SimulationConfig +from conmech.helpers.config import Config from conmech.scenarios import scenarios -from conmech.scenarios.scenarios import bunny_obstacles from conmech.simulations import simulation_runner from deep_conmech.graph.model_jax import RMSE from deep_conmech.run_model import get_newest_checkpoint_path @@ -17,15 +16,22 @@ def main(): load_dotenv() - cmh.print_jax_configuration() + modes = [ + "skinning_backwards", + "skinning", + "net", + "pca", + ] # Do not use "normal" - does not contain reduced + run_all_simulations(modes=modes) + compare_latest(modes=modes) + input("Press Enter to continue...") - # for mode in ["skinning_backwards", "skinning", "net"]: # Do not use normal - # run_simulation(mode) - training_config = TrainingConfig(shell=False) - checkpoint_path = get_newest_checkpoint_path(training_config) - compare_latest(checkpoint_path.split("/")[-1]) - input("Press Enter to continue...") +def run_all_simulations(modes): + cmh.print_jax_configuration() + + for mode in modes: + run_simulation(mode) def run_simulation(mode): @@ -35,7 +41,8 @@ def run_simulation(mode): config = get_train_config(shell=False, mode=mode) all_scenarios = scenarios.all_validation(config.td, config.sc) - all_scenarios = all_scenarios[-1] + all_scenarios = all_scenarios[-1] #-1] + # all_scenarios[0].schedule.final_time = 1. simulation_runner.run_examples( all_scenarios=all_scenarios, @@ -50,63 +57,61 @@ def get_error(simulation_1, simulation_2, index, key): return RMSE(simulation_1[index][key], simulation_2[index][key]) -def compare_latest(label=None): - current_time: str = datetime.now().strftime("%m.%d-%H.%M.%S") +def compare_latest(modes): + training_config = TrainingConfig(shell=False) + checkpoint_path = get_newest_checkpoint_path(training_config) input_path = "output" all_scene_files = cmh.find_files_by_extension(input_path, "scenes_comparer") + current_time: str = datetime.now().strftime("%m.%d-%H.%M.%S") + label = checkpoint_path.split("/")[-1] + main_path = f"output/{current_time}_COMPARE_{label}" dense = True path_id = "/scenarios/" if dense else "/scenarios_reduced/" scene_files = [f for f in all_scene_files if path_id in f] # all_arrays_path = max(scene_files, key=os.path.getctime) - normal = cmh.get_simulation(scene_files, "skinning_backwards") - skinning = cmh.get_simulation(scene_files, "skinning") - net = cmh.get_simulation(scene_files, "net") + def add_centered(simulation): + for step in simulation: + step["centered_new_displacement"] = step["new_displacement"] - step[ + "new_displacement" + ].mean(axis=0) + return simulation + + base = cmh.get_simulation(scene_files, "skinning_backwards") + # base = add_centered(base) - simulation_len = min(len(normal), len(skinning), len(net)) + cmh.create_folder(main_path) for key in [ "norm_lifted_new_displacement", "recentered_norm_lifted_new_displacement", - "displacement_old", - "exact_acceleration", "normalized_nodes", - "lifted_acceleration", + # "centered_new_displacement", + # "new_displacement", + # "displacement_old", + # "exact_acceleration", + # "normalized_nodes", + # # "lifted_acceleration", ]: - errors_skinning = [] - errors_net = [] - for index in tqdm(range(simulation_len)): - errors_skinning.append(get_error(skinning, normal, index, key)) - errors_net.append(get_error(net, normal, index, key)) - - print("Error net: ", np.mean(errors_net)) - print("Error skinning: ", np.mean(errors_skinning)) - all_errorrs = np.array([errors_skinning, errors_net]) - errors_df = pd.DataFrame(all_errorrs.T, columns=["skinning", "net"]) + errors_df = pd.DataFrame() + for mode in modes: + if mode == "skinning_backwards": + continue + pretendent = cmh.get_simulation(scene_files, mode) + # pretendent = add_centered(pretendent) + + simulation_len = min(len(base), len(pretendent)) + + errors = [] + for index in tqdm(range(simulation_len)): + errors.append(get_error(base, pretendent, index=index, key=key)) + errors_df[mode] = np.array(errors) + print(f"Error {mode} {key}: ", np.mean(errors)) + print() + plot = errors_df.plot() fig = plot.get_figure() - fig.savefig(f"output/{current_time}_{label}_dense:{dense}_{key}.png") - - #### - - # dense = False - # path_id = "/scenarios/" if dense else "/scenarios_reduced/" - # scene_files = [f for f in scene_files if path_id in f] - - # skinning = cmh.get_simulation(all_scene_files, 'skinning') - # net = cmh.get_simulation(all_scene_files, 'net') - - # simulation_len = min(len(normal), len(skinning), len(net)) - # for key in ['displacement_old', 'exact_acceleration', 'normalized_nodes']: - # errors_reduced = [] - # for index in tqdm(range(simulation_len)): - # errors_reduced.append(get_error(skinning, net, index, key)) - - # errors_df = pd.DataFrame(np.array([errors_reduced]).T, columns=['reduced']) - # plot = errors_df.plot() - # fig = plot.get_figure() - # fig.savefig(f"output/{current_time}_{label}_dense:{dense}_errors_{key}.png") - # print(errors_df) + fig.savefig(f"{main_path}/{dense}_{key}.png") if __name__ == "__main__": diff --git a/examples/draft_examples.py b/examples/draft_examples.py index bec39d2..86f0a25 100644 --- a/examples/draft_examples.py +++ b/examples/draft_examples.py @@ -32,7 +32,7 @@ def main(): load_dotenv() cmh.print_jax_configuration() - def get_simulation_config(mode="normal", use_pca=False): + def get_simulation_config(mode="normal"): return SimulationConfig( use_normalization=False, use_linear_solver=False, @@ -41,7 +41,6 @@ def get_simulation_config(mode="normal", use_pca=False): use_constant_contact_integral=False, use_lhs_preconditioner=False, with_self_collisions=True, - use_pca=use_pca, mesh_layer_proportion=4, mode=mode, ) diff --git a/examples/draft_presentation.py b/examples/draft_presentation.py index 5fccd5c..8d6e449 100644 --- a/examples/draft_presentation.py +++ b/examples/draft_presentation.py @@ -19,17 +19,16 @@ use_constant_contact_integral=False, use_lhs_preconditioner=False, with_self_collisions=False, - use_pca=False, ) -name = "linear_16" -simulation_config.use_green_strain = False -simulation_config.use_linear_solver = True -simulation_config.use_constant_contact_integral = True +name = "test_green" +# simulation_config.use_green_strain = False +# simulation_config.use_linear_solver = True +# simulation_config.use_constant_contact_integral = True animation_backend = "three" # None -def main(mesh_density=16, final_time=3, plot_animation=True): # 100 +def main(mesh_density=16, final_time=2, plot_animation=True): # 100 schedule = Schedule(final_time=final_time, time_step=0.01) body_prop = TimeDependentBodyProperties( mu=8, diff --git a/examples/examples_2d.py b/examples/examples_2d.py index a9e2792..7e24ace 100644 --- a/examples/examples_2d.py +++ b/examples/examples_2d.py @@ -18,7 +18,6 @@ use_constant_contact_integral=False, use_lhs_preconditioner=False, with_self_collisions=False, - use_pca=False, ) @@ -48,68 +47,6 @@ def main(mesh_density=20, final_time=5, plot_animation=True): # 40 ), ] all_scenarios = [ - # scenarios.polygon_mesh_obstacles( - # mesh_density=mesh_density, - # scale=1, - # final_time=final_time, - # simulation_config=simulation_config, - # ), - # Scenario( - # name="circle_slide_roll", - # mesh_prop=MeshProperties( - # dimension=2, - # mesh_type=scenarios.M_CIRCLE, - # scale=[1], - # mesh_density=[mesh_density], - # ), - # body_prop=scenarios.default_body_prop, - # schedule=schedule, - # forces_function=np.array([0.0, -0.5]), - # obstacle=obstacles[0], - # simulation_config=simulation_config, - # ), - # Scenario( - # name="circle_flat_A_roll", - # mesh_prop=MeshProperties( - # dimension=2, - # mesh_type=scenarios.M_CIRCLE, - # scale=[1], - # mesh_density=[mesh_density], - # ), - # body_prop=scenarios.default_body_prop, - # schedule=schedule, - # forces_function=np.array([2.0, -0.5]), - # obstacle=obstacles[1], - # simulation_config=simulation_config, - # ), - # Scenario( - # name="circle_flat_B_roll", - # mesh_prop=MeshProperties( - # dimension=2, - # mesh_type=scenarios.M_CIRCLE, - # scale=[1], - # mesh_density=[mesh_density], - # ), - # body_prop=scenarios.default_body_prop, - # schedule=schedule, - # forces_function=np.array([2.0, -0.5]), - # obstacle=obstacles[2], - # simulation_config=simulation_config, - # ), - # Scenario( - # name="circle_flat_C_roll", - # mesh_prop=MeshProperties( - # dimension=2, - # mesh_type=scenarios.M_CIRCLE, - # scale=[1], - # mesh_density=[mesh_density], - # ), - # body_prop=scenarios.default_body_prop, - # schedule=schedule, - # forces_function=np.array([2.0, -0.5]), - # obstacle=obstacles[3], - # simulation_config=simulation_config, - # ), Scenario( name="rectangle_flat_roll", mesh_prop=MeshProperties( @@ -124,6 +61,68 @@ def main(mesh_density=20, final_time=5, plot_animation=True): # 40 obstacle=obstacles[4], simulation_config=simulation_config, ), + scenarios.polygon_mesh_obstacles( + mesh_density=mesh_density, + scale=1, + final_time=final_time, + simulation_config=simulation_config, + ), + Scenario( + name="circle_slide_roll", + mesh_prop=MeshProperties( + dimension=2, + mesh_type=scenarios.M_CIRCLE, + scale=[1], + mesh_density=[mesh_density], + ), + body_prop=scenarios.default_body_prop, + schedule=schedule, + forces_function=np.array([0.0, -0.5]), + obstacle=obstacles[0], + simulation_config=simulation_config, + ), + Scenario( + name="circle_flat_A_roll", + mesh_prop=MeshProperties( + dimension=2, + mesh_type=scenarios.M_CIRCLE, + scale=[1], + mesh_density=[mesh_density], + ), + body_prop=scenarios.default_body_prop, + schedule=schedule, + forces_function=np.array([2.0, -0.5]), + obstacle=obstacles[1], + simulation_config=simulation_config, + ), + Scenario( + name="circle_flat_B_roll", + mesh_prop=MeshProperties( + dimension=2, + mesh_type=scenarios.M_CIRCLE, + scale=[1], + mesh_density=[mesh_density], + ), + body_prop=scenarios.default_body_prop, + schedule=schedule, + forces_function=np.array([2.0, -0.5]), + obstacle=obstacles[2], + simulation_config=simulation_config, + ), + Scenario( + name="circle_flat_C_roll", + mesh_prop=MeshProperties( + dimension=2, + mesh_type=scenarios.M_CIRCLE, + scale=[1], + mesh_density=[mesh_density], + ), + body_prop=scenarios.default_body_prop, + schedule=schedule, + forces_function=np.array([2.0, -0.5]), + obstacle=obstacles[3], + simulation_config=simulation_config, + ), ] return simulation_runner.run_examples( diff --git a/examples/examples_3d.py b/examples/examples_3d.py index 647d739..33a552b 100644 --- a/examples/examples_3d.py +++ b/examples/examples_3d.py @@ -27,7 +27,6 @@ use_constant_contact_integral=False, use_lhs_preconditioner=False, with_self_collisions=False, - use_pca=False, ) diff --git a/examples/examples_temperature_2d.py b/examples/examples_temperature_2d.py index fafe48f..2dfe619 100644 --- a/examples/examples_temperature_2d.py +++ b/examples/examples_temperature_2d.py @@ -19,14 +19,13 @@ from conmech.state.obstacle import Obstacle simulation_config = SimulationConfig( - use_normalization=False, + use_normalization=True, use_linear_solver=False, use_green_strain=False, use_nonconvex_friction_law=False, use_constant_contact_integral=False, use_lhs_preconditioner=False, with_self_collisions=False, - use_pca=False, ) diff --git a/examples/examples_temperature_3d.py b/examples/examples_temperature_3d.py index b3bbf08..80aa1e4 100644 --- a/examples/examples_temperature_3d.py +++ b/examples/examples_temperature_3d.py @@ -30,7 +30,6 @@ use_constant_contact_integral=False, use_lhs_preconditioner=False, with_self_collisions=False, - use_pca=False, ) diff --git a/examples/temperature_2023.py b/examples/temperature_2023.py index e71a80d..4924bda 100644 --- a/examples/temperature_2023.py +++ b/examples/temperature_2023.py @@ -30,7 +30,6 @@ def main(mesh_density=32, final_time=2.5, plot_animation=True, shell=False): use_constant_contact_integral=False, use_lhs_preconditioner=False, with_self_collisions=False, - use_pca=False, ) schedule = Schedule(final_time=final_time, time_step=0.01) diff --git a/scripts/examples.sh b/scripts/examples.sh index 4462c91..b7bde9e 100755 --- a/scripts/examples.sh +++ b/scripts/examples.sh @@ -2,4 +2,4 @@ # Using .venv screen -XS examples quit -screen -dmS examples bash -c 'cd ..; source .venv/bin/activate; clear; PYTHONPATH=. python examples/draft_examples.py --shell; exec bash' +screen -dmS examples bash -c 'cd ..; source .venv/bin/activate; clear; PYTHONPATH=. python examples/draft_compare.py --shell; exec bash' diff --git a/tests/test_conmech/regression/test_schedule_temperature_2d.py b/tests/test_conmech/regression/test_schedule_temperature_2d.py index 40af56e..83b4303 100644 --- a/tests/test_conmech/regression/test_schedule_temperature_2d.py +++ b/tests/test_conmech/regression/test_schedule_temperature_2d.py @@ -16,7 +16,6 @@ use_constant_contact_integral=False, use_lhs_preconditioner=False, with_self_collisions=False, - use_pca=False, ) diff --git a/tests/test_conmech/regression/test_shedule_roll_2d.py b/tests/test_conmech/regression/test_shedule_roll_2d.py index 7b349d5..6f53fd3 100644 --- a/tests/test_conmech/regression/test_shedule_roll_2d.py +++ b/tests/test_conmech/regression/test_shedule_roll_2d.py @@ -16,7 +16,6 @@ use_constant_contact_integral=False, use_lhs_preconditioner=False, with_self_collisions=False, - use_pca=False, ) diff --git a/tests/test_conmech/smoke/test_examples.py b/tests/test_conmech/smoke/test_examples.py index df95bb4..1d6dd99 100644 --- a/tests/test_conmech/smoke/test_examples.py +++ b/tests/test_conmech/smoke/test_examples.py @@ -26,7 +26,6 @@ use_constant_contact_integral=False, use_lhs_preconditioner=False, with_self_collisions=False, - use_pca=False, ) default_args = dict(show=False, save=False)