Bugfix Model class incomaptible with inhomogeneous arrays

[bch0w]

Related to #214, tested and contributed by @casarotti

Changelog:

• Model class now saves internal model representation as 'object' arrays to deal with chunks of differing lengths
• Model .npz files are now saved as pickle files
• Update and improve tests to cover Model class

[raulleoncz]

Hellow @bch0w,

Could this bug be related to the mpirun issue #194 using specfem2d?

[bch0w]

@raulleoncz It may be! Your original error message in #194 matches the one seen in #214. I was unaware of this issue as I was always using uniformly distributed meshes. Are you able to re-test your problem with the latest devel branch of the code to see if you still have any issues?

[raulleoncz]

Hello @bch0w, thanks for the answer.

I tried to simulate my problem using the new version of seisflows and this is what I got. First of all, I tried using python 3.12 and numpy 1.26.4 and I got the following error:

2024-06-05 17:37:53 [INFO] | checking true/target model parameters: 2024-06-05 17:37:53 [DEBU] | checkpointing workflow to seisflows state file Traceback (most recent call last): File "/Users/raul/anaconda3/envs/seisflows/bin/seisflows", line 8, in sys.exit(main()) File "/Users/raul/seisflows/seisflows/seisflows.py", line 1456, in main sf() File "/Users/raul/seisflows/seisflows/seisflows.py", line 456, in call getattr(self, self._args.command)(vars(self._args)) File "/Users/raul/seisflows/seisflows/seisflows.py", line 813, in submit system.submit(workdir=self._args.workdir, File "/Users/raul/seisflows/seisflows/system/workstation.py", line 213, in submit workflow.run() File "/Users/raul/seisflows/seisflows/workflow/inversion.py", line 239, in run super().run() # Runs task list File "/Users/raul/seisflows/seisflows/workflow/forward.py", line 283, in run func() File "/Users/raul/seisflows/seisflows/workflow/inversion.py", line 281, in generate_synthetic_data super().generate_synthetic_data(kwargs) File "/Users/raul/seisflows/seisflows/workflow/forward.py", line 313, in generate_synthetic_data self.solver.check_model_values(path=self.path.model_true) File "/Users/raul/seisflows/seisflows/solver/specfem.py", line 329, in check_model_values _model.check() File "/Users/raul/seisflows/seisflows/tools/model.py", line 402, in check assert(val.any()), ( File "/Users/raul/anaconda3/envs/seisflows/lib/python3.10/site-packages/numpy/core/_methods.py", line 57, in _any return umr_any(a, axis, dtype, out, keepdims) ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

At first, I thought this problem was related to the numpy version but after I downgraded it to numpy=1.23.5 the above error persisted. The original code in model.py was this:

        assert(val.any()), (
                f"SPECFEM_{self.flavor} model '{key}' has no values, please "
                f"check your input model `path_model_init` and the chosen "
                f"`material` which controls the expected parameters"
                )

I tried using val.all() instead of val.any() but the error was still there. Since 'val' is a numpy array with more than one element, I thought we were not working with this variable correctly. After trying several ideas, I thought we could evaluate each of the individual arrays in the val list using:

for key, val in self.model.items():

        # Modified
        print("val:", val) # Just to see if val was an array with more than one element
        for array in val:
            assert array.all(), (
        # Original part
        # assert(val.any()), (
                f"SPECFEM_{self.flavor} model '{key}' has no values, please "
                f"check your input model `path_model_init` and the chosen "
                f"`material` which controls the expected parameters"
                )

        # Make sure none of the values are NaNs
        assert(not np.isnan(np.hstack(val).astype(float)).any()), (
             f"SPECFEM_{self.flavor} model '{key}' contains NaN values and "
             f"should not, please check your model construction"
             )

    # Check the physicality of the parameters
    if self.flavor in ["2D", "3D"]:
        self._check_2d3d_parameters(min_pr, max_pr)
    elif self.flavor == "3DGLOBE":
        self._check_3dglobe_parameters(min_pr, max_pr)

The above modification worked for my particular problem. I haven't analyzed if this can affect other parts of the workflow, but for now I was able to perform the simulation using NPROC > 1 when each processor has a variable workload.

Do you think this change can affect other parts of the seisflows workflow? If not, do you think it could be a solution?

UPDATE

1. After finalizing the evaluate_initial_misfit tasks, I got the following message:

/Users/raul/seisflows/seisflows/tools/model.py:362: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray. model[key] = np.array(model[key])

This warning is because I'm working with lists of different sizes. As far as I know, older version of numpy had no problem dealing with this but now this warning appears. Maybe modifying model[key] = np.array(model[key]) to model[key] = np.array(model[key], dtype=object) could work.

2. After 3 iterations I got the following error before evaluating line search misfit:

Traceback (most recent call last): File "/Users/raul/anaconda3/envs/seisflows/bin/seisflows", line 8, in sys.exit(main()) File "/Users/raul/seisflows/seisflows/seisflows.py", line 1456, in main sf() File "/Users/raul/seisflows/seisflows/seisflows.py", line 456, in call getattr(self, self._args.command)(**vars(self._args)) File "/Users/raul/seisflows/seisflows/seisflows.py", line 813, in submit system.submit(workdir=self._args.workdir, File "/Users/raul/seisflows/seisflows/system/workstation.py", line 213, in submit workflow.run() File "/Users/raul/seisflows/seisflows/workflow/inversion.py", line 239, in run super().run() # Runs task list File "/Users/raul/seisflows/seisflows/workflow/forward.py", line 283, in run func() File "/Users/raul/seisflows/seisflows/workflow/inversion.py", line 561, in evaluate_line_search_misfit self.solver.check_model_values(path=os.path.join(self.path.eval_func, File "/Users/raul/seisflows/seisflows/solver/specfem.py", line 329, in check_model_values _model.check() File "/Users/raul/seisflows/seisflows/tools/model.py", line 419, in check self._check_2d3d_parameters(min_pr, max_pr) File "/Users/raul/seisflows/seisflows/tools/model.py", line 443, in _check_2d3d_parameters logger.warning(f"Vp minimum is negative {self.model.vp.min()}") File "/Users/raul/anaconda3/envs/seisflows/lib/python3.10/site-packages/numpy/core/_methods.py", line 44, in _amin return umr_minimum(a, axis, None, out, keepdims, initial, where) ValueError: operands could not be broadcast together with shapes (149750,) (202375,)

Maybe a reshaping function could avoid this problem.

[bch0w]

Hi @raulleoncz, thanks for the detailed error description. I'll share some thoughts here but would you mind opening this as a new issue as well so that we can continue the discussion outside of the closed PR? Thanks!

1. It might help if you are able to share your model file with me, I'm not sure how big the files are but I think it would help a lot in diagnosing the issue. It seems like there is still some lingering issues with non-uniform model arrays that are not being accounted for.
2. Looks like I missed adding an 'object' dtype to get around this ragged array here (https://github.com/adjtomo/seisflows/blob/devel/seisflows/tools/model.py#L362). That should be a simple fix, thanks for catching it, I will try and get to that asap.
3. In your last issue it looks like a warning is preceding the actual error ( logger.warning(f"Vp minimum is negative {self.model.vp.min()}")), which suggests your model is becoming unphysical (negative velocities). I can fix the ValueError but you will likely still get errors from having unphysical values.