Repository: IndEcol/pymrio

Issue:
downstream sales coefficient calculations (and Ghosh)
@Beckebanze 
It seems that the Ghosh implementation in pymrio is transposed to normal? As is transposed (G is calculated correctly). As should be sales coefficients – so along a row, how much of the sale of commodity A to total sales of commodity A. If a fix is generated, great if you can use conventional nomenclature (B instead of As).
 
See p543 and so on of miller and blair

sorry, formatting of py code is off...
`
import pymrio 
import numpy as np

# load the test system
io = pymrio.load_test()

# calculate the basic parts of the system
io.calc_all()

# we need value added, but the test system doesn't have an explicit value added - so we generate it as the different between total output(=total input) and intermediate inputs
io.v=io.x.T - io.Z.sum(axis=0)

# standard Ghosh
x=io.v@io.G
print(x.T-io.x)
# alt standard method is
x=io.G.T@io.v.T
print(x-io.x)

# why transpose? this works
x=io.v@io.G.T
print(x.T-io.x)

# why transpose? this works
x=io.G@io.v.T
print(x-io.x)
`

Hints from issue comments:
@Beckebanze: Happy to change As to B in my next pull request (it's in fact already on my to-do).

Regarding the definition of the Ghosh inverse: I am aware that it is the transpose compared to the definition in several papers in literature, including Miller and Blair. For value-added (v) and final-demand (y) methods, it is convenient to work with L.dot(y) and v.dot(G). However, for the combined upstream and downstream footprint of some economic activity it is somewhat cumbersome to have the upstream footprint vector vertically and downstream horizontally. The total footprint(s) of satellite account(s) F of any economic activity expressed in terms of vertical revenue-vector r can now be expressed as F.dot(L+G-I).dot(r). 
Personally, I prefer to call the L and G in this expression in 'Leontief and Ghosh inverses'. Using the transpose of this for G makes the expression less elegant, because then the transpose of the Ghosh inverse has to be used at all times. If you still think it is better to implement the transpose of the current implemented G, then happy to do so as well with my next PR.
@rich-wood: Hi, I'd say it is best to keep consistency to established methods. It gets too confusing if there are slightly different implementations. Not everyone will use the sales coefficients and Ghosh in the way you do. For your implementation, it is still pretty straightforward, as you just need to transpose G when using it: F.dot(L+G.T-I).dot(r). That'd have to be easier right?@Beckebanze: I changed the implementation for Ghosh accordingly with PR https://github.com/IndEcol/pymrio/pull/146. @konstantinstadler: Thanks a lot for that, I will release that soon.
Just a heads up: I will make the gosh a optional calculation. It doubles the storage requirment and most people are not using it (AFAIK)

Solution patch:
diff --git a/doc/source/math.rst b/doc/source/math.rst
index e9532e5..0c6e2c2 100644
--- a/doc/source/math.rst
+++ b/doc/source/math.rst
@@ -181,13 +181,15 @@ Upstream and Downstream scope 3
 ------------
 
 In the context of impact analyses the factors of production are often categorized into scope 1, 2 and 3, with scope 3
+
 sub-divided into upstream and downstream.
 For a MRIO the scope 1 is the direct impact of the industries. The factors of production scope 1 associated
 with some product or service in sector 'i' of monetary value 'r' is given by :math:`S e_i r',
 where :math:`e_i' is the 'i^{th}' unit vector.
-Scope 2 is the indirect impact through directly consumed energy (mostly electricity). The precise defintion of scope 2 in an MRIO depends
- on the list of MRIO sectors that are classified as scope 2 energy suppliers. Scope 2 is therefore included in the
-upstream scope 3, which we refer to as upstream indirect impacts. The upstream multipliers are
+Scope 2 is the indirect impact through directly consumed energy (mostly electricity). The precise definition of scope 2
+in an MRIO depends on the list of MRIO sectors that are classified as scope 2 energy suppliers. Scope 2 is therefore
+included in the upstream scope 3, which we refer to as upstream indirect impacts. The upstream multipliers are
+
 .. math::
 
     \begin{equation}
@@ -201,25 +203,28 @@ The downstream attribution according to Ghosh is done by the input share matrix
 .. math::
 
     \begin{equation}
-        A^{*} = Z^{T} \hat{x}^{-1}
+        B = \hat{x}^{-1} Z
     \end{equation}
 
 Note that we have defined this matrix in analogy with :math:`A`, meaning that the factors of production coefficient
-are applied from the right-hand side. The full downstream multiplier (including scope 1) is given
-by :math:`S G` where
+are given by :math:`S G^{T}` where
 
 .. math::
 
     \begin{equation}
-        G = (\mathrm{I}- Z^{T}\hat{x}^{-1})^{-1}
+        G = (\mathrm{I}- B)^{-1} = \hat{x} L \hat{x}^{-1}
     \end{equation}
-is the transpose of the Ghosh inverse matrix.
+
+is the Ghosh inverse matrix.
 The pure downstream multiplier (excluding scope 1) is given by
 
+
 .. math::
 
     \begin{equation}
-        M_{down} = S((\mathrm{I}- Z^{T}\hat{x}^{-1})^{-1} -I) = S(G-I) = S(\hat{x}^{-1} L^{T} \hat{x} - I)
+
+        M_{down} = S(G^{T}-I)
+
     \end{equation}
 
 The sector's total impact multiplier is simply the sum of :math:`M_{up}`, :math:`S` and :math:`M_{down}`.
diff --git a/pymrio/core/mriosystem.py b/pymrio/core/mriosystem.py
index bf83f40..8492021 100644
--- a/pymrio/core/mriosystem.py
+++ b/pymrio/core/mriosystem.py
@@ -31,7 +31,7 @@ from pymrio.core.constants import (
 from pymrio.tools.iomath import (
     calc_A,
     calc_accounts,
-    calc_As,
+    calc_B,
     calc_F,
     calc_F_Y,
     calc_G,
@@ -2144,7 +2144,7 @@ class IOSystem(BaseSystem):
         Z=None,
         Y=None,
         A=None,
-        As=None,
+        B=None,
         x=None,
         L=None,
         G=None,
@@ -2164,7 +2164,7 @@ class IOSystem(BaseSystem):
         self.Y = Y
         self.x = x
         self.A = A
-        self.As = As
+        self.B = B
         self.L = L
         self.G = G
         self.unit = unit
@@ -2294,8 +2294,8 @@ class IOSystem(BaseSystem):
             self.A = calc_A(self.Z, self.x)
             self.meta._add_modify("Coefficient matrix A calculated")
 
-        if self.As is None:
-            self.As = calc_As(self.Z, self.x)
+        if self.B is None:
+            self.B = calc_B(self.Z, self.x)
             self.meta._add_modify("Coefficient matrix As calculated")
 
         if self.L is None:
@@ -2303,7 +2303,7 @@ class IOSystem(BaseSystem):
             self.meta._add_modify("Leontief matrix L calculated")
 
         if self.G is None:
-            self.G = calc_G(self.As)
+            self.G = calc_G(self.B)
             self.meta._add_modify("Ghosh matrix G calculated")
 
         return self
diff --git a/pymrio/tools/iomath.py b/pymrio/tools/iomath.py
index 2a41b58..daca2c8 100644
--- a/pymrio/tools/iomath.py
+++ b/pymrio/tools/iomath.py
@@ -152,10 +152,10 @@ def calc_A(Z, x):
         return Z * recix
 
 
-def calc_As(Z, x):
-    """Calculate the As matrix (coefficients) from Z and x
+def calc_B(Z, x):
+    """Calculate the B matrix (coefficients) from Z and x
 
-    As is a normalized version of the industrial flows of the transpose of Z, which quantifies the input to downstream
+    B is a normalized version of the industrial flows of the transpose of Z, which quantifies the input to downstream
     sectors.
 
     Parameters
@@ -187,13 +187,15 @@ def calc_As(Z, x):
         recix = recix.reshape((1, -1))
     # use numpy broadcasting - factor ten faster
     # Mathematical form - slow
-    # return Z.dot(np.diagflat(recix))
+    # return np.diagflat(recix).dot(Z)
     if type(Z) is pd.DataFrame:
         return pd.DataFrame(
-            np.transpose(Z.values) * recix, index=Z.index, columns=Z.columns
+            np.transpose(np.transpose(Z.values) * recix),
+            index=Z.index,
+            columns=Z.columns,
         )
     else:
-        return np.transpose(Z) * recix
+        return np.transpose(np.transpose(Z) * recix)
 
 
 def calc_L(A):
@@ -232,30 +234,26 @@ def calc_L(A):
         return np.linalg.inv(I - A)
 
 
-def calc_G(As, L=None, x=None):
-    """Calculate the Ghosh inverse matrix G either from As (high computation effor) or from Leontief matrix L and x
+def calc_G(B, L=None, x=None):
+    """Calculate the Ghosh inverse matrix G either from B (high computation effort) or from Leontief matrix L and x
     (low computation effort)
 
-    G = inverse matrix of (I - As) = hat(x)^{-1} *  L^T * hat(x)
+    G = inverse matrix of (I - B) = hat(x) *  L * hat(x)^{-1}
 
-    where I is an identity matrix of same shape as As, and hat(x) is the diagonal matrix with values of x on the
+    where I is an identity matrix of same shape as B, and hat(x) is the diagonal matrix with values of x on the
     diagonal.
 
-    Note that we define G as the transpose of the Ghosh inverse matrix, so that we can apply the factors of
-    production intensities from the left-hand-side for both Leontief and Ghosh attribution. In this way the
-    multipliers have the same (vector) dimensions and can be added.
-
     Parameters
     ----------
-    As : pandas.DataFrame or numpy.array
+    B : pandas.DataFrame or numpy.array
         Symmetric input output table (coefficients)
 
     Returns
     -------
     pandas.DataFrame or numpy.array
         Ghosh input output table G
-        The type is determined by the type of As.
-        If DataFrame index/columns as As
+        The type is determined by the type of B.
+        If DataFrame index/columns as B
 
     """
     # if L has already been calculated, then G can be derived from it with low computational cost.
@@ -275,20 +273,20 @@ def calc_G(As, L=None, x=None):
 
         if type(L) is pd.DataFrame:
             return pd.DataFrame(
-                recix * np.transpose(L.values * x), index=Z.index, columns=Z.columns
+                np.transpose(recix * np.transpose(L.values * x)),
+                index=Z.index,
+                columns=Z.columns,
             )
         else:
-            # G = hat(x)^{-1} *  L^T * hat(x) in mathematical form np.linalg.inv(hatx).dot(L.transpose()).dot(hatx).
+            # G = hat(x) *  L * hat(x)^{-1} in mathematical form hatx.dot(L.transpose()).dot(np.linalg.inv(hatx)).
             # it is computationally much faster to multiply element-wise because hatx is a diagonal matrix.
-            return recix * np.transpose(L * x)
+            return np.transpose(recix * np.transpose(L * x))
     else:  # calculation of the inverse of I-As has a high computational cost.
-        I = np.eye(As.shape[0])  # noqa
-        if type(As) is pd.DataFrame:
-            return pd.DataFrame(
-                np.linalg.inv(I - As), index=As.index, columns=As.columns
-            )
+        I = np.eye(B.shape[0])  # noqa
+        if type(B) is pd.DataFrame:
+            return pd.DataFrame(np.linalg.inv(I - B), index=B.index, columns=B.columns)
         else:
-            return np.linalg.inv(I - As)  # G = inverse matrix of (I - As)
+            return np.linalg.inv(I - B)  # G = inverse matrix of (I - B)
 
 
 def calc_S(F, x):
@@ -421,7 +419,7 @@ def calc_M(S, L):
 def calc_M_down(S, G):
     """Calculate downstream multipliers of the extensions
 
-    M_down = S * ( G - I )
+    M_down = S * ( G^T - I )
 
     Where I is an identity matrix of same shape as G
 
@@ -440,7 +438,7 @@ def calc_M_down(S, G):
         If DataFrame index/columns as S
 
     """
-    return S.dot(G - np.eye(G.shape[0]))
+    return S.dot(np.transpose(G) - np.eye(G.shape[0]))
 
 
 def calc_e(M, Y):


Test patch:
diff --git a/tests/test_math.py b/tests/test_math.py
index 48a6a5e..f92c3e2 100644
--- a/tests/test_math.py
+++ b/tests/test_math.py
@@ -16,7 +16,7 @@ sys.path.append(os.path.join(TESTPATH, ".."))
 # the function which should be tested here
 from pymrio.tools.iomath import calc_A  # noqa
 from pymrio.tools.iomath import calc_accounts  # noqa
-from pymrio.tools.iomath import calc_As  # noqa
+from pymrio.tools.iomath import calc_B  # noqa
 from pymrio.tools.iomath import calc_e  # noqa
 from pymrio.tools.iomath import calc_F  # noqa
 from pymrio.tools.iomath import calc_F_Y  # noqa
@@ -235,21 +235,28 @@ def td_small_MRIO():
             columns=_Z_multiindex,
         )
 
-        As = pd.DataFrame(
+        B = pd.DataFrame(
             data=[
-                [0.19607843, 0.0, 0.17241379, 0.1, 0.0, 0.12820513],  # noqa
-                [0.09803922, 0.13333333, 0.05172414, 0.0, 0.19230769, 0.0],  # noqa
-                [0.01960784, 0.0, 0.34482759, 0.0, 0.01923077, 0.0],  # noqa
-                [0.11764706, 0.0, 0.06896552, 0.02, 0.0, 0.02564103],  # noqa
                 [
+                    0.19607843,
                     0.09803922,
-                    0.33333333,
+                    0.01960784,
+                    0.11764706,
+                    0.09803922,
+                    0.1372549,
+                ],  # noqa
+                [0.0, 0.13333333, 0.0, 0.0, 0.33333333, 0.2],  # noqa
+                [
+                    0.17241379,
+                    0.05172414,
+                    0.34482759,
+                    0.06896552,
                     0.03448276,
-                    0.2,
-                    0.38461538,
-                    0.25641026,
+                    0.0,
                 ],  # noqa
-                [0.1372549, 0.2, 0.0, 0.18, 0.01923077, 0.25641026],  # noqa
+                [0.1, 0.0, 0.0, 0.02, 0.2, 0.18],  # noqa
+                [0.0, 0.19230769, 0.01923077, 0.0, 0.38461538, 0.01923077],  # noqa
+                [0.12820513, 0.0, 0.0, 0.02564103, 0.25641026, 0.25641026],
             ],
             index=_Z_multiindex,
             columns=_Z_multiindex,
@@ -314,52 +321,52 @@ def td_small_MRIO():
             data=[
                 [
                     1.33871463,
-                    0.07482651,
-                    0.38065717,
-                    0.19175489,
-                    0.04316348,
-                    0.25230906,
+                    0.28499301,
+                    0.05727924,
+                    0.1747153,
+                    0.58648041,
+                    0.38121958,
                 ],  # noqa
                 [
-                    0.28499301,
+                    0.07482651,
                     1.37164729,
-                    0.22311379,
-                    0.15732254,
-                    0.44208094,
-                    0.20700334,
+                    0.02969614,
+                    0.02185805,
+                    0.93542302,
+                    0.41222074,
                 ],  # noqa
                 [
-                    0.05727924,
-                    0.02969614,
+                    0.38065717,
+                    0.22311379,
                     1.54929428,
-                    0.02349465,
-                    0.05866155,
-                    0.03091401,
+                    0.15941084,
+                    0.39473223,
+                    0.17907022,
                 ],  # noqa
                 [
-                    0.1747153,
-                    0.02185805,
-                    0.15941084,
+                    0.19175489,
+                    0.15732254,
+                    0.02349465,
                     1.05420074,
-                    0.01404088,
-                    0.07131676,
+                    0.60492361,
+                    0.34854312,
                 ],  # noqa
                 [
-                    0.58648041,
-                    0.93542302,
-                    0.39473223,
-                    0.60492361,
+                    0.04316348,
+                    0.44208094,
+                    0.05866155,
+                    0.01404088,
                     1.95452858,
-                    0.79595212,
+                    0.18081884,
                 ],  # noqa
                 [
-                    0.38121958,
-                    0.41222074,
-                    0.17907022,
-                    0.34854312,
-                    0.18081884,
+                    0.25230906,
+                    0.20700334,
+                    0.03091401,
+                    0.07131676,
+                    0.79595212,
                     1.48492515,
-                ],  # noqa
+                ],
             ],
             index=_Z_multiindex,
             columns=_Z_multiindex,
@@ -637,15 +644,15 @@ def test_calc_A_MRIO(td_small_MRIO):
     pdt.assert_frame_equal(td_small_MRIO.A, calc_A(td_small_MRIO.Z, x_Tvalues))
 
 
-def test_calc_As_MRIO(td_small_MRIO):
-    pdt.assert_frame_equal(td_small_MRIO.As, calc_As(td_small_MRIO.Z, td_small_MRIO.x))
+def test_calc_B_MRIO(td_small_MRIO):
+    pdt.assert_frame_equal(td_small_MRIO.B, calc_B(td_small_MRIO.Z, td_small_MRIO.x))
     # we also test the different methods to provide x:
     x_values = td_small_MRIO.x.values
     x_Tvalues = td_small_MRIO.x.T.values
     x_series = pd.Series(td_small_MRIO.x.iloc[:, 0])
-    pdt.assert_frame_equal(td_small_MRIO.As, calc_As(td_small_MRIO.Z, x_series))
-    pdt.assert_frame_equal(td_small_MRIO.As, calc_As(td_small_MRIO.Z, x_values))
-    pdt.assert_frame_equal(td_small_MRIO.As, calc_As(td_small_MRIO.Z, x_Tvalues))
+    pdt.assert_frame_equal(td_small_MRIO.B, calc_B(td_small_MRIO.Z, x_series))
+    pdt.assert_frame_equal(td_small_MRIO.B, calc_B(td_small_MRIO.Z, x_values))
+    pdt.assert_frame_equal(td_small_MRIO.B, calc_B(td_small_MRIO.Z, x_Tvalues))
 
 
 def test_calc_Z_MRIO(td_small_MRIO):
@@ -664,7 +671,7 @@ def test_calc_L_MRIO(td_small_MRIO):
 
 
 def test_calc_G_MRIO(td_small_MRIO):
-    pdt.assert_frame_equal(td_small_MRIO.G, calc_G(td_small_MRIO.As))
+    pdt.assert_frame_equal(td_small_MRIO.G, calc_G(td_small_MRIO.B))
 
 
 def test_calc_S_MRIO(td_small_MRIO):
diff --git a/tests/test_parsers.py b/tests/test_parsers.py
index 9d70fb5..224830c 100644
--- a/tests/test_parsers.py
+++ b/tests/test_parsers.py
@@ -138,10 +138,14 @@ def test_parse_exio_ext():
         col1, index_col=1, name="col1", sep="\t", drop_compartment=False
     )
 
-    pdt.assert_frame_equal(ext1.F, ext2.F)
-    pdt.assert_frame_equal(ext1.F, ext3.F)
-    pdt.assert_frame_equal(ext1.F, ext5.F)
-    pdt.assert_frame_equal(ext1.F, ext5_wcmp.F.droplevel("compartment"))
+
+    assert np.allclose(ext5.F, ext3.F)
+    assert np.allclose(ext3.F, ext2.F)
+    assert np.allclose(ext2.F, ext1.F)
+    assert np.allclose(ext1.F, ext5_wcmp.F)
+
+    # pdt.assert_frame_equal(ext1.F, ext5_wcmp.F.droplevel("compartment"))
+
 
     assert ext5.F.iloc[1, 1] == 102
     assert ext5.F.iloc[-1, -1] == 148