Merge branch 'mpi'

Makefile

@@ -15,12 +15,15 @@ test: lib
 
 # LIBRARY
 
-lib: alg.o $(HEADERS) 
-	$(AR) $(AROPT) libaduprop.a alg.o
+lib: alg.o parallel.o $(HEADERS) 
+	$(AR) $(AROPT) libaduprop.a alg.o parallel.o
 
 alg.o: $(SRC_DIR)/alg.cpp $(HEADERS)
 	$(CXX) $(CFLAGS) -c $(SRC_DIR)/alg.cpp
 
+parallel.o: $(SRC_DIR)/parallel.cpp $(HEADERS)
+	$(CXX) $(CFLAGS) -c $(SRC_DIR)/parallel.cpp
+
 # DOCUMENTATION
 
 doc: $(HEADERS) doc/Doxyfile.in doc/mainpage.md
@@ -33,7 +36,7 @@ clean:
 
 cleanall:
 	- $(RM) *.o libaduprop.a
-	- cd test ; for i in `ls -d */` ; do cd $$i ; make clean ; cd .. ; done ; rm *.log
+	- cd test ; for i in `ls -d */` ; do cd $$i ; make clean ; cd .. ; done ; $(RM) *.log
 	- cd examples ; for i in `ls -d */` ; do cd $$i ; make clean ; cd .. ; done
 
 # BUILD EXAMPLES

Makefile.sample

@@ -1,10 +1,11 @@
-CXX = clang++
-CFLAGS = -g -O0 -Wall 
+CXX = mpicxx
+CFLAGS = -DAUPROP_MPI -g -O0 -Wall 
 # DEBUG
 # CFLAGS = -g -O0 -Wall -DDBUG
 MATH_LIBS = -llapack -lblas
 CODI_DIR = $(HOME)/git/CoDiPack
 ADUPROP_DIR = $(HOME)/git/aduprop
+EIGEN_DIR = $(ADUPROP_DIR)/eigen/install/include/eigen3
 
 # DON'T EDIT BELOW THIS LINE
 
@@ -14,7 +15,7 @@ endif
 
 INCLUDE_DIR = $(ADUPROP_DIR)/include
 SRC_DIR = $(ADUPROP_DIR)/src
-CFLAGS += -std=c++11 -I$(CODI_DIR)/include -I$(INCLUDE_DIR)
+CFLAGS += -std=c++11 -I$(CODI_DIR)/include -I$(EIGEN_DIR) -I$(INCLUDE_DIR) 
 CFLAGS += -I$(ADUPROP_DIR)/CONTRIB
 LDLIBS = -L$(ADUPROP_DIR) -laduprop $(MATH_LIBS)
 AR = ar

examples/plorenz/Makefile

@@ -0,0 +1,14 @@
+include ../../Makefile.inc
+
+CFLAGS+= -I$(PWD)
+
+all: user.hpp main.o
+	$(CXX) -o plorenz main.o $(LDLIBS)
+
+main.o: $(HEADERS) main.cpp
+	$(CXX) $(CFLAGS) -c main.cpp
+
+clean:
+	$(RM) *.o plorenz datas data_cov1 data_cov2 data_cov3 
+
+

examples/plorenz/lorenz96.py

@@ -0,0 +1,41 @@
+from scipy.integrate import odeint
+import matplotlib.pyplot as plt
+import numpy as np
+
+# these are our constants
+N = 4  # number of variables
+F = 2.0  # forcing
+
+def Lorenz96(x,t):
+
+  # compute state derivatives
+  d = np.zeros(N)
+  # first the 3 edge cases: i=1,2,N
+  d[0] = (x[1] - x[N-2]) * x[N-1] - x[0]
+  d[1] = (x[2] - x[N-1]) * x[0]- x[1]
+  d[N-1] = (x[0] - x[N-3]) * x[N-2] - x[N-1]
+  # then the general case
+  for i in range(2, N-1):
+      d[i] = (x[i+1] - x[i-2]) * x[i-1] - x[i]
+  # add the forcing term
+  d = d + F
+
+  # return the state derivatives
+  return d
+
+x0 = F*np.ones(N) # initial state (equilibrium)
+x0[0] += 0.05 # add small perturbation to 20th variable
+t = np.arange(0.0, 10.0, 0.01)
+
+x = odeint(Lorenz96, x0, t)
+print(x)
+
+# plot first three variables
+from mpl_toolkits.mplot3d import Axes3D
+fig = plt.figure()
+ax = fig.gca(projection='3d')
+ax.plot(x[:,0],x[:,1],x[:,2])
+ax.set_xlabel('$x_1$')
+ax.set_ylabel('$x_2$')
+ax.set_zlabel('$x_3$')
+plt.show()

