CODEWITHSUNDEEP
pythonarraysnumpy
8556732
8556732

Reputation: 311

Python, numpy and arrays - array multiplication/division in my function

I've edited this post to make things simpler, but the original is below.

I basically want to do some math on an array which has the same number of columns but varying rows for a series of inputs. So do the calculation using an input array of columns (1 row), where the calculation is applied on a row by row basis for the larger array. The larger array has the same number of columns as the input. The number of columns will vary depending on what calculation I am doing, but the two inputs will always have the same number of columns; for example array 1 is (2000L,3L), array 2 would be (1L,3L) and if array 1 is (3050L,5L) array 2 (1L,5L) and so on.

Here's the example I am trying at the moment: 

# Make data (for 3 components):
K1, K2, K3 = 30., 16., 36.

G1, G2, G3 = 12., 15., 44.

G = np.array([[G1, G2, G3]])
K = np.array([[K1, K2, K3]])
# Fake data:
fracs = np.array([[0.3,0.6,0.1],[0.5,0.3,0.2],[0.6,0.0,0.4],[0.2,0.8,0.0],[0.1,0.3,0.6]])

def reuss(fracs, K, G):
    rK = 1.0 / (np.sum(fracs/K))
    rG = 1.0 / (np.sum(fracs/G))
    return rK, rG
# Do the function
Kr, Gr = reuss(fracs, K, G)

Which gives me: 

In[]: fracs
Out[]: 
array([[ 0.3,  0.6,  0.1],
       [ 0.5,  0.3,  0.2],
       [ 0.6,  0. ,  0.4],
       [ 0.2,  0.8,  0. ],
       [ 0.1,  0.3,  0.6]])
In[]: K, G
Out[]: (array([[  30., 16., 36.]]), array([[  12., 15., 44.]]))
In[]: Kr, Gr
Out[]: (nan, nan)

I'm aware you can broadcast scalars, but what I want to do is have a single function that will let me broadcast an array of inputs.

ORIGINAL POST HERE FOR COMPLETENESS:

So I'm trying to improve my existing function that implements the Reuss average of a mixture of minerals in a rock. For now I've been writing this very very long windedly for a system of nfractions and nmineral phases. This is my pasted code (which is very bad):

# Reuss isostress bound

def reuss3(K1, G1, K2=None, G2=None, K3=None, G3=None, PHI, f1=1.0, f2=None, f3=None, Kf, Gf=None):
'''
Calculates the Reuss lower bound of Bulk and Shear Moduli. 
Also known as the isostress bound. Represents a suspension of grains.
Default assumes 1 rock phase, this function can be expanded to 3. 

Form taken from Mavko et.al (2009) p.174 

Arguments:
K# = Effective Bulk modulus of phase # in GPa
G# = Effective Shear modulus of phase # in GPa
PHI = Porosity fraction of mixture
f# = Fraction of phase #, default is 100% of phase 1
Kf, Gf = Moduli of fluid, default assumes fluid has 0 shear modulus

Results: 
Kr = Reuss bound of Bulk modulus in GPa
Gr = Reuss bound of Shear modulus in GPa
'''
Kr = 1. /( (((1.-PHI)/K1)*f1) + (((1.-PHI)/K2)*f2) + (((1.-PHI)/K3)*f3) + (PHI/Kf) )
Gr = 1. /( (((1.-PHI)/G1)*f1) + (((1.-PHI)/G2)*f2) + (((1.-PHI)/G3)*f3) + (PHI/Gf) )
return Kr, Gr

I've been trying to improve this, and actually forward model the Reuss average of real rocks from their known composition, porosity and water saturation. My inputs for this are a series of mineral fractions that sum to 1 sampled in the depth dimension, so basically a 2D array with (nsamples, nPhases). I've no idea how to take an array of mineral properties in a 1D array of (nphases) and do the calculation for each dept sample. So far I've tried this, which outputs nan's:

# Make data (for 3 components):
K1, K2, K3 = 30., 16., 36.

G1, G2, G3 = 12., 15., 44.

G = np.array([[G1, G2, G3]])
K = np.array([[K1, K2, K3]])
# Fake data:
fracs = np.array([[0.3,0.6,0.1],[0.5,0.3,0.2],[0.6,0.0,0.4],[0.2,0.8,0.0],[0.1,0.3,0.6]])

def reuss(fracs, K, G):
    rK = 1.0 / (np.sum(fracs/K))
    rG = 1.0 / (np.sum(fracs/G))
    return rK, rG
# Do the function
rK, rG = reuss(fracs, K, G)

I imagine there's something wrong with the shape of the arrays, though I need some advice on how to move forward with this.

Upvotes: 0

Views: 178

Answers (2)

kevinkayaks
kevinkayaks

Reputation: 2726

I interpreted your question as the number of rows is variable, so that your output array should be a 1D array with the length the same as the number of rows. You're essentially doing a dot product between each row fracs[i] (3d) and the vector 1/K ( or 1/G) (3d) for each of the 5 rows, then doing 5 sums. Is this correct?

If so, I think your code as written will sum over the wrong axis. It should be this: 

def reuss(fracs, K, G): # You want 5 outputs for rK and rG, right? 
    rK = 1 / np.sum(fracs / K, axis = 1) # I interpret this as what you want. 
    rG = 1 / np.sum(fracs / G, axis = 1)
    return rK, rG

# alternatively you can stack K and G and do it all at once 
rK,rG = (1/((fracs[...,np.newaxis]/np.stack((K,G),-1)[np.newaxis,...]).sum(axis=1))).T

Either way returns

>>> reuss(fracs,K,G)
(array([19.88950276, 24.40677966, 32.14285714, 17.64705882, 25.80645161]),
 array([14.86486486, 15.10297483, 16.92307692, 14.28571429, 23.8267148 ]))

edit: this is with no extra axis inside your definition of K and G. You write np.array([[K1,K2,K3]]). I write np.array([K1,K2,K3]).

Upvotes: 0

steph
steph

Reputation: 86

When I run both your snippets, it works properly. Can you try it again? If it still does not work, give us your python and numpy versions.

Side note: the only piece of code that DOES NOT WORK is the one the you say you use, because PHI has no default and follows arguments with defaults.

Upvotes: 1

Related Questions