PBS job chaining

Often times we want to run a large job (or a few smaller jobs) on an HPC system but we cannot make one large submission due to constraints in the maximum number of cores we can use at a time. One way around this issue is to submit smaller jobs, some of which will be sent to the queue. However, this often leads to using more cores than we can or should at a time due to the queue management system fitting as many of our smaller jobs as possible in the available cores, which may prevent other users from using the system for a long time or lead to higher system usage than allowed.

To avoid this, one solution is to submit one smaller job at a time, waiting until the job(s) that is(are) currently running to be finished. This solution requires constant queue monitoring and is only as efficient as we can monitor the queue.

A better solution is the use of job chaining (qsub -W depend), which automatically submits a new job after the currently running one is finished. Below is an example of how to use job chaining for three jobs with one submission script:

#!/bin/bash
FIRST=$(qsub job1.sh)
echo $FIRST
SECOND=$(qsub -W depend=afterany:$FIRST job2.sh)
echo $SECOND
THIRD=$(qsub -W depend=afterany:$SECOND job3.sh)
echo $THIRD

In the example above, variables FIRST, SECOND and THIRD will receive the job ID, and line qsub -W depend=afterany:$FIRST job2.sh can be read as: submit job2.sh after job with ID $FIRST is finished with or without errors (afterany). It is also possible to set a job to be submitted only if the previous job ends with ok status or with an error.

Command	Rule
after	Execute current job after listed jobs have begun.
afterany	Execute current job after job has terminated with any status.
afterok	Execute current job after job has terminated without error.
afternotok	Execute current job after job has terminated with an error.

With this, you should be able to not lock everyone else out of the HPC system while not having to spend a week checking on the status of runs every few hours.

Root finding in MATLAB, R, Python and C++

In dynamical systems, we are often interested in finding stable points, or equilibria. Some systems have multiple equilibria. As an example, take the lake problem, which is modeled by the equation below where X_t is the lake P concentration, a_t are the anthropogenic P inputs, Y_t~LN(μ,σ²)  are random natural P inputs, b is the P loss rate, and q is a shape parameter controlling the rate of P recycling from the sediment. The first three terms on the right hand side make up the “Inputs” in the figure, while the last term represents the “Outputs.” A lake is in equilibrium when the inputs are equal to the outputs and the lake P concentration therefore is not changing over time.

For irreversible lakes this occurs at three locations, even in the absence of anthropogenic and natural inputs: an oligotrophic equilibrium, an unstable equilibrium (called the critical P threshold) and a eutrophic equilibrium (see figure below).

The unstable equilibrium in this case is called the critical P threshold because once it is crossed, it is impossible to return to an oligotrophic equilibrium by reducing anthropogenic and natural P inputs alone. In irreversible lakes like this, we would therefore like to keep the lake P concentration below the critical P threshold. How do we find the critical P threshold? With a root finding algorithm!

As stated earlier, the system above will be in equilibrium when the inputs are equal to the outputs and the P concentration is not changing over time, i.e. when

Therefore we simply need to find the zero, or “root” of the above equation.  Most of the methods for this require either an initial estimate or upper and lower bounds on the location of the root. These are important, since an irreversible lake will have three roots. If we are only interested in the critical P threshold, we have to make sure that we provide an estimate which leads to the unstable equilibrium, not either of the stable equilibria. If possible, you should plot the function whose root you are finding to make sure you are giving a good initial estimate or bounds, and check afterward to ensure the root that was found is the one you want! Here are several examples of root-finding methods in different programming languages.

In MATLAB, roots can be found with the function fzero(fun,x0) where ‘fun’ is the function whose root you want to find, and x0 is an initial estimate. This function uses Brent’s method, which combines several root-finding methods: bisection, secant, and inverse quadratic interpolation. Below is an example using the lake problem.

myfun = @(x,b,q) x^q/(1+x^q)-b*x;
b = 0.42;
q = 2.0;
fun = @(x) myfun(x,b,q);
pcrit = fzero(fun,0.75);

