Merge pull request #29 from sandialabs/25-uqpc-workflow

25 uqpc workflow

Author: ksargsyan

apps/awkies/getrange.x

@@ -0,0 +1,48 @@
+#!/bin/bash
+#=====================================================================================
+
+# Given a Nxd matrix of samples, compute dx2 ranges
+# Example: getrange.x samples.dat > ranges.dat
+
+shopt -s expand_aliases
+alias awk="awk -v OFMT='%.15e'"
+
+if [ $# -lt 1 ]; then
+    echo "Number of arguments can not be less than 1"
+    echo "Syntax: $0 <filename> [<cushion_fraction>]"
+    exit
+elif [ $# -gt 2 ]; then
+    echo "Number of arguments can not be greater than 2"
+    echo "Syntax: $0 <filename> [<cushion_fraction>]"
+    exit
+elif [ $# -eq 1 ]; then
+    fr=0.0
+elif [ $# -eq 2 ]; then
+    fr=$2
+fi
+
+filename=$1
+
+DIM=`awk 'END{print NF}' $filename`
+
+for (( COL=1; COL<=$DIM ; COL++ )); do
+
+awk 'BEGIN {
+valmin=1e+100;
+valmax=-1e+100;
+line=1;
+}
+{
+if( $(col) < valmin )
+{
+  valmin=$(col);
+}
+if( $(col) > valmax )
+{
+  valmax=$(col);
+}
+}
+END{
+print valmin-fr*(valmax-valmin), valmax+fr*(valmax-valmin)
+}' col=$COL $filename
+done

apps/awkies/scale.x

@@ -0,0 +1,35 @@
+#!/bin/bash -e
+# Scales matrix data to or from a given parameter domain to [-1,1]^d
+# scale.x <input> <to or from> <domain> <output>
+
+shopt -s expand_aliases
+alias awk="awk -v OFMT='%.15e'"
+
+IN_FILE=$1
+TO_FROM=$2
+DOM_FILE=$3
+OUT_FILE=$4
+
+N=`awk 'END{print NR}' $IN_FILE`
+D=`awk 'END{print NR}' $DOM_FILE`
+
+DD=`awk 'END{print NF}' $IN_FILE`
+
+## check that DD=D
+
+echo "" > $OUT_FILE
+for (( i=1; i<=$D ; i++ )); do
+A=`awk 'NR==i{print $1}' i=$i $DOM_FILE`
+B=`awk 'NR==i{print $2}' i=$i $DOM_FILE`
+
+if [ "$TO_FROM" = "from" ]; then
+    awk '{print -1.+2.*($i-a)/(b-a)}' i=$i a=$A b=$B $IN_FILE > tmp
+elif [ "$TO_FROM" = "to" ]; then
+    awk '{print (a+b)/2.+$i*(b-a)/2.}' i=$i a=$A b=$B $IN_FILE > tmp
+else
+    echo "Second argument has to be 'to' or 'from'"
+    exit
+fi
+paste $OUT_FILE tmp > tmpp; mv tmpp $OUT_FILE
+
+done

apps/awkies/transpose.x

@@ -0,0 +1,26 @@
+#!/bin/bash -e
+# Transpose a matrix: assumes all lines have same number
+# of fields
+
+# Example: transpose_file.x file_in > file_out
+
+shopt -s expand_aliases
+alias awk="awk -v OFMT='%.15e'"
+
+exec awk '
+NR == 1 {
+	n = NF
+	for (i = 1; i <= NF; i++)
+		row[i] = $i
+	next
+}
+{
+	if (NF > n)
+		n = NF
+	for (i = 1; i <= NF; i++)
+		row[i] = row[i] " " $i
+}
+END {
+	for (i = 1; i <= n; i++)
+		print row[i]
+}' ${1+"$@"}

apps/pc_fit.py

@@ -87,7 +87,7 @@
 ################################################################################
 
 
-pcrv = pc_fit(x[indtrn], y[indtrn], order=order, pctype=pctype, method=method, eta=1.e-3)
+pcrv, lregs = pc_fit(x[indtrn], y[indtrn], order=order, pctype=pctype, method=method, eta=1.e-3)
 ypred=pcrv.function(x)
 
 if nout>=1000:
@@ -113,7 +113,7 @@
     for isam in range(0, nsam, nevery):
         f = plt.figure(figsize=(12,4))
         plt.plot(xc, y[isam,:], 'bo-', ms=8, label='Model')
-        plt.plot(xc, ypred[isam,:], 'go-', ms=8, label='PC apprx.')
+        plt.plot(xc, ypred[isam,:], 'go-', ms=8, label='PC Apprx.')
         plt.title(f'Sample #{isam+1}')
         plt.xlabel(xlabel)
         plt.ylabel(ylabel)
@@ -128,13 +128,6 @@
 
 plot_sens(mainsens,pars,cases,vis="bar", ncol=5, par_labels=pnames, case_labels=[str(i) for i in cases], lbl_size=25, xticklabel_size=18, legend_size=18, figname='sens_pc.png')
 
-# plot_sens(sensdata, pars, cases,
-#               vis="bar", reverse=False, topsens=None,
-#               par_labels=None, case_labels=None, colors=None,
-#               xlbl='', title='', grid_show=True,
-#               legend_show=2, legend_size=10, ncol=4,
-#               lbl_size=22, yoffset=0.1,
-#               xdatatick=None, xticklabel_size=None, xticklabel_rotation=0,
-#               figname='sens.png')
 
 savepk(pcrv, 'pcrv')
+savepk(lregs, 'lregs')

apps/pc_prep.py

@@ -1,34 +1,45 @@
 #!/usr/bin/env python
 
-################################################################################
-# Input PC generation given mvn, samples or marginal pc
-################################################################################
-
-
-
-# Input
-# Usage  : pc_prep.py <format> <filename> <input_pcorder>
-# e.g.   : pc_prep.py marg marg_pc.txt 3
-#        : pc_prep.py sam inp_sam.txt 3
-#        : pc_prep.py mvn mean.txt cov.txt
-# Output : pcf.txt
-
-
 import os
 import sys
+import argparse
 import numpy as np
 
 from pytuq.utils.mindex import get_mi
 from pytuq.utils.xutils import safe_cholesky
+from pytuq.workflows.fits import pc_ros
+
+################################################################################
+################################################################################
+################################################################################
+
+usage_str = 'Input PC generation given mvn, samples or marginal PC.'
+parser = argparse.ArgumentParser(description=usage_str)
+parser.add_argument("-f", "--fmt", dest="input_format", type=str, default='sam', help="Input format", choices=['marg', 'sam', 'mvn'])
+parser.add_argument("-i", "--inp", dest="filename", type=str, default='xsam.txt', help="Input filename: marginal coefficients (if format is marg), samples (if format is sam), mean (if format is mvn).")
+parser.add_argument("-c", "--cov", dest="cov_filename", type=str, default='cov.txt', help="Covariance filename (relevant if format is mvn).")
+parser.add_argument("-p", "--pco", dest="pcorder", type=int, default=1, help="PC order (relevant if format is marg or sam).")
+parser.add_argument("-t", "--pct", dest="pctype", type=str, default='HG', help="PC type (relevant if format is sam).")
+args = parser.parse_args()
+
+################################################################################
+################################################################################
+################################################################################
 
-input_format=sys.argv[1]
-filename=sys.argv[2]
+input_format=args.input_format
+filename=args.filename
 
 if input_format == "marg" or input_format == "sam":
-	input_pcorder=int(sys.argv[3])
-elif input_format == "mvn":
-	filename2=sys.argv[3]
+	pcorder=args.pcorder
+if input_format == "mvn":
+	cov_filename=args.cov_filename
+if input_format == "sam":
+	pctype = args.pctype
+
 
+################################################################################
+################################################################################
+################################################################################
 
 if input_format=="marg":
 
@@ -50,8 +61,8 @@
 		margpc_all.append(margpc_cur)
 
 
-	assert(input_pcorder >= maxord)
-	mindex_totalorder=get_mi(input_pcorder, dim)
+	assert(pcorder >= maxord)
+	mindex_totalorder=get_mi(pcorder, dim)
 
 	mindex=np.zeros((1,dim),dtype=int)
 	cfs=np.zeros((1,dim))
@@ -81,18 +92,15 @@
 
 
 elif input_format=="sam":
-	sam = np.loadtxt(filename)
-	ros = Rosenblatt(sam)
-	print("Not implemented yet")
-	sys.exit()
-	# TODO: make ex_ros_pc.py a function and use it here
+	sams = np.loadtxt(filename)
 
-	# cmd=uqtkbin+'pce_quad -o '+str(input_pcorder)+' -w HG -f '+filename+ ' > pcq.log; mv PCcoeff.dat pcf.txt'
-	# os.system(cmd)
+	nsam, dim = sams.shape
+	pcrv = pc_ros(sams, pctype=pctype, order=pcorder, nreg=nsam, bwfactor=1.0)
+	np.savetxt('pcf.txt', np.array(pcrv.coefs).T)
 
 elif input_format=="mvn":
 	mean = np.loadtxt(filename)
-	cov = np.loadtxt(filename2)
+	cov = np.loadtxt(cov_filename)
 
 	dim = mean.shape[0]
 
@@ -106,4 +114,4 @@
 	np.savetxt('pcf.txt', param_pcf)
 
 else:
-	print("pc_prep.py : Input format not recognized. Must be marg or sam.")
+	print("pc_prep.py : Input format not recognized. Must be marg or sam or mvn.")

apps/pc_sam.py

@@ -2,17 +2,26 @@
 
 import os
 import sys
+import argparse
 import numpy as np
 
 from pytuq.rv.pcrv import PCRV
 from pytuq.utils.mindex import get_mi
 
 
-pccf_file = sys.argv[1]
-pc_type = sys.argv[2]
-nsam = int(sys.argv[3])
+usage_str = 'Sampling multivariate PC given coefficient file.'
+parser = argparse.ArgumentParser(description=usage_str)
+parser.add_argument("-f", "--pcf", dest="pcf_file", type=str, default='pcf.txt', help="PC coefficient file: each column is PC coefficient vector for the corresponding dimension.")
+parser.add_argument("-t", "--pct", dest="pc_type", type=str, default='HG', help="PC type", choices=['LU', 'HG'])
+parser.add_argument("-n", "--nsam", dest="nsam", type=int, default=111, help="Number of requested samples.")
+args = parser.parse_args()
 
-pccf = np.loadtxt(pccf_file)
+
+pcf_file = args.pcf_file
+pc_type = args.pc_type
+nsam = args.nsam
+
+pccf = np.loadtxt(pcf_file)
 if len(pccf.shape)==1:
 	pccf = pccf[:, np.newaxis]
 

apps/plot_cov.py

@@ -7,7 +7,7 @@
 import scipy.stats as ss
 import matplotlib.pyplot as plt
 
-from pytuq.utils.plotting import myrc
+from pytuq.utils.plotting import myrc, plot_cov, plot_cov_tri
 
 myrc()
 
@@ -36,25 +36,14 @@
 dim_show = len(ind_show)
 
 
-f=3.
-
 for ii in range(dim_show):
     for jj in range(ii+1,dim_show):
         i, j = ind_show[ii], ind_show[jj]
-        x = np.linspace(mean[i]-f*np.sqrt(cov[i,i]), mean[i]+f*np.sqrt(cov[i,i]), 100)
-        y = np.linspace(mean[j]-f*np.sqrt(cov[j,j]), mean[j]+f*np.sqrt(cov[j,j]), 100)
-        X, Y = np.meshgrid(x, y)
-
-        try:
-            rv = ss.multivariate_normal([mean[i], mean[j]], [[cov[i,i], cov[i,j]],[cov[j,i], cov[j,j]]], allow_singular=True)
-            XY = np.dstack((X, Y))
-
-            Z = rv.pdf(XY)
-            plt.contour(X,Y,Z)
-            plt.xlabel('p'+str(i+1))
-            plt.ylabel('p'+str(j+1))
-            plt.savefig('cov_'+str(i)+'_'+str(j)+'.png')
-            plt.clf()
-
-        except ValueError:
-            print(f"Covariance for pair ({i},{j}) is not positive-semidefinite.")
+
+        mm = np.array([mean[i], mean[j]])
+        cc = np.array([[cov[i,i], cov[i,j]],[cov[j,i], cov[j,j]]])
+        plot_cov(mm, cc, f=3., pnames=[f'p{i}', f'p{j}'], ngr=100, savefig=True)
+        plt.clf()
+
+plot_cov_tri(mean[ind_show], cov[np.ix_(ind_show, ind_show)], names=[f'p{i}' for i in ind_show])
+

apps/plot_ens.py

@@ -5,28 +5,31 @@
 import numpy as np
 import matplotlib.pyplot as plt
 
+from pytuq.utils.plotting import myrc
 
+myrc()
 
 usage_str = 'Script to plot ensemble.'
 parser = argparse.ArgumentParser(description=usage_str)
-#parser.add_argument('ind_show', type=int, nargs='*',
-#                    help="indices of requested parameters (count from 0)")
-
 parser.add_argument("-y", "--ydata", dest="ydata", type=str, default='ytrain.dat',
                     help="Ydata file")
 
 args = parser.parse_args()
 
 ydata = np.loadtxt(args.ydata)
+if len(ydata.shape)==1:
+    ydata = ydata[:, np.newaxis]
 
 
 nsam, nout = ydata.shape
 
-ind_plot=np.arange(nout)
 
-nout_plot=len(ind_plot)
 
 for i in range(nsam):
-    plt.plot(np.arange(1, nout+1)[ind_plot], ydata[i, ind_plot], 'b-', lw=0.1)
+    plt.plot(np.arange(1, nout+1), ydata[i, :], 'bo-', markeredgecolor='w', lw=0.5)
 
+plt.xticks(np.arange(1, nout+1))
+plt.xtickslabels = [str(i) for i in range(1, nout+1)]
+plt.xlabel('Output Id')
+plt.ylabel('Output Value')
 plt.savefig('ensemble.png')

apps/plot_pcoord.py

@@ -72,7 +72,7 @@
 labels=labels[::args.every]
 
 for i in range(ndg):
-    print("Plotting %d / %d " % (i+1,ndg))
+    print(f"Plotting parallel coordinates for output {i+1} / {ndg}")
     #names=range(1+i*ndcut,min(1+(i+1)*ndcut,ndim+1))
 
     values=xdata[::args.every,i*ndcut:min((i+1)*ndcut, ndim)].T
@@ -90,7 +90,7 @@
         parallel_coordinates(pnames_this,
                              values[:, np.array(labels)==lab],
                              list(labels_),
-                             'pcoord_'+str(i+1)+'_lab'+lab+'.png')
+                             'pcoord_'+str(i+1)+'_lab_'+lab+'.png')
 
 
     #labels_only=labels[labels==True]