joheinze / PriNCe

Propagation including Nuclear Cascade equations
BSD 3-Clause "New" or "Revised" License
10 stars 3 forks source link

Idea how to make the cross section swapping for PANDORA fast #19

Open afedynitch opened 4 years ago

afedynitch commented 4 years ago

I'm having an idea how to make a simple replacement of individual cross sections possible without putting all the time in the optimizations.

Roughly the idea would be that single cross sections would be replaced in the completed batch matrix. Of course, new channels can not be added because of system size and sparsity pattern, but this is not what they intend to do in the first place.

joheinze commented 4 years ago

Yes, I think this will work and would be easy to implement.

The internal data structure consists just of a IntRate._batch_matrix with shape $(E_{cr} \cdot nchannels ) \times n{ph}$ and row and column vectors of size $(E_{cr} \cdot n_channels)$. The row and column vectors uniquely map the elements to the interactions matrix. So all that needs to be done is use this mapping to replace the right elements in self_batch_matrix. The only minor complication here is that due to the CSC format of the sparse matrix, the different energy bins of one interaction channel are not necessarily adjacent.

I am pretty sure, that it will even be easy to add additional channels.

One just needs to append them to IntRate._batch_matrix, IntRate._batch_rows, IntRate._batch_cols and then resort them to create a new coupling matrix. The functions for this are already in IntRate._init_coupling_mat, since IntRate._init_matrices does not created the batch vectors in the correct order. Sorting takes less than 2s on my machine, so performance is not an issue here

afedynitch commented 4 years ago

Actually I even have a simpler idea:

If one wants to swap a cross section, then one can take a cross section. Calculate the f and g functions "manually" and then swap the elements in the batch matrix. This may be even done without modifying the package on the level of a notebook or so. I don't see any evident limitation to this approach.

joheinze commented 4 years ago

Yes, that is basically what I meant in my last comment. One only need to overwrite IntRate._batch_matrix and IntRate._batch_rows and IntRate._batch_cols. This can in principle done in a Notebook if one really know what they are doing, but that would be more of a hack. The intention was however for IntRate._batch_matrix to be a private attribute (therefore the underscore). The complication is that the CSR format requires the the batch vector (and by extension the batch matrix) to be odered by row and then by column. Since the interaction rates are on the diagonal of the interaction matrix, the elements of one mother-daughter channel will not be adjacent in the batch vector. So one really needs to understand the structure, which is not very user friendly...

I would suggest to add simple getter and setter functions for the batch matrix, which take care of the mapping. If only the values in the batch matrix are changed, this will be enough.

One could even add new channels by apending to IntRate._batch_matrix and IntRate._batch_rows and IntRate._batch_cols (or delete them by removing elements). Then the coupling matrix in CSR format has to be resorted. probably it is enought to call IntRate._init_coupling_mat again. The sorting is done in these lines: https://github.com/joheinze/PriNCe/blob/862bec85354186151a04cad76a7916f895eafcf8/prince/interaction_rates.py#L308-L321

afedynitch commented 4 years ago

Haha true, yes. It's the same :) So with what you're saying one doesn't need to rely on the same sparsity structure if the conversion to CSR/CSC is repeated.

The only thing I added was that one has to simply do the averaging eps_r -> y externally. So some part of the cross section class needs to be exposed or implemented in external functions.

joheinze commented 4 years ago

Yes, true just repeat the conversion to CSR. This should be no problem since the batch vector will already be mostly sorted, so np.lexsort should be fast.

True, we basically need only need to factorize these two parts from _init_matrices into separate function. (This would also improve the readability of _init_matrices anyway.)

https://github.com/joheinze/PriNCe/blob/68b809a58ad931206dc2dffb67707967b21ced46/prince_cr/interaction_rates.py#L166-L190

https://github.com/joheinze/PriNCe/blob/68b809a58ad931206dc2dffb67707967b21ced46/prince_cr/interaction_rates.py#L205-L262