This returns pcrit = 0.5445, which is correct. If we had provided an initial estimate of 0.25 instead of 0.75, we would get pcrit = 2.6617E-19, basically 0, which is the oligotrophic equilibrium in the absence of anthropogenic and natural P inputs. If we had used 1.5 as an initial estimate, we would get pcrit = 1.8364, the eutrophic equilibrium.

In R, roots can be found with the function uniroot, which also uses Brent’s method. Dave uses this on line 10 of the function lake.eval in his OpenMORDM example. Instead of taking in an initial estimate of the root, this function takes in a lower and upper bound. This is safer, as you at least know that the root estimate will lie within these bounds. Providing an initial estimate that is close to the true value should do well, but is less predictable; the root finding algorithm may head in the opposite direction from what is desired.

b <- 0.42
q <- 2.0
pcrit <- uniroot(function(x) x^q/(1+x^q) - b*x, c(0.01, 1.5))$root

This returns pcrit = 0.5445145. Good, we got the same answer as we did with MATLAB! If we had used bounds of c(0.75, 2.0) we would have gotten 1.836426, the eutrophic equilibrium.

What if we had given bounds that included both of these equilibria, say c(0.5, 2.0)? In that case, R returns an error: ‘f() values at end points not of opposite sign’. That is, if the value returned by f(x) is greater than 0 for the lower bound, it must be less than 0 for the upper bound and vice versa. In this case both f(0.5) and f(2.0) are greater than 0, so the algorithm fails. What if we gave bounds for which one is greater than 0 and another less, but within which there are multiple roots, say c(-0.5,2.0)? Then R just reports the first one it finds, in this case pcrit = 0.836437, the eutrophic equilibrium. So it’s important to make sure you pick narrow enough bounds that include the root you want, but not roots you don’t!

In Python, you can use either scipy.optimize.root or scipy.optimize.brentq, which is what Jon uses on line 14 here. scipy.optimize.root can be used with several different algorithms, but the default is Powell’s hybrid method, also called Powell’s dogleg method. This function only requires an initial estimate of the root.

from scipy.optimize import root
b = 0.42
q = 2.0
pcrit = root(lambda x: x**(1+x**q) - b*x, 0.75)

scipy.optimize.root returns an object with several attributes. The attribute of interest to us is the root, represented by x, so we want pcrit.x. In this case, we get the correct value of 0.54454. You can play around with initial estimates to see how pcrit.x changes.

Not surprisingly, scipy.optimize.brentq uses Brent’s method and requires bounds as an input.

from scipy.optimize import brentq as root
b = 0.42
q = 2.0
pcrit = root(lambda x: x**(1+x**q) - b*x, 0.01, 1.5)

This just returns the root itself, pcrit = 0.5445. Again, you can play around with the bounds to see how this estimate changes.

In C++, Dave again shows how this can be done in the function ‘main-lake.cpp’ provided in the Supplementary Material to OpenMORDM linked from this page under the “Publications” section. On lines 165-168 he uses the bisect tool to find the root of the function given on lines 112-114. I’ve copied the relevant sections of his code into the function ‘find_Pcrit.cpp’ below.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <boost/math/tools/roots.hpp>

namespace tools = boost::math::tools;
using namespace std;

double b, q, pcrit;

double root_function(double x) {
  return pow(x,q)/(1+pow(x,q)) - b*x;
}

bool root_termination(double min, double max) {
  return abs(max - min) <= 0.000001;
}

int main(int argc, char* argv[])
{
  b = 0.42;
  q = 2.0;

  std::pair<double, double> result = tools::bisect(root_function, 0.01, 1.0, root_termination);
  pcrit = (result.first + result.second)/2;
  cout << pcrit << endl;
}

This yields the desired root of pcrit = 0.54454, but of course, changing the bounds may result in different estimates. In case you missed it, the take home message is to be careful about your initial estimate and bounds ;).

 

A quick example code to write data to a csv file in C++

Writing data output files in C++ is a little more complicated than doing so in languages like Python or Matlab. I’ve been learning some basic C++ this summer and I thought I’d share a simple example code for writing data to a CSV file.  See the code below, as well as some resources I’ve found helpful.

 