examples/plorenz/main.cpp

@@ -0,0 +1,254 @@
+/*!
+   \file "power.cpp"
+   \brief "Jacobian, Hessian and Tensor accumulation using Automatic
+           Differentiation. This is the main function where everything comes 
+           together"
+   \author "Adrian Maldonado and Michel Schanen"
+   \date 21/11/2017
+*/
+
+#include <cmath>
+#include <codi.hpp>
+#include <iostream>
+#include <fstream>
+// AD library used
+// Rudimentary linear solver interface
+#include "ad.hpp"
+#include "cxxopts.hpp"
+#include "parallel.hpp"
+
+int main(int argc, char* argv[]) {
+  // Options parser
+
+
+  // activate timer in propagateAD
+
+  global_prof.activate("propagateAD");
+  global_prof.activate("propagateCOV");
+  global_prof.activate("propagateH");
+  global_prof.activate("propagateMU");
+  global_prof.activate("reduction");
+  cxxopts::Options options("UQ Power", "Perform UQ on power system with AD");
+
+  bool test_jac = false;
+  bool test_tensor1  = false;
+  bool test_tensor2  = false;
+  bool test_tensor3  = false;
+  bool tensor1  = false;
+  bool tensor2  = false;
+  bool tensor3  = false;
+
+  options.add_options()
+    ("test_jac", "Test inner jacobian", cxxopts::value<bool>(test_jac))
+    ("test_tensor1", "Compute AD and HC jacobian for 1 ts", cxxopts::value<bool>(test_tensor1))
+    ("test_tensor2", "Compute AD and HC hessian for 1 ts", cxxopts::value<bool>(test_tensor2))
+    ("test_tensor3", "Compute AD and HC tensor of order 3 for 1 ts",
+     cxxopts::value<bool>(test_tensor3))
+    ("tensor1", "Compute AD and HC jacobian", cxxopts::value<bool>(tensor1))
+    ("tensor2", "Compute AD and HC hessian", cxxopts::value<bool>(tensor2))
+    ("tensor3", "Compute AD and HC tensor of order 3",
+     cxxopts::value<bool>(tensor3));
+
+  auto result = options.parse(argc, argv);
+  // Variable declaration
+  // pVector<double> xold, x, F;
+  // pVector<double> x0;
+  // pMatrix<double> TMAT;
+  System sys;
+  if (argc < 3)
+  {
+    std::cout << "plorenz [options] [dimension N] [forcing F] [timesteps]" << std::endl;
+    std::cout << "Example: mpiexec -n 2 ./plorenz --tensor1 7 4.4 40000" << std::endl;
+    std::cout << "--tensorX: Use x-th (= 1,2 or 3) order derivatives" << std::endl;
+    return 1;
+  }
+  sys.dimension = atoi(argv[1]);
+  sys.F = atof(argv[2]);
+  size_t tsteps = atoi(argv[3]);
+  size_t dim = sys.dim();
+
+  pVector<double> xold(dim);
+  pVector<double> x(dim);
+  // F.alloc(sys.dimension);
+
+  // Read initial conditions from external file. Check for dimension
+  // consistency.
+
+  sys.ic(xold);
+  ad drivers(sys);
+  dim = sys.dimension;
+  size_t chunk = paduprop_getend(dim) - paduprop_getstart(dim);
+  size_t start = paduprop_getstart(dim);
+  size_t end = paduprop_getend(dim);
+   x = xold;
+  // jacobian test
+   if (test_jac)
+   {
+     drivers.jactest(xold);
+     return 0;
+  }
+  if (test_tensor1) {
+    pMatrix<double> J(dim, dim);
+    pMatrix<double> J_fd(dim, dim);
+    J.zeros();
+    drivers.t1s_driver(x, J);
+    std::cout << "Jacobian using 1st order AD" << endl;
+    std::cout << "-----------------" << std::endl;
+    std::cout << J << std::endl;
+
+    J_fd.zeros();
+    drivers.fdJ_driver(x, J_fd);
+    std::cout << "Jacobian using FD" << std::endl;
+    std::cout << "-----------------" << std::endl;
+    std::cout << "J" << std::endl << J << std::endl;
+    double diff = (J - J_fd).norm();
+    cout << "norm(J - J_fd)" << endl << diff << endl;
+    if(diff > 1e-4) {
+      cerr << "ERROR: HC and AD Jacobian differ." << endl;
+      return 1;
+    }
+    return 0;
+  }
+  if (test_tensor2) {
+    pMatrix<double> J(dim, dim);
+    pTensor3<double> H(dim, dim, chunk);
+    pTensor3<double> H_fd(dim, dim, chunk);
+    J.zeros();
+    H.zeros();
+    drivers.t2s_t1s_driver(x, J, H, start, end);
+    std::cout << "Hessian and Jacobian using 2nd order AD" << std::endl;
+    std::cout << "-----------------" << std::endl;
+    // std::cout << "H" << std::endl << H << std::endl;
+
+    H_fd.zeros();
+    drivers.fdH_driver(x, H_fd, start, end);
+
+    double diff = (H - H_fd).norm();
+    diff = diff*diff;
+    paduprop_sum(diff);
+    cout << "norm(H - H_fd)" << endl << sqrt(diff) << endl;
+    if(diff > 1e-4) {
+      cerr << "ERROR: HC and AD Hessian differ." << endl;
+      return 1;
+    }
+    return 0;
+  }
+  if (test_tensor3) {
+    pMatrix<double>  J(dim, dim);
+    pTensor3<double> H(dim, dim, dim);
+    pTensor4<double> T(dim, dim, dim, chunk);
+    pTensor4<double> T_fd(dim, dim, dim, chunk);
+    J.zeros();
+    H.zeros();
+    T.zeros();
+    std::cout << "Tensor, Hessian and Jacobian using 3rd order AD" << std::endl;
+    std::cout << "-----------------" << std::endl;
+    drivers.t3s_t2s_t1s_driver(x, J, H, T, start, end);
+    // std::cout << "T" << std::endl << T << std::endl;
+    // std::cout << std::endl;
+
+    T_fd.zeros();
+    drivers.fdT_driver(x, T_fd, start, end);
+    double diff = (T - T_fd).norm();
+    diff = diff*diff;
+    paduprop_sum(diff);
+    cout << "norm(T - T_fd)" << endl << sqrt(diff) << endl;
+    if(diff > 1e-5) {
+      cerr << "ERROR: HC and AD tensor differ." << endl;
+      return 1;
+    }
+    return 0;
+  }
+  // Co-variance
+  pMatrix<double> cv0(sys.dimension, sys.dimension);
+  cv0.zeros();
+  // Strings
+  std::string fcov;
+
+  // Where do we put this???
+  // We should only alocate if we're using all of these.
+  pTensor4<double> T;
+  pTensor3<double> H;
+  pMatrix<double>  J(dim, dim);
+  if (tensor2)
+  {
+    H=pTensor3<double>(dim, dim, chunk);
+  }
+  if (tensor3)
+  {
+    T=pTensor4<double>(dim, dim, dim, chunk);
+    H=pTensor3<double>(dim, dim, dim);
+  }
+
+  for (size_t i = 0; i < sys.dimension; ++i) 
+    cv0[i][i] = 0.0000001;
+
+  cv0[sys.dim()][sys.dimension] = 0.000001;
+
+  std::ofstream xfile;
+  std::ofstream covfile;
+  if (paduprop_getrank() == 0)
+  {
+    xfile.open("datas");
+    if (tensor2)
+    {
+      covfile.open("data_cov2");
+    }
+    else
+    {
+      if (tensor3)
+      {
+        covfile.open("data_cov3");
+      }
+      else
+      {
+        covfile.open("data_cov1");
+      }
+    }
+  }
+  for (size_t i = 0; i < tsteps; ++i) {
+    // Save mean in trajectory matrix
+    // for (size_t j = 0; j < sys.dimension; ++j) {
+    //   TMAT.set(j, i, x[j]);
+    // }
+
+    // Propagate
+    if (paduprop_getrank() == 0)
+    {
+      std::cout << "Step: " << i << "." << std::endl;
+    }
+    int degree = 1;
+    if(tensor1) degree = 1;
+    if(tensor2) degree = 2;
+    if(tensor3) degree = 3;
+    propagateAD(x, cv0, sys, J, H, T, drivers, degree);
+    if (paduprop_getrank() == 0)
+    {
+      // xfile << "X: " << i * sys.h << " " << x << std::endl;
+      // covfile << "COV: " << i * sys.h << " ";
+      for (size_t j = 0; j < cv0.nrows(); j++)
+      {
+        xfile << x[j] << " ";
+        covfile << cv0[j][j] << " ";
+      }
+      xfile << std::endl;
+      covfile << std::endl;
+    }
+  }
+  // cout << cv0;
+  if(paduprop_getrank() == 0) {  
+    xfile.close();
+    covfile.close();
+    double compare, diff;
+    std::ifstream infile;
+    infile.open("solution.txt");
+    infile >> compare;
+    diff = compare - cv0.norm();
+    std::cout << std::setprecision(16) << "Sol: " << compare << ". " << std::endl;
+    std::cout << std::setprecision(16) << "cv0 norm: " << cv0.norm() << "." << std::endl;
+    std::cout << std::setprecision(16) << "Diff: " << diff << "." << std::endl;
+    infile.close();
+  }
+  std::cout << global_prof << std::endl;
+}
+

