Merge pull request #1 from mpimd-csc/prep-for-zenodo

Prep for zenodo

.gitignore

@@ -3,6 +3,9 @@
 *.xml.gz
 _tmp_*_chunk*
 *.json
+*.tikz
+*.xopp
+scripts/mh-syncit.sh
 
 pics/
 logs/
@@ -12,6 +15,9 @@ mechthild-logs/
 gen_pod_uq.egg-info/
 __pycache__
 scripts/job-files
+scripts/cached-data
+scripts/mh-data
+scripts/deprecated
 
 ChangeLog
 build/

README.md

@@ -8,7 +8,7 @@ This is the code of the numerical experiments in our paper
 > Benner, Heiland (2020): *Space and Chaos-expansion Galerkin POD Low-order
 > Discretization of PDEs for Uncertainty Quantification*
 
-in the second version from July 2022.
+in the third version from December 2022.
 
 ## Installation
 
@@ -21,16 +21,16 @@ pip install -e .  # make sure you use Python 3
 ```
 
 if the installation of `multim-galerkin-pod` fails because of `scikit-sparse`
-try `pip install --no-deps multidim-galerkin-pod==1.0.2` instead. 
+try `pip install --no-deps multidim-galerkin-pod==1.1.0` instead. 
 
 The source are in `gen_pod_uq` and the files for the simulations are in `scripts`.
 
 ## Rerun the simulations
 
-**NOTE**: For reproduction of the results, use version `1.1.1` of the package to be installed like
+**NOTE**: For reproduction of the results, use version `1.1.4` of the package to be installed like
 
 ```sh
-pip install gen-pod-uq==1.1.1
+pip install gen-pod-uq==1.1.4
 ```
 
 from the [`pypi repo`](https://pypi.org/project/gen-pod-uq/)
@@ -53,6 +53,8 @@ source runitall.sh
 
 You may want to comment out some parts.
 
+### Post Processing
+
 The raw data of our simulations is provided in the folder `rawdata`. In order to postprocess copy it to the `scripts/cached-data` folder
 
 ```sh
@@ -61,10 +63,29 @@ The raw data of our simulations is provided in the folder `rawdata`. In order to
 # ## caution: this may replace computed data
 ```
 
-### Post Processing
 
 ```
 cd scripts
 source postprocess.sh
 ```
 
+### Evaluating the Kolmogorov Metric
+
+```sh
+# ## caution: this may replace computed data
+# cp rawdata/ysoltens-for-kolmogorov-metric-evaluation/*npy scripts/cached-data/
+# ## caution: this may replace computed data
+```
+
+```sh
+cd scripts
+python3 kolmogorov-metrix.py
+```
+
+In order to (only) compute the plots, one may run a reduced experiment by setting 
+
+```py
+onlyplots = True
+```
+
+in `kolmogorov-metrix.py`.

gen_pod_uq/mc_pce_utils.py

@@ -40,16 +40,16 @@ def sfvariance(para):
     return sfvariance
 
 
-def get_nu_sample(distribution='uniform', uncdims=1, nulb=0, nuub=1):
+def get_nu_sample(distribution='uniform', uncdims=1, nulb=0., nuub=1.):
     if distribution == 'uniform':
         def nusample(nsamples):
             return nulb + (nuub-nulb)*np.random.rand(nsamples, uncdims)
     else:
         dstp = distribution.split('-')
         wghtfnc = dstp[0]
         if wghtfnc == 'beta':
-            alpha = np.float(dstp[1])
-            beta = np.float(dstp[2])
+            alpha = float(dstp[1])
+            beta = float(dstp[2])
 
             def nusample(nsamples):
                 return nulb + \
@@ -219,8 +219,8 @@ def setup_pce(distribution='uniform', distrpars={}, pcedim=None, uncdims=None):
         dstp = distribution.split('-')
         wghtfnc = dstp[0]
         if wghtfnc == 'beta':
-            alpha = np.float(dstp[1])
-            beta = np.float(dstp[2])
+            alpha = float(dstp[1])
+            beta = float(dstp[2])
             abscissae, weights = ceu.\
                 get_weighted_gaussqr(N=pcedim, weightfunction=wghtfnc,
                                      wfpdict=dict(alpha=alpha, beta=beta),
@@ -229,16 +229,33 @@ def setup_pce(distribution='uniform', distrpars={}, pcedim=None, uncdims=None):
             raise NotImplementedError()
 
     scalefac = (1./weights.sum())**uncdims
+    # this seems to be right
+
+    wprod = weights
+    for _ in range(uncdims-1):
+        wprod = np.kron(wprod, weights)
+
+    # scalefactwo = (1./wprod.sum())
+    # if not np.allclose(scalefac, scalefactwo):
+    #     raise UserWarning('need to check this again')
 
     def comp_expv(ytens):
         ydim = ytens.shape[0]
-        expv = 0
         maty = ytens.reshape((ydim, -1))
-    # for kk, cw in enumerate(weights):
+        expvone = wprod @ maty.T
+
+        expvtwo = 0
         for idx, wtpl in enumerate(product(weights, repeat=uncdims)):
             cw = (np.array(wtpl)).prod()
-            expv += cw*maty[:, idx]
-        return scalefac*expv
+            expvtwo += cw*maty[:, idx]
+            # print(wtpl)
+            # print(idx)
+
+        if ydim == 1:
+            logging.info(f'diff in expvs {expvone.item()-expvtwo.item():.4e}')
+        else:
+            pass
+        return scalefac*expvone
 
     def comp_vrnc(ytens, expv):
         ydim = ytens.shape[0]
@@ -263,3 +280,38 @@ def pce_comp_vrnc(ytens, expv, weights=None, uncdims=None, scalefac=None):
         cw = (np.array(wtpl)).prod()
         vrnc += cw*matysqrd[:, idx]
     return scalefac*vrnc - expv**2
+
+
+def eva_all_lgrngs(abscissae, x):
+    '''by ChatGPT
+    '''
+    n = len(abscissae)
+    polynomials = []
+    for i in range(n):
+        # Initialize the Lagrange polynomial as L_i(x) = 1
+        polynomial = 1
+        for j in range(n):
+            if i == j:
+                # Skip the term corresponding to the current abscissa
+                continue
+            polynomial *= (x - abscissae[j]) / (abscissae[i] - abscissae[j])
+        polynomials.append(polynomial)
+    return polynomials
+
+
+def empirical_cdf(samples):
+    ''' by ChatGPT'''
+    # Sort the samples in ascending order
+    sorted_samples = sorted(samples)
+    N = len(samples)
+
+    # Create a list to store the CDF values
+    cdf_values = []
+
+    # Iterate over the sorted samples
+    for i, _ in enumerate(sorted_samples):
+        # Compute the CDF value for the current sample
+        cdf_value = (i + 1) / N
+        cdf_values.append(cdf_value)
+
+    return cdf_values