// write_csv.cpp
// updated 8/12/16

// the purpose of this code is to demonstrate how to write data
// to a csv file in C++

// inlcude iostream and string libraries
#include <iostream>
#include <fstream>
#include <string>

using namespace std;

// create an ofstream for the file output (see the link on streams for
// more info)
ofstream outputFile;
ofstream fs
// create a name for the file output
std::string filename = "exampleOutput.csv";

// create some variables for demonstration
int i;
int A;
int B;
int C;

int main(int argc, char** argv)
{
    // create and open the .csv file
    fs.open(outputFile,filename);
    
    // write the file headers
    outputFile << "Column A" << "," << "Column B" << "Column C" << std::endl;
    
    // i is just a variable to make numbers with in the file
    i = 1;
  
    // write data to the file
    for (int counter = 0; counter <  10; counter++)
    {
        // make some data
        A = i + 5;
        B = i + 10;
        C = i + 20;
        
        // write the data to the output file
        outputFile << A << "," << B << "," << C << std::endl;
        
        // increase i
        i = i * 5;
    }
    
    // close the output file
    outputFile.close();
    return 0;   
}

Helpful resources on writing data outputs in C++:

Basic info on input and output in C++: http://www.cplusplus.com/doc/tutorial/basic_io/

Short tutorial on writing txt file output: http://www.cplusplus.com/doc/tutorial/files/

Background on streams: http://www.cprogramming.com/tutorial/c++-iostreams.html

A nice tutorial from BU: https://www.cs.bu.edu/teaching/cs111/spring-2000/file-io/

 

 

 

 

 

 

Basic implementation of the parallel Borg MOEA

This post walks through the implementation of the different parallel versions of the Borg MOEA: master-slave and multi-master slave.  The implementation is demonstrated  using the DTLZ2 example provided with the multi-master source code. We will break this post in four sections, the fist one describes the Main file were the problem is defined and the required libraries are specified.  The second part describes the Makefile to compile and create the executables for the dtlz2 problem, the third section describes the Submission file to manage the distribution of jobs in a cluster, and finally we’ll cover a  Job submission example in the fourth section.

Both of the parallel Borg implementations are described in detail in the following papers:

Hadka, D., and Reed, P.M., “Large-scale Parallelization of the Borg MOEA for Many-Objective Optimization of Complex Environmental Systems”, Environmental Modelling & Software, v69, 353-369, 2015.

Reed, P.M. and Hadka, D., “Evolving Many-Objective Water Management to Exploit Exascale Computing”, Water Resources Research, v50, n10, 8367–8373, 2014.

1. Main file

1.1. Required headers

The main file refers to the file were we specify the main function called : dtlz2_mm.c,  dtlz2 is the test problem we are solving,  the mm refers to multi-master implementation.  First, you will need the mpi header (line 6), this is a message passing library required for parallel computers and clusters.  You will also need to provide the Borg multi-master header file: borgmm.h (line 7),  if your working  directory is different from that of the multi-master directory, you will need to specify the path, for intstance: “./mm-borg-moea/borgmm.h”.

#include <stdio.h>;
#include <stdlib.h>;
#include <math.h>;
#include <time.h>;
#include <math.h>;
#include <mpi.h>;
#include "borgmm.h";

1.2. Problem definition

In the following lines the DTLZ2 problem is defined as it would be in the serial version.  The following funciton is responsible for reading the decision variables and evaluating the problem.  We will be using the 2 objective version of this problem; however, you can scale it up.  The rule is nvars=nobjs+9.  Hence, if you want to try up to 5 objectives you can simply change lines 2 and 3 to be: nvars= 14 and nobjs=5.

#define PI 3.14159265358979323846
int nvars = 11;
int nobjs = 2;
void dtlz2(double* vars, double* objs, double* consts) {
	int i;
	int j;
	int k = nvars - nobjs + 1;
	double g = 0.0;

	for (i=nvars-k; i<nvars; i++) {
		g += pow(vars[i] - 0.5, 2.0);
	}

	for (i=0; i<nobjs; i++) {
		objs[i] = 1.0 + g;

		for (j=0; j<nobjs-i-1; j++) {
			objs[i] *= cos(0.5*PI*vars[j]);
		}

		if (i != 0) {
			objs[i] *= sin(0.5*PI*vars[nobjs-i-1]);
		}
	}
}

