chmarti1
Array Write Protected Error
Ok, I warn up front that this one is probably going to be very hard to track down.
pyromat version = bleeding edge np version = 2.1.0 python version = 3.12.3
Steps to reproduce:
from pyromat.igtools import IgtMix
air = IgtMix('0.01289563 Ar + 0.00047711 CO2 + 0.75520558 N2 + 0.23142168 O2')
air.p(p=1, T=300)
Result:
ValueError: output array is read-only
For what it's worth, it's being caused by the array of p being marked as read only, and you're trying to edit its value in place with np.multiply(). If you at an intermediate point within igtmix.p() re-assign the value with p = p.copy() it resolves the issue, but that doesn't seem ideal. I'm not up on the specifics of numpy data typing sufficiently to be able to trace this much further. If I have a chance to test in a different virtual environment with different numpy version, etc. I will do so.
Edit, I'm also able to achieve it on the d function call with the same setup as before and:
air.d(p=1, d=1.1)
Ok I did a little messing around on a different computer, previous version of this message was in error because I had been using the wrong branch of the pyromat code.
in python 3.10.4 numpy versions 1.26.4 - code sample above fails numpy version 2.1.0 - code sample above fails numpy version 2.2.0 - code sample above fails
in python 3.12.8 numpy version 2.2.0 - code sample above fails
I'm not totally clear on whether this is something that's caused by pyromat or the underlying numpy compiled code. Searching for the error suggests that it has to do with when an array occupies a non-contiguous area in the memory. So it's possible that in one step along the way, the inplace calls are resulting in something affecting the memory.
Great find! This is definitely not the first time I've seen this error. It came up quite a bit in development, and I thought I had chased down all of its gremlins! In most of my use of Numpy, I try to do in-place operations wherever I can, but for precisely this reason, that can be difficult. For example, if the user passes an array that's a view of an array, things can start breaking.
You're right, the simplest answer is usually making copies, but as a rule, I try very hard not to do that for large-array efficiency. The compromise is often to disallow in-place operations on anything I didn't initialize myself, so I can be certain that it isn't a view or some other problematic instance.
This issue convinces me that it is time to comb through the IgtMix again looking for consistency. Frankly, array broadcasting with the mixtures was the absolute biggest headache to implementing IgtMix, so it isn't shocking that there may still be some improvements to be made.
If I change the temperature units, the same bug happens on T as well. However, the problem is isolated to T(), p(), and d().
The source of the problem comes from the way IgtMix._argparse() handles broadcasting versus the way ig._argparse() or ig2._argparse() handles it. In addition to the normal broadcasting rules, IgtMix instances also need to handle broadcasting with a (potentially multidimensional) mixture array. The result is that the IgtMix class has a method called _broadcast_shape that determines the dimensions of the final property result, and all intermediates are broadcast into that shape. By contrast, the other classes just use Numpy's broadcast_arrays(), which automatically detects array shape and "just handles it."
The Problem To do this job, IgtMix._argparse() uses Numpy's broadcast_to() function, which lets users specify the broadcast shape explicitly. Just like the reshape() function, the broadcast_to() function returns a view that is read-only. That was a deliberate choice to avoid huge redundant copies of simple data, which can be a real problem in these high dimensional arrays.
This behavior can be overwritten by setting p.flags.writeable = True, but writing to a broadcast view can give unexpected results that will frustrate most users. This should not be visible to the user.
The Solution If T, p, and d are only used as intermediate values in calculations, then this works precisely as intended. It is efficient, elegant, and compact. We don't want to abandon that behavior.
Your solution of making a copy is actually ideal.
if p.base is not None:
p = np.copy(p)
pm.units.pressure(p, from_units='Pa', inplace=True)
return p
The ensures that only when pressure is being passed to the user, it will be promoted to a fully populated array of the correct shape.