examples/plorenz/plot.cfg

@@ -0,0 +1,7 @@
+set term png
+set output output_file
+
+set pm3d map
+set xlabel "Variable"
+set ylabel "Time"
+plot filename matrix with image 

examples/plorenz/plot.sh

@@ -0,0 +1,7 @@
+#!/bin/bash
+
+gnuplot -e "filename='datas' ; output_file='datas.png'" -p -c plot.cfg
+gnuplot -e "filename='data_cov1' ; output_file='data_cov1.png'" -p -c plot.cfg
+gnuplot -e "filename='data_cov2' ; output_file='data_cov2.png'" -p -c plot.cfg
+gnuplot -e "filename='data_cov3' ; output_file='data_cov3.png'" -p -c plot.cfg
+

examples/plorenz/run.sh

@@ -0,0 +1,15 @@
+#!/bin/bash
+
+# Run Lorentz 96 and compute covariance with 1st and 2nd order derivatives
+
+# ARG1: dimension
+# ARG2: forcing
+# ARG3: timesteps
+
+./plorenz --tensor1 $1 $2 $3 | grep 'COV: ' | sed 's/[COV:|\[|,]//g' | sed 's/\]//' | awk '{$1=""; print $0}' > data_cov1
+
+./plorenz --tensor2 $1 $2 $3 | grep 'COV: ' | sed 's/[COV:|\[|,]//g' | sed 's/\]//' | awk '{$1=""; print $0}' > data_cov2
+
+./plorenz --tensor3 $1 $2 $3 | grep 'COV: ' | sed 's/[COV:|\[|,]//g' | sed 's/\]//' | awk '{$1=""; print $0}' > data_cov3
+
+./plorenz $1 $2 $3 | grep 'X: ' | sed 's/[X:|\[|,]//g' | sed 's/\]//' | awk '{$1=""; print $0}' > datas