1.3. Enabling multiple seeds

The following lines enable to use random seeds as external arguments.  This links your main file to your submission file where multiple seeds are submitted to the cluster.  This segment is coupled with the submission file discussed in section 3.

int main(int argc, char* argv[]) {
unsigned int seed = atoi(argv[1]);

1.4. Variable declarations

In lines 1-2 of the following code, a couple of loop variables are declared, j and rank, which will be used later in section 1.7 and 1.9.  Also, the maximum number of function evaluations are declared. We will print out the runtime file and the file with the Pareto set; hence, the size of the output file names are also specified in lines 4-7.

int j;
int rank;
int NFE = 100000;
char outputFilename[256];
FILE* outputFile = NULL;
char runtime[256];
char timing[256];

1.5. Parallel borg parameters and core allocation

In the following segment, we specify some key parameters for the parallel Borg.  First, all multi-master runs need to call startup (in line 1).  The name of the variable argc stands for “argument count”; argc contains the number of arguments we pass to the program. The name of the variable argv stands for “argument vector”, which are the arguments that we pass to the program.  Then, the number of islands is specified in line 2.  If you want to use the master-slave configuration, you can simply assign one island.  The only difference with the pure master-slave code is that it uses uniform population samples, whereas the multi-master code uses latin hypercube population samples.  Line 3 specifies the maximum wallclock time estimated for your job.  Finally, global latin hypercube initialization is specified in line 5 to ensure each island gets a well sampled distribution of solutions.

BORG_Algorithm_ms_startup(&argc, &argv);
BORG_Algorithm_ms_islands(2);
BORG_Algorithm_ms_max_time(0.1);
BORG_Algorithm_ms_max_evaluations(NFE);
BORG_Algorithm_ms_initialization(INITIALIZATION_LATIN_GLOBAL);

Keep in mind the core allocation, if you have 32 available cores for a computation, when using a master-slave configuration, one core will be allocated to the master and the remaining  31 cores will be available for the function evaluations.  If you use a 2 multi-master configuration, one core will be allocated to the controller and two will be allocated to the masters (1 core per master), leaving 29 cores available for the function evaluations.  If you are using a small cluster, I wouldn’t recommend going higher than 4 masters.

1.6. Problem definition

The next segment creates the DTLZ2 problem.  In line 1, the number of decision variables, objectives and constraints are specified, the last argument, dtlz2, references the function that evaluates the DTLZ2 problem shown in section 1.2.  In lines 3-5 the lower and upper bounds for each decision variable are set to 0 and 1.  In lines 7-9 the epsilon values used by the Borg MOEA are specified to define the problem resolution for each objective to 0.01.

BORG_Problem problem = BORG_Problem_create(nvars, nobjs, 0, dtlz2);

for (j=0; j<nvars; j++) {
BORG_Problem_set_bounds(problem, j, 0.0, 1.0);
}

for (j=0; j<nobjs; j++) {
BORG_Problem_set_epsilon(problem, j, 0.01);
}

1.7. Printing runtime output

We specify the output frequency in line 1, this will print 100 generations, since we specified 100,000 maximum function evaluations, Borg will provide runtime output every 1000 NFE.  Lines 2 and 3 save the Pareto sets and the runtime dynamics to a file.  The %d gets replaced by the index of the seed and the %%d gets replaced by the index of the master, make sure that the sets and runtime folders exist in your working directory.

	BORG_Algorithm_output_frequency((int)NFE/100);
	sprintf(outputFilename, "./sets/DTLZ2_S%d.set", seed);
	sprintf(runtime, "./runtime/DTLZ2_S%d_M%%d.runtime", seed);
	BORG_Algorithm_output_runtime(runtime);

1.8. Parallelizing seeds

The MPI_Comm_rank routine gets the rank of this process. The rank is used to ensure each parallel process uses a different random seed.

    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    BORG_Random_seed(37*seed*(rank+1));

1.9. Printing the Pareto set

The optimization is performed by the multi-master Borg MOEA on the problem in line 1.  Then, only the controller process will return a non-NULL result.  The controller aggregates all of the Pareto optimal solutions generated by each master. Then print the Pareto optimal solutions to a separate file and frees any allocated memory.  Lines 13-15 shutdown the parallel processes and exit the program.

	BORG_Archive result = BORG_Algorithm_ms_run(problem);

	if (result != NULL) {
		outputFile = fopen(outputFilename, "w");
		if (!outputFile) {
			BORG_Debug("Unable to open final output file\n");
		}
		BORG_Archive_print(result, outputFile);
		BORG_Archive_destroy(result);
		fclose(outputFile);
	}

	BORG_Algorithm_ms_shutdown();
	BORG_Problem_destroy(problem);
	return EXIT_SUCCESS;
}

