added function IA3

source/python/mods/disaggregation.py

@@ -39,6 +39,7 @@
 # @Module Content:
 #    - dapoly
 #    - darain
+#    - IA3
 #
 #*******************************************************************************
 
@@ -139,3 +140,263 @@ def darain(alist):
     nfield = xac + xbd
 
     return nfield
+
+def IA3(g):
+    """
+
+    ***************************************************************************
+    * Copyright 2017                                                          *
+    * Sabine Hittmeir, Anne Philipp, Petra Seibert                            *
+    *                                                                         *
+    * This work is licensed under the Creative Commons Attribution 4.0        *
+    * International License. To view a copy of this license, visit            *
+    * http://creativecommons.org/licenses/by/4.0/ or send a letter to         *
+    * Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.            *
+    ***************************************************************************
+
+    @Description
+       The given data series will be interpolated with a non-negative geometric
+       mean based algorithm. The original grid is reconstructed by adding two
+       sampling points in each data series interval. This subgrid is used to
+       keep all information during the interpolation within the associated
+       interval. Additionally, an advanced monotonicity filter is applied to
+       improve the monotonicity properties of the series.
+
+       For more information see article:
+       Hittmeir, S.; Philipp, A.; Seibert, P. (2017): A conservative
+       interpolation scheme for extensive quantities with application to the
+       Lagrangian particle dispersion model FLEXPART.,
+       Geoscientific Model Development
+
+    @Input
+       g: list of float values
+          A list of float values which represents the complete data series that
+          will be interpolated having the dimension of the original raw series.
+
+    @Return
+       f: list of float values
+          The interpolated data series with additional subgrid points.
+          Its dimension is equal to the length of the input data series
+          times three.
+    """
+
+    #######################  variable description #############################
+    #                                                                         #
+    # i      - index variable for looping over the data series                #
+    # g      - input data series                                              #
+    # f      - interpolated and filtered data series with additional          #
+    #          grid points                                                    #
+    # fi     - function value at position i, f_i                              #
+    # fi1    - first  sub-grid function value f_i^1                           #
+    # fi2    - second sub-grid function value f_i^2                           #
+    # fip1   - next function value at position i+1, f_(i+1)                   #
+    # dt     - time step                                                      #
+    # fmon   - monotonicity filter                                            #
+    #                                                                         #
+    ###########################################################################
+
+
+    import numpy as np
+
+    # time step
+    dt=1.0
+
+    ############### Non-negative Geometric Mean Based Algorithm ###############
+
+    # for the left boundary the following boundary condition is valid:
+    # the value at t=0 of the interpolation algorithm coincides with the
+    # first data value according to the persistence hypothesis
+    f=[g[0]]
+
+    # compute two first sub-grid intervals without monotonicity check
+    # go through the data series and extend each interval by two sub-grid
+    # points and interpolate the corresponding data values
+    # except for the last interval due to boundary conditions
+    for i in range(0,2):
+
+        # as a requirement:
+        # if there is a zero data value such that g[i]=0, then the whole
+        # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0
+        # according to Eq. (6)
+        if g[i]==0.:
+            f.extend([0.,0.,0.])
+
+        # otherwise the sub-grid values are calculated and added to the list
+        else:
+            # temporal save of last value in interpolated list
+            # since it is the left boundary and hence the new (fi) value
+            fi = f[-1]
+
+            # the value at the end of the interval (fip1) is prescribed by the
+            # geometric mean, restricted such that non-negativity is guaranteed
+            # according to Eq. (25)
+            fip1=min( 3.*g[i] , 3.*g[i+1] , np.sqrt(g[i+1]*g[i]) )
+
+            # the function value at the first sub-grid point (fi1) is determined
+            # according to the equal area condition with Eq. (19)
+            fi1=3./2.*g[i]-5./12.*fip1-1./12.*fi
+
+            # the function value at the second sub-grid point (fi2) is determined
+            # according Eq. (18)
+            fi2=fi1+1./3.*(fip1-fi)
+
+            # add next interval of interpolated (sub-)grid values
+            f.append(fi1)
+            f.append(fi2)
+            f.append(fip1)
+
+    # compute rest of the data series intervals
+    # go through the data series and extend each interval by two sub-grid
+    # points and interpolate the corresponding data values
+    # except for the last interval due to boundary conditions
+    for i in range(2,len(g)-1):
+
+        # as a requirement:
+        # if there is a zero data value such that g[i]=0, then the whole
+        # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0
+        # according to Eq. (6)
+        if g[i]==0.:
+            # apply monotonicity filter for interval before
+            # check if there is "M" or "W" shape
+            if     np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5])==-1 \
+               and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4])==-1 \
+               and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3])==-1:
+
+                # the monotonicity filter corrects the value at (fim1) by
+                # substituting (fim1) with (fmon), see Eq. (27), (28) and (29)
+                fmon = min(3.*g[i-2], \
+                           3.*g[i-1], \
+                           np.sqrt(max(0,(18./13.*g[i-2] - 5./13.*f[-7]) *
+                                         (18./13.*g[i-1] - 5./13.*f[-1]))))
+
+                # recomputation of the sub-grid interval values while the
+                # interval boundaries (fi) and (fip2) remains unchanged
+                # see Eq. (18) and (19)
+                f[-4]=fmon
+                f[-6]=3./2.*g[i-2]-5./12.*fmon-1./12.*f[-7]
+                f[-5]=f[-6]+(fmon-f[-7])/3.
+                f[-3]=3./2.*g[i-1]-5./12.*f[-1]-1./12.*fmon
+                f[-2]=f[-3]+(f[-1]-fmon)/3.
+
+            f.extend([0.,0.,0.])
+
+        # otherwise the sub-grid values are calculated and added to the list
+        else:
+            # temporal save of last value in interpolated list
+            # since it is the left boundary and hence the new (fi) value
+            fi = f[-1]
+
+            # the value at the end of the interval (fip1) is prescribed by the
+            # geometric mean, restricted such that non-negativity is guaranteed
+            # according to Eq. (25)
+            fip1=min( 3.*g[i] , 3.*g[i+1] , np.sqrt(g[i+1]*g[i]) )
+
+            # the function value at the first sub-grid point (fi1) is determined
+            # according to the equal area condition with Eq. (19)
+            fi1=3./2.*g[i]-5./12.*fip1-1./12.*fi
+
+            # the function value at the second sub-grid point (fi2) is determined
+            # according Eq. (18)
+            fi2=fi1+1./3.*(fip1-fi)
+
+            # apply monotonicity filter for interval before
+            # check if there is "M" or "W" shape
+            if     np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5])==-1 \
+               and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4])==-1 \
+               and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3])==-1:
+
+                # the monotonicity filter corrects the value at (fim1) by
+                # substituting (fim1) with fmon, see Eq. (27), (28) and (29)
+                fmon = min(3.*g[i-2], \
+                           3.*g[i-1], \
+                           np.sqrt(max(0,(18./13.*g[i-2] - 5./13.*f[-7]) *
+                                         (18./13.*g[i-1] - 5./13.*f[-1]))))
+
+                # recomputation of the sub-grid interval values while the
+                # interval boundaries (fi) and (fip2) remains unchanged
+                # see Eq. (18) and (19)
+                f[-4]=fmon
+                f[-6]=3./2.*g[i-2]-5./12.*fmon-1./12.*f[-7]
+                f[-5]=f[-6]+(fmon-f[-7])/3.
+                f[-3]=3./2.*g[i-1]-5./12.*f[-1]-1./12.*fmon
+                f[-2]=f[-3]+(f[-1]-fmon)/3.
+
+            # add next interval of interpolated (sub-)grid values
+            f.append(fi1)
+            f.append(fi2)
+            f.append(fip1)
+
+    # separate treatment of the final interval
+
+    # as a requirement:
+    # if there is a zero data value such that g[i]=0, then the whole
+    # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0
+    # according to Eq. (6)
+    if g[-1]==0.:
+        # apply monotonicity filter for interval before
+        # check if there is "M" or "W" shape
+        if     np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5])==-1 \
+           and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4])==-1 \
+           and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3])==-1:
+
+            # the monotonicity filter corrects the value at (fim1) by
+            # substituting (fim1) with (fmon), see Eq. (27), (28) and (29)
+            fmon = min(3.*g[-3], \
+                       3.*g[-2], \
+                       np.sqrt(max(0,(18./13.*g[-3] - 5./13.*f[-7]) *
+                                     (18./13.*g[-2] - 5./13.*f[-1]))))
+
+            # recomputation of the sub-grid interval values while the
+            # interval boundaries (fi) and (fip2) remains unchanged
+            # see Eq. (18) and (19)
+            f[-4]=fmon
+            f[-6]=3./2.*g[-3]-5./12.*fmon-1./12.*f[-7]
+            f[-5]=f[-6]+(fmon-f[-7])/3.
+            f[-3]=3./2.*g[-2]-5./12.*f[-1]-1./12.*fmon
+            f[-2]=f[-3]+(f[-1]-fmon)/3.
+
+        f.extend([0.,0.,0.])
+
+    # otherwise the sub-grid values are calculated and added to the list
+    # using the persistence hypothesis as boundary condition
+    else:
+        # temporal save of last value in interpolated list
+        # since it is the left boundary and hence the new (fi) value
+        fi = f[-1]
+        # since last interval in series, last value is also fip1
+        fip1 = g[-1]
+        # the function value at the first sub-grid point (fi1) is determined
+        # according to the equal area condition with Eq. (19)
+        fi1 = 3./2.*g[-1]-5./12.*fip1-1./12.*fi
+        # the function value at the second sub-grid point (fi2) is determined
+        # according Eq. (18)
+        fi2 = fi1+dt/3.*(fip1-fi)
+
+        # apply monotonicity filter for interval before
+        # check if there is "M" or "W" shape
+        if     np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5])==-1 \
+           and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4])==-1 \
+           and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3])==-1:
+
+            # the monotonicity filter corrects the value at (fim1) by
+            # substituting (fim1) with (fmon), see Eq. (27), (28) and (29)
+            fmon = min(3.*g[-3], \
+                       3.*g[-2], \
+                       np.sqrt(max(0,(18./13.*g[-3] - 5./13.*f[-7]) *
+                                     (18./13.*g[-2] - 5./13.*f[-1]))))
+
+            # recomputation of the sub-grid interval values while the
+            # interval boundaries (fi) and (fip2) remains unchanged
+            # see Eq. (18) and (19)
+            f[-4]=fmon
+            f[-6]=3./2.*g[-3]-5./12.*fmon-1./12.*f[-7]
+            f[-5]=f[-6]+(fmon-f[-7])/3.
+            f[-3]=3./2.*g[-2]-5./12.*f[-1]-1./12.*fmon
+            f[-2]=f[-3]+(f[-1]-fmon)/3.
+
+        # add next interval of interpolated (sub-)grid values
+        f.append(fi1)
+        f.append(fi2)
+        f.append(fip1)
+
+    return f