Multistatereporter variable pos frequency#767
Conversation
…rent frequency to energies data.
currently no pos test is failing...
a0e7dbb to
a262c2b
Compare
|
This PR will take over #712. |
Accessing data from netCDF returns a masked array, even if no data is masked. For valid data, extract it and give back just the basic numpy array. For invalid data, give back a nan array with the proper shape.
Codecov ReportAttention: Patch coverage is
🚀 New features to boost your workflow:
|
|
@ijpulidos @IAlibay @mikemhenry I can't request reviewers, but I think this is ready for review. |
ijpulidos
left a comment
There was a problem hiding this comment.
This looks good to me. Thanks a lot! Great work. Just a couple of non-blocking comments.
| raise IndexError | ||
| # pull the valid data out of masked array | ||
| positions = unit.Quantity(x.data, unit.nanometers) | ||
| except (IndexError, KeyError): |
There was a problem hiding this comment.
Non-blocking comment. Should we try outputting something here when an error is caught? Maybe something warning users that there's missing data or something similar. What are the reasons for the array to be masked?
There was a problem hiding this comment.
On the topic of how MaskedArrays work here. Assuming that a particular coordinate exists in a dimension (e.g. at least some data exists for a given iteration) we always get back a MaskedArray. The MaskedArray will have a data attribute with a numpy array that will be masked by the MaskedArrays mask attribute. If the data was all valid, then we just pull that information out and expose it to the user. If no data exists at that dimension coordinate, then we'll just get an IndexError so we give back a zero array in place of the nonsense data inside the MaskedArray. By doing it this way, a user never has to deal with a MaskedArray.
I think that if the data was masked it's basically the same behavior as if the data didn't exist and we would have expected a zero array anyways. I don't really see much of a reason to give a warning about this since having a giant array of zeros is kind of a warning already. I'd even be careful saying the data is missing since that would imply it should have been there in the first place.
| information, 0 would prevent information being written | ||
| velocity_interval : int, default 1 | ||
| the frequency at which to write positions relative to analysis | ||
| information, 0 would prevent information being written |
There was a problem hiding this comment.
Non-blocking question. This might be too late now to ask, but why is that we want positions and velocities to be stored at different intervals? Does it make sense to define a state (i.e. with positions and velocities) from different intervals? Just thinking about ways of not having to deal with the "complexity" of having to check for pos and vel intervals independently.
There was a problem hiding this comment.
I think this would depend on the intent of the simulation. I could see that someone interested in doing VAC calculations wouldn't care as much about the positions.
There was a problem hiding this comment.
So the one case I can think of is that you might not want to store the velocities but you do want to store the positions.
I might be wrong, but as far as I remember the "state" is stored in the checkpoint? So you could recover a simulation from the checkpoint, and just store positions in the simulation file.
|
I'll be testing this in a Protocol before I review :) |
| analysis_particle_indices : tuple of ints, Optional. Default: () (empty tuple) | ||
| If specified, it will serialize positions and velocities for the specified particles, at every iteration, in the | ||
| reporter storage (.nc) file. If empty, no positions or velocities will be stored in this file for any atoms. | ||
| position_interval : int, default 1 |
There was a problem hiding this comment.
Could we store these values inside of the NetCDF file as additional variables? This would make it a lot easier to know what the spacing is between frames, rather than having to iterate through the positions.
There was a problem hiding this comment.
I'm being blind sorry 🙈 - it's below
| if 'box_vectors' in storage.variables: | ||
| # Restore box vectors. | ||
| x = storage.variables['box_vectors'][read_iteration, replica_index, :, :].astype(np.float64) | ||
| # TODO: Are box vectors also variably saved? |
There was a problem hiding this comment.
From wwhat I can tell, they seem to be saved at the same rate as positions - worth double checking though!
There was a problem hiding this comment.
It is saved at the same rate of positions assuming that the sampler_state's box_vectors is not None (line 1669). I find this line a little odd though since there is ais_periodic property of the reporter (line 204) that follows slightly different logic. I'm not sure which is best to use. @ijpulidos do you have any thoughts on this (I've pinged you directly at the offending line)?
There was a problem hiding this comment.
To me it makes sense that these are stored at the same rate of positions, but again, there might be a use case where we want different rates for it? The logic here being that to maintain the thermodynamic states (especially pressure) you need to vary the box size and that's dependent on the positions. I'd argue that having positions without box sizes, or vice versa, is probably something we don't want to deal with.
Try and retrieve values from the netcdf file. If the data is masked, assume it hasn't been written at that iteration and raise an `IndexError`. If the box_vectors variable doesn't exist, then a `KeyError` is raised. Both execptions will bind None to the `box_vectors` variable.
| storage = self._storage_dict[storage_file] | ||
| # Check if the schema must be initialized, do this regardless of the checkpoint_interval for consistency | ||
| is_periodic = True if (sampler_states[0].box_vectors is not None) else False | ||
| is_periodic = sampler_states[0].box_vectors is not None |
There was a problem hiding this comment.
@ijpulidos is there a reason we aren't using the MultiStateReporter.is_periodic property here?
There was a problem hiding this comment.
I cannot say for sure. I guess this means that one can have a box_vectors variable in the reporter but not necessarily have box_vectors attribute in the sampler states, and vice versa. I don't know if we ever want those NOT to match. Is there a use case for such scenario?
There was a problem hiding this comment.
I suppose I'd be inclined to keep it the same in this PR, it just struck me as strange. Both lines were introduced in the same commit so I'm going to go on the assumption it was deliberate rather than an oversight.
|
@IAlibay @ijpulidos unless there is more to do here, I think this can be merged! |
IAlibay
left a comment
There was a problem hiding this comment.
I can't approve on this repo, but lgtm!
Continuation of #712.
There are multiple avenues to address the current failing tests. I would like to address those concerns here before continuing in that PR.