2. Makefile

The following makefile compiles the different versions of the Borg MOEA (serial, master-slave and multi-master slave) as well as the DTLZ2 examples for each version and generates the executables.  The CC compiler is required for the serial version while the mpicc compiler is required for the parallel versions.  From your terminal, access the mm-borg-moea directory and type make to compile the dtlz2 examples.

CC = gcc
MPICC = mpicc
CFLAGS = -O3
LDFLAGS = -Wl,-R,\.
LIBS = -lm
UNAME_S = $(shell uname -s)

ifneq (, $(findstring SunOS, $(UNAME_S)))
	LIBS += -lnsl -lsocket -lresolv
else ifneq (, $(findstring MINGW, $(UNAME_S)))
	# MinGW is not POSIX compliant
else
	POSIX = yes
endif

compile:
	$(CC) $(CFLAGS) $(LDFLAGS) -o dtlz2_serial.exe dtlz2_serial.c borg.c mt19937ar.c $(LIBS)

ifdef POSIX
	$(CC) $(CFLAGS) $(LDFLAGS) -o borg.exe frontend.c borg.c mt19937ar.c $(LIBS)
endif

	$(MPICC) $(CFLAGS) $(LDFLAGS) -o dtlz2_ms.exe dtlz2_ms.c borgms.c mt19937ar.c $(LIBS)
	$(MPICC) $(CFLAGS) $(LDFLAGS) -o dtlz2_mm.exe dtlz2_mm.c borgmm.c mt19937ar.c $(LIBS)
.PHONY: compile

 

3. Submission file

This is an example of a submission file.  Parallel Borg uses portable batch system (PBS) to manage the distribution of batch jobs across the available nodes in the cluster.  In the following script, the -N flag is the name of the job, the -l flag is the list of required resources, in this file we are requesting 1 node with 16 processors per node.  We are also requesting for 5 hours of wallclock time.  Index of the output and error files and the path of the output stream.  Line 6 specifies the path of the current working directory.  In lines 10 through 16 we run multiple seeds in a loop and call mpirun with  the name of the executable that was generated compiling the examples in section 2.  The following script is saved as a bash file, for instance: mpi-dtlz2.sh.

#!/bin/bash
#PBS -N dtlz2
#PBS -l nodes=1:ppn=16
#PBS -l walltime=5:00:00
#PBS -j oe
#PBS -o output

cd $PBS_O_WORKDIR

NSEEDS=9
SEEDS=$(seq 0 ${NSEEDS})

for SEED in ${SEEDS}
do
  mpirun ./dtlz2_mm.exe ${SEED}
done

4.  Submitting and managing jobs in a cluster

Once our submission file is ready, we can grant it permission using the following command:

chmod -x ./mpi-dtlz2.sh

Submit it to the cluster as such:

qsub mpi-dtlz2.sh

You can also delete a job using the qdel command and the job identifier:

qdel 24590

You can also hold a job:

qhold

Show status of the jobs:

qstat

Displays information about active, queued or recently completed jobs:

showq