# Profiling C++ code with Callgrind

Often times, we have to write code to perform tasks whose complexity vary from mundane, such as retrieving and organizing data, to highly complex, such as simulations CFD simulations comprising the spine of a project. In either case, depending on the complexity of the task and amount of data to be processed, it may happen for the newborn code to leave us staring at an underscore marker blinking gracefully for hours on a command prompt during its execution until the results are ready, leading to project schedule delays and shortages of patience. Two standard and preferred approaches to the problem of time intensive codes are to simplify the algorithm and to make the code more efficient. In order to better select the parts of the code to work on, it is often useful to first find the parts of the code in which more time by profiling the code.

In this post, I will show how to use Callgrind, part of Valgrind, and KCachegrind to profile C/C++ codes on Linux — unfortunately, Valgrind is not available for Windows or Mac, although it can be ran on cluster from which results can be downloaded and visualized on Windows with QCachegrind. The first step is to install Valgrind and KCachegrind by typing the following commands in the terminal of a Debian based distribution, such as Ubuntu (equivalent yum commands area available for Red Hat based distributions):

$sudo apt-get install valgrind$ sudo apt-get install kcachegrind


Now that the required tools are installed, the next step is to compile your code with GCC/G++ (with a make file, cmake, IDE or by running the compiler directly from the terminal) and then run the following command in a terminal (type ctrl+shift+T to open the terminal):

$valgrind --tool=callgrind path/to/your/compiled/program program_arguments  Callgrind will then run your program with some instrumentation added to its execution to measure time expenditures and cache use by each function in your code. Because of the instrumentation, Your code will take considerably longer to run under Callgrind than it typically would, so be sure to run a representative task that is as small as possible when profiling your code. During its execution, Callgrind will output a report similar to the one below on terminal itself: ==12345== Callgrind, a call-graph generating cache profiler ==12345== Copyright (C) 2002-2015, and GNU GPL'd, by Josef Weidendorfer et al. ==12345== Using Valgrind-3.11.0 and LibVEX; rerun with -h for copyright info ==12345== Command: path/to/your/compiled/program program_arguments ==12345== ==12345== For interactive control, run 'callgrind_control -h'. IF YOUR CODE OUTPUTS TO THE TERMINAL, THE OUTPUT WILL BE SHOWN HERE. ==12345== ==12345== Events : Ir ==12345== Collected : 4171789731 ==12345== ==12345== I refs: 4,171,789,731  The report above shows that it collected 4 billion events in order to generate the comprehensive report saved in the file callgrind.out.12345 — 12345 is here your process id, shown in the report above. Instead of submerging your soul into a sea of despair by trying to read the output file in a text editor, you should load the file into KCachegrind by typing: $ kcachegrind calgrind.out.12345


You should now see a screen like the one below:

The screenshot above shows the profiling results for my code. The left panel shows the functions called by my code sorted by total time spent inside each function. Because functions call each other, callgrind shows two cost metrics as proxies for time spent in each function: Incl., showing the total cost of a function, and self, showing the time spent in each function itself discounting the callees. By clicking on “Self” to order to functions by the cost of the function itself, we sort the functions by the costs of their own codes, as shown below:

Callgrind includes functions that are native to C/C++ in its analysis. If one of them appears in the highest positions of the left panel, it may be the case to try to use a different function or data structure that performs a similar task in a more efficient way. Most of the time, however, our functions are the ones in most of the top positions in the list. In the example above, we can see that a possible first step I can take to improve the time performance of my code is to make function “ContinuityModelROF::shiftStorage” more efficient. A few weeks ago, however, the function “ContinuityModel::continuityStep” was ranked first with over 30% of the cost, followed by a C++ map related function. I then replaced a map inside that function by a pointer vector, resulting in the drop of my function’s cost to less than 5% of the total cost of the code.

In case KCachegrind shows that a given function that is called from multiple places in the code is costly, you may want to know which function is the main culprit behind the costly calls. To do this, click on the function of interest (in this case, “_memcpy_sse2_unalight”) in the left panel, and then click on “Callers” in the right upper panel and on “Call Graph” in the lower right panel. This will show in list and graph forms the calls made to the function by other functions, and the asociated percent costs. Unfortunately, I have only the function “ContinuityModelROF::calculateROF” calling “_memcpy_sse2_unalight,” hence the simple graph, but the graph would be more complex if multiple functions made calls to “_memcpy_sse2_unalight.”

I hope this saves you at least the time spend reading this post!

# 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 Xt is the lake P concentration, at are the anthropogenic P inputs, Yt~LN(μ,σ2)  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

$X_{t+1} - X_t = \frac{X^q_t}{1+X^q_t} - bX_t = 0$

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 ;). # Getting started with C and C++ I’ve been learning C and C++ recently and I thought I’d share my experience learning these languages through a post. Prior to learning C and C++, I had experience in Python and Matlab, but this was my first foray into lower level languages. In my attempts to learn each language I made my way through several courses and books available online; some I found very helpful, others not so much. In this post I’ll detail my experiences with each resource and provide some tips that may help those planning on learning these languages. ### Main takeaways To learn the languages, I used four main resources. Online courses from Lynda.com, a book titled Learn C the Hard Way, a book commonly known as K&R 2 and tutorials from cplusplus.com. For those who do not have the time or desire to read this whole post, I found the following resources to be the most useful: #### For C: Learning C the Hard Way K&R 2 #### For C++: Foundations of Programming: Object – Oriented Design (you need a lynda.com login to access) Up and Running with C++ (again, you need a lynda.com login to access) cplusplus.com Everyone’s learning style is different, but I found that I learned the languages much faster by coding examples myself, rather than watching someone walk through a script. I also found that courses that taught straight from the command line were more effective than courses that taught through an IDE. When using an IDE, I often found myself spending more time working through glitches or nuances within the IDE than learning the languages themselves. I’ll detail my experiences with each resource below. I’ll start with C resources, then discuss C++. ### Resources for Learning C #### C Essential Training – From Lynda.com I started my training by taking course on Lynda.com titled “C Essential Training”. Lynda.com is an online educational website with thousands of videos, many of which focus on programming. The service is free to Cornell students and graduate students (though I checked and unfortunately neither PSU, UC Davis nor CU Boulder have agreements with the site). I found the course to be well structured and I felt that the instructor presented the material clearly and understandably. Despite this, I do not feel that the course did an effective job teaching me C. The main problem I had with the course was its reliance on the Eclipse IDE. Eclipse seems like a fine IDE, but I don’t plan to use an IDE when writing code and I didn’t want to spend the time taking a separate course to learn its intricacies (though Lynda.com does have a full course devoted to Eclipse). Throughout the course, I kept finding myself having small Eclipse problems (e.g. not being able to change the project I was working on or having compiler errors) that were not hard to solve, but were never discussed in the lectures. I was able to solve each problem by doing some research online, but each little problem took me time to resolve and was mentally taxing. After spending 30 minutes looking up an Eclipse fix, I was not in the mood to go troubleshooting interesting C questions . Another problem with using Eclipse is that the user is never forced to write their own makefiles, an omission that seems like it could really hurt someone who plans to run C programs through the command line. In summary, I would not recommend taking this course unless you are either thoroughly versed in Eclipse or plan to write all of your code through Eclipse. #### Learning C the Hard Way The next resource I used to learn C was a book that Jazmin pointed me to called Learning C the Hard Way by Zed A. Shaw (after some poking around I found this had been mentioned previously on this blog). The book is laid out as a tutorial, where each chapter teaches a new C concept (it’s really a C course in book form).The author takes a slightly nontraditional teaching approach in that he makes you write the code first, then explains in detail what you just wrote. I found this very hands on teaching method extremely helpful. When I wrote the code myself, I was forced to focus on every detail of the code (something that is very important in a language like C). I also was able to learn which concepts were genuinely challenging for me and which concepts I needed more work on. When I watched the Lynda.com lectures, I’d often believe I understood a concept, only to find out later that I had misinterpreted the instructors lesson. The book does not use an IDE, but rather writes code in a text editor (I used Sublime Text) and runs them on the Unix command line. The author provides a succinct introduction to makefiles and how to use them, which was refreshing after the Eclipse based course that never mention makefiles or compilers. Overall I found the teaching method employed by the book to be very effective, and I would highly recommend it to new C users. I should note however, that there seems to be some controversy surrounding the book. If you google “Learning C the hard way” you’ll find some very heated exchanges between the author and a blogger who criticized the book’s teaching methodology. The blogger had two main criticisms of the book; first that it over simplified and inaccurately presented key C concepts, and second, that the author failed to teach accepted programming standards for the C language. Mr. Shaw’s rebuttal was that the book’s purpose was to teach people get people comfortable with C and begin actually coding with it, then once they are literate, have them go back and learn more about the language’s nuances. I personally agree with Mr. Shaw on this point, though I don’t have a background in computer science so my opinion is only that of an beginner. Many of the criticisms of the book seemed to come from the perspective of an advanced coder who is unable to see the language through the eyes of a beginner. Mr. Shaw’s explanations might be over simplified, but they do a good job demystifying many of the most foreign aspects of C. I think that use of this book should be supplemented with other sources, especially documents on accepted C coding standards, but if you’re looking for a quick way to get on your feet with C and gain some confidence, then the book is a great resource. I used a free beta version of the book which can be found here: http://c.learncodethehardway.org/book/ but you can also purchase the book from the author here: https://www.amazon.com/Learn-Hard-Way-Practical-Computational/dp/0321884922 I found the beta version to be just fine, but there were some minor errors and some sections were clearly under construction. The blog post criticizing the book can be found here: http://hentenaar.com/dont-learn-c-the-wrong-way #### K&R 2 A resource that I discovered through reading the exchanges between the Shaw and his critics was “The C Programming Language” by Brian W. Kernighan and Dennis M. Ritchie (commonly referred to as K&R 2 which is what I’ll call it for the rest of the post). One of the Authors of this book, Dennis Ritchie, actually coauthored the C language and this book is talked of as the go to authority of all matters C. Mr. Shaw devoted a whole chapter of “Learning C the Hard way” to bashing this book, but I found its layout and explanations quite accessible and useful. I did not find the tutorials as direct as “Learning C the Hard Way”, but I found it to be a helpful supplement. ### Resources for Learning C++ #### Foundations of Programming: Object-Oriented Design – From Lynda.com A main difference between C and C++ is that C++ is an object oriented language. I had some very basic experience in using object oriented programming, but was looking for a refresher before learning C++. “Foundations of Programming: Object-Oriented Design” was an excellent course that taught me all I wanted to know about object-oriented programming and more. The course is purely conceptual and does not teach any actual code or contain any practice problems. It presents the concepts in a simple yet comprehensive manner that I found very helpful. I would highly recommend this course to anyone hoping to learn or brush up their knowledge of how object-oriented programming works. #### Up and Running with C++ – From Lynda.com This course was very similar in layout to the C course from Lynda.com, and I have the same criticisms. The entire course used Eclipse, and I kept having minor problems that were never addressed by the lectures but prevented me from executing my code. I did feel like I was able to learn the basic tools I needed from the lectures, but I would have gotten much more out of the class if it had been taught through the command line. I also felt that the course was sparse on exercises and heavy on lectures. I found that I got much less out of watching the instructor write code than being forced to write out the code myself (as Learning C the Hard Way forces you to do). #### cplusplus.com This resource is linked often in older posts on this blog, and I found it helpful in answering C++ questions I had after finishing the Lynda.com courses. I did not find that tutorial had the most helpful narration of how one may learn C++ from scratch, but it has very succinct definitions of many C++ components and was helpful as a reference. I think this site is the one I will look to most when troubleshooting future C++ code. ### Final thoughts I’d like to note that I found WRASEMAN’s post on makefiles a few weeks back to be quite helpful. From my limited experience, ensuring that your code compiles correctly can be one of the most challenging parts of using a lower level language and the post has some excellent resources that explain makefiles are and how they can be used. I know there are a lot of contributors and readers of this blog who are much more versed in C and C++ than I am, so if you’d like to chime in on useful resources, please do so in the comments. # Using HDF5/zlib Compression in NetCDF4 Not too long ago, I posted an entry on writing NetCDF files in C and loading them in R. In that post, I mentioned that the latest and greatest version of NetCDF includes HDF5/zlib compression, but I didn’t say much more beyond that. In this post, I’ll explain briefly how to use this compression feature in your NetCDF4 files. Disclaimer: I’m not an expert in any sense on the details of compression algorithms. For more details on how HDF5/zlib compression is integrated into NetCDF, check out the NetCDF Documentation. Also, I’ll be assuming that the NetCDF4 library was compiled on your machine to enable HDF5/zlib compression. Details on building and installing NetCDF from source code can be found in the documentation too. I will be using code similar to what was in my previous post. The code generates three variables (x, y, z) each with 3 dimensions. I’ve increased the size of the dimensions by an order of magnitude to better accentuate the compression capabilities.  // Loop control variables int i, j, k; // Define the dimension sizes for // the example data. int dim1_size = 100; int dim2_size = 50; int dim3_size = 200; // Define the number of dimensions int ndims = 3; // Allocate the 3D vectors of example data float x[dim1_size][dim2_size][dim3_size]; float y[dim1_size][dim2_size][dim3_size]; float z[dim1_size][dim2_size][dim3_size]; // Generate some example data for(i = 0; i < dim1_size; i++) { for(j = 0; j < dim2_size; j++) { for(k = 0; k < dim3_size; k++) { x[i][j][k] = (i+j+k) * 0.2; y[i][j][k] = (i+j+k) * 1.7; z[i][j][k] = (i+j+k) * 2.4; } } } Next is to setup the various IDs, create the NetCDF file, and apply the dimensions to the NetCDF file. This has not changed since the last post.  // Allocate space for netCDF dimension ids int dim1id, dim2id, dim3id; // Allocate space for the netcdf file id int ncid; // Allocate space for the data variable ids int xid, yid, zid; // Setup the netcdf file int retval; if((retval = nc_create(ncfile, NC_NETCDF4, &ncid))) { ncError(retval); } // Define the dimensions in the netcdf file if((retval = nc_def_dim(ncid, "dim1_size", dim1_size, &dim1id))) { ncError(retval); } if((retval = nc_def_dim(ncid, "dim2_size", dim2_size, &dim2id))) { ncError(retval); } if((retval = nc_def_dim(ncid, "dim3_size", dim3_size, &dim3id))) { ncError(retval); } // Gather the dimids into an array for defining variables in the netcdf file int dimids[ndims]; dimids[0] = dim1id; dimids[1] = dim2id; dimids[2] = dim3id; Here’s where the magic happens. The next step is to define the variables in the NetCDF file. The variables must be defined in the file before you tag it for compression.  // Define the netcdf variables if((retval = nc_def_var(ncid, "x", NC_FLOAT, ndims, dimids, &xid))) { ncError(retval); } if((retval = nc_def_var(ncid, "y", NC_FLOAT, ndims, dimids, &yid))) { ncError(retval); } if((retval = nc_def_var(ncid, "z", NC_FLOAT, ndims, dimids, &zid))) { ncError(retval); } Now that we’ve defined the variables in the NetCDF file, let’s tag them for compression.  // OPTIONAL: Compress the variables int shuffle = 1; int deflate = 1; int deflate_level = 4; if((retval = nc_def_var_deflate(ncid, xid, shuffle, deflate, deflate_level))) { ncError(retval); } if((retval = nc_def_var_deflate(ncid, yid, shuffle, deflate, deflate_level))) { ncError(retval); } if((retval = nc_def_var_deflate(ncid, zid, shuffle, deflate, deflate_level))) { ncError(retval); } The function nc_def_var_deflate() performs this. It takes the following parameters: • int ncid – The NetCDF file ID returned from the nc_create() function • int varid – The variable ID associated with the variable you would like to compress. This is returned from the nc_def_var() function • int shuffle – Enables the shuffle filter before compression. Any non-zero integer enables the filter. Zero disables the filter. The shuffle filter rearranges the byte order in the data stream to enable more efficient compression. See this performance evaluation from the HDF group on integrating a shuffle filter into the HDF5 algorithm. • int deflate – Enable compression at the compression level indicated in the deflate_level parameter. Any non-zero integer enables compression. • int deflate_level – The level to which the data should be compressed. Levels are integers in the range [0-9]. Zero results in no compression whereas nine results in maximum compression. The rest of the code doesn’t change from the previous post.  // OPTIONAL: Give these variables units if((retval = nc_put_att_text(ncid, xid, "units", 2, "cm"))) { ncError(retval); } if((retval = nc_put_att_text(ncid, yid, "units", 4, "degC"))) { ncError(retval); } if((retval = nc_put_att_text(ncid, zid, "units", 1, "s"))) { ncError(retval); } // End "Metadata" mode if((retval = nc_enddef(ncid))) { ncError(retval); } // Write the data to the file if((retval = nc_put_var(ncid, xid, &x[0][0][0]))) { ncError(retval); } if((retval = nc_put_var(ncid, yid, &y[0][0][0]))) { ncError(retval); } if((retval = nc_put_var(ncid, zid, &z[0][0][0]))) { ncError(retval); } // Close the netcdf file if((retval = nc_close(ncid))) { ncError(retval); } So the question now is whether or not it’s worth compressing your data. I performed a simple experiment with the code presented here and the resulting NetCDF files: 1. Generate the example NetCDF file from the code above using each of the available compression levels. 2. Time how long the code takes to generate the file. 3. Note the final file size of the NetCDF. 4. Time how long it takes to load and extract data from the compressed NetCDF file. Below is a figure illustrating the results of the experiment (points 1-3). Before I say anything about these results, note that individual results may vary. I used a highly stylized data set to produce the NetCDF file which likely benefits greatly from the shuffle filtering and compression. These results show a compression of 97% – 99% of the original file size. While the run time did increase, it barely made a difference until hitting the highest compression levels (8,9). As for point 4, there was only a small difference in load/read times (0.2 seconds) between the uncompressed and any of the compressed files (using ncdump and the ncdf4 package in R). There’s no noticeable difference among the load/read times for any of the compressed NetCDF files. Again, this could be a result of the highly stylized data set used as an example in this post. For something more practical, I can only offer anecdotal evidence about the compression performance. I recently included compression in my current project due to the large possible number of multiobjective solutions and states-of-the-world (SOW). The uncompressed file my code produced was on the order of 17.5 GB (for 300 time steps, 1000 SOW, and about 3000 solutions). I enabled compression of all variables (11 variables – 5 with three dimensions and 6 with two dimensions – compression level 4). The next run produced just over 7000 solutions, but the compressed file size was 9.3 GB. The down side is that it took nearly 45 minutes to produce the compressed file, as opposed to 10 minutes with the previous run. There are many things that can factor into these differences that I did not control for, but the results are promising…if you’ve got the computer time. I hope you found this post useful in some fashion. I’ve been told that compression performance can be increased if you also “chunk” your data properly. I’m not too familiar with chunking data for writing in NetCDF files…perhaps someone more clever than I can write about this? Acknowledgement: I would like to acknowledge Jared Oyler for his insight and helpful advice on some of the more intricate aspects of the NetCDF library. # From Writing NetCDF Files in C to Loading NetCDF Files in R # So much data from such little models… It’s been my experience that even simple models can generate lots of data. If you’re a regular reader of this blog, I can imagine you’ve had similar experiences as well. My most recent experience with this is the work I’ve done with the Dynamic Integrated Climate-Economic model (DICE). I had inherited a port of the 2007 version of the model, which would print relevant output to the screen. During my initial runs with the model, I would simply redirect the output to ascii files for post-processing. I knew that eventually I would be adding all sorts of complexity to this model, ultimately leading to high-dimensional model output and rendering the use of ascii files as impractical. I knew that I would need a better way to handle all this data. So in updating the model to the 2013 version, I decided to incorporate support for netCDF file generation. You can find details about the netCDF file format through Unidata (a University Cooperation for Atmospheric Research [UCAR] Community Program) and through some of our previous blog posts (here, here, and here). What’s important to note here is that netCDF is a self-describing file format designed to manage high-dimensional hierarchical data sets. I had become accustomed to netCDF files in my previous life as a meteorologist. Output from complex numerical weather prediction models would often come in netCDF format. While I had never needed to generate my own netCDF output files, I found it incredibly easy and convenient to process them in R (my preferred post-processing and visualization software). Trying to incorporate netCDF output support in my simple model seemed daunting at first, but after a few examples I found online and a little persistence, I had netCDF support incorporated into the DICE model. The goal of this post is to guide you through the steps to generate and process a netCDF file. Some of our earlier posts go through a similar process using the Python and Matlab interfaces to the netCDF library. While I use R for post-processing, I generally use C/C++ for the modeling; thus I’ll step through generating a netCDF file in C and processing the generated netCDF file in R on a Linux machine. Edit: I originally put a link to following code at the bottom of this post. For convenience, here’s a link to the bitbucket repository that contains the code examples below. # Writing a netCDF file in C… ## Confirm netCDF installation First, be sure that netCDF is installed on your computing platform. Most scientific computing clusters will have the netCDF library already installed. If not, contact your system administrator to install the library as a module. If you would like to install it yourself, Unidata provides the source code and great documentation to step you through the process. The example I provide here isn’t all that complex, so any recent version (4.0+) should be able to handle this with no problem. ## Setup and allocation ### Include the header files With the netCDF libraries installed, you can now begin to code netCDF support into your model. Again, I’ll be using C for this example. Begin by including the netCDF header file with your other include statements: #include <stdlib.h> #include <stdio.h> #include <netcdf.h>  ### Setup an error handler The netCDF library includes a nice way of handling possible errors from the various netCDF functions. I recommend writing a simple wrapper function that can take the returned values of the netCDF functions and produce the appropriate error message if necessary: void ncError(int val) { printf("Error: %s\n", nc_strerror(val)); exit(2); } ### Generate some example data Normally, your model will have generated important data at this point. For the sake of the example, let’s generate some data to put into a netCDF file:  // Loop control variables int i, j, k; // Define the dimension sizes for // the example data. int dim1_size = 10; int dim2_size = 5; int dim3_size = 20; // Define the number of dimensions int ndims = 3; // Allocate the 3D vectors of example data float x[dim1_size][dim2_size][dim3_size]; float y[dim1_size][dim2_size][dim3_size]; float z[dim1_size][dim2_size][dim3_size]; // Generate some example data for(i = 0; i < dim1_size; i++) { for(j = 0; j < dim2_size; j++) { for(k = 0; k < dim3_size; k++) { x[i][j][k] = (i+j+k) * 0.2; y[i][j][k] = (i+j+k) * 1.7; z[i][j][k] = (i+j+k) * 2.4; } } } This generates three variables, each with three different size dimensions. Think of this, for example, as variables on a 3-D map with dimensions of [latitude, longitude, height]. In my modeling application, my dimensions were [uncertain state-of-the-world, BORG archive solution, time]. ### Allocate variables for IDs Everything needed in creating a netCDF file depends on integer IDs, so the next step is to allocate variables for the netCDF file id, the dimension ids, and the variable ids: // Allocate space for netCDF dimension ids int dim1id, dim2id, dim3id; // Allocate space for the netcdf file id int ncid; // Allocate space for the data variable ids int xid, yid, zid; Each one of these IDs will be returned through reference by the netCDF functions. While we’re at it, let’s make a variable to hold the return status of the netCDF function calls: // Allocate return status variable int retval; ## Define the meta-data Now we will start to build the netCDF file. This is a two-part process. The first part is defining the meta-data for the file and the second part is assigning the data. ### Create an empty netCDF file First, create the file: // Setup the netcdf file if((retval = nc_create("example.nc", NC_NETCDF4, &ncid))) { ncError(retval); } Note that we store the return status of the function call in retval and test the return status for an error. If there’s an error, we pass retval to our error handler. The first parameter to the function call is the name of the netCDF file. The second parameter is a flag that determines the type of netCDF file. Here we use the latest-and-greatest type of NETCDF4, which includes the HDF5/zlib compression features. If you don’t need these features, or you need a version compatible with older versions of netCDF libraries, then use the default or 64-bit offset (NC_64BIT_OFFSET) versions. The third parameter is the netCDF integer ID used for assigning variables to this file. ### Add the dimensions Now that we have a clean netCDF file to work with, let’s add the dimensions we’ll be using:  // Define the dimensions in the netcdf file if((retval = nc_def_dim(ncid, "dim1_size", dim1_size, &dim1id))) { ncError(retval); } if((retval = nc_def_dim(ncid, "dim2_size", dim2_size, &dim2id))) { ncError(retval); } if((retval = nc_def_dim(ncid, "dim3_size", dim3_size, &dim3id))) { ncError(retval); } // Gather the dimids into an array for defining variables in the netcdf file int dimids[ndims]; dimids[0] = dim1id; dimids[1] = dim2id; dimids[2] = dim3id; Just as before, we catch and test the function return status for any errors. The function nc_def_dim() takes four parameters. First is the netCDF file ID returned when we created the file. The second parameter is the name of the dimension. Here we’re using “dimX_size” – you would want to use something descriptive of this dimension (i.e. latitude, time, solution, etc.). The third parameter is the size of this dimension (i.e. number of latitude, number of solutions, etc.). The last is the ID for this dimension, which will be used in the next step of assigning variables. Note that we create an array of the dimension IDs to use in the next step. ### Add the variables The last step in defining the meta-data for the netCDF file is to add the variables: // Define the netcdf variables if((retval = nc_def_var(ncid, "x", NC_FLOAT, ndims, dimids, &xid))) { ncError(retval); } if((retval = nc_def_var(ncid, "y", NC_FLOAT, ndims, dimids, &yid))) { ncError(retval); } if((retval = nc_def_var(ncid, "z", NC_FLOAT, ndims, dimids, &zid))) { ncError(retval); } The nc_def_var() function takes 6 parameters. These include (in order) the netCDF file ID, the variable name to be displayed in the file, the type of data the variable contains, the number of dimensions of the variable, the IDs for each of the dimensions, and the variable ID (which is returned through reference). The type of data in our example is NC_FLOAT, which is a 32-bit floating point. The netCDF documentation describes the full set of data types covered. The IDs for each dimension are passed as that combined array of dimension IDs we made earlier. ### Optional: Add variable attributes This part is optional, but is incredibly useful and true to the spirit of making a netCDF file. When sharing a netCDF file, the person receiving the file should have all the information they need about the data within the file itself. This can be done by adding “attributes”. For example, let’s add a “units” attribute to each of the variables:  // OPTIONAL: Give these variables units if((retval = nc_put_att_text(ncid, xid, "units", 2, "cm"))) { ncError(retval); } if((retval = nc_put_att_text(ncid, yid, "units", 4, "degC"))) { ncError(retval); } if((retval = nc_put_att_text(ncid, zid, "units", 1, "s"))) { ncError(retval); } The function nc_put_att_text() puts a text-based attribute onto a variable. The function takes the netCDF ID, the variable ID, the name of the attribute, the length of the string of characters for the attribute, and the text associated with the attribute. In this case, we’re adding an attribute called “units”. Variable ‘x’ has units of “cm”, which has a length of 2. Variable ‘y’ has units of “degC”, which has a length of 4 (and so on). You can apply text-based attributes as shown here or numeric-based attributes using the appropriate nc_put_att_X() function (see documentation for the full list of numeric attribute functions). You can also apply attributes to dimensions by using the appropriate dimension ID or set a global attribute using the ID “0” (zero). ### End the meta-data definition portion At this point, we’ve successfully created a netCDF file and defined the necessary meta-data. We can now end the meta-data portion:  // End "Metadata" mode if((retval = nc_enddef(ncid))) { ncError(retval); } …and move on to the part 2 of the netCDF file creation process. ## Populate the file with data ### Put your data into the netCDF file Here, all we do is put data into the variables we defined in the file:  // Write the data to the file if((retval = nc_put_var(ncid, xid, &x[0][0][0]))) { ncError(retval); } if((retval = nc_put_var(ncid, yid, &y[0][0][0]))) { ncError(retval); } if((retval = nc_put_var(ncid, zid, &z[0][0][0]))) { ncError(retval); } The function nc_put_var() takes three parameters: the netCDF file ID, the variable ID, and the memory address of the start of the multi-dimensional data array. At this point, the data will be written to the variable in the netCDF file. There is a way to write to the netCDF file in data chunks, which can help with memory management, and a way to use parallel I/O for writing data in parallel to the file, but I have no experience with that (yet). I refer those interested in these features to the netCDF documentation. ### Finalize the netCDF file That’s it! We’re done writing to the netCDF file. Time to close it completely:  // Close the netcdf file if((retval = nc_close(ncid))) { ncError(retval); } ## Compile and run the code Let’s compile and run the code to generate the example netCDF file: gcc -o netcdf_example netcdf_write_example.c -lnetcdf Some common problems people run into here are not including the netCDF library flag at the end of the compilation call, not having the header files in the include-path, and/or not having the netCDF library in the library-path. Check your user environment to make sure the netCDF paths are included in your C_INCLUDE_PATH and LIBRARY_PATH: env | grep –i netcdf Once the code compiles, run it to generate the example netCDF file: ./netcdf_example If everything goes according to plan, there should be a file called “example.nc” in the same directory as your compiled code. Let’s load this up in R for some post-processing. # Reading a netCDF file in R… ## Install and load the “ncdf4” package To start using netCDF files in R, be sure to install the netCDF package “ncdf4”: install.packages("ncdf4") library(ncdf4) Note that there’s also an “ncdf” package. The “ncdf” package reads and writes the classic (default) and 64-bit offset versions of netCDF file. I recommend against using this package as the new package “ncdf4” can handle the old file versions as well as the new netCDF4 version. Turns out the “ncdf” package has been removed from the CRAN repository. It’s just as well since the new “ncdf4” package obsoletes the “ncdf” package. ## Open the netCDF file With the library installed and sourced, let’s open the example netCDF file we just created:  nc <- nc_open("example.nc") This stores an open file handle to the netCDF file. ## View summary of netCDF file Calling or printing the open file handle will produce a quick summary of the contents of the netCDF file:  print(nc) This summary produces the names of the available variables, the appropriate dimensions, and any global/dimension/variable attributes. ## Extract variables from the netCDF file To extract those variables, use the command: x <- ncvar_get(nc, "x") y <- ncvar_get(nc, "y") z <- ncvar_get(nc, "z") At this point, the data you extracted from the netCDF file are loaded into your R environment as 3-dimensional arrays. You can treat these the same as you would any multi-dimensional array of data (i.e. subsetting, plotting, etc.). Note that the dimensions are reported in reverse order from which you created the variables. dim(x) ## Close the netCDF file When you’re done, close the netCDF file: nc_close(nc) And there you have it! Hopefully this step-by-step tutorial has helped you incorporate netCDF support into your project. The code I described here is available through bitbucket. Happy computing! ~Greg # Porting GCAM onto a UNIX cluster Hi All, One of my recent research tasks has been to port GCAM (Global Change Assessment Model) on to the Cube, our groups HPC cluster, and some of the TACC resources. This turned out to be more challenging than I had anticipated, so I wanted to share my experiences in the hopes that it will save you all some time in the future. This post might be a bit pedestrian for some readers, but I hope it’s helpful to folks that are new to C++ and Linux. Before getting started, some background on GCAM. GCAM was developed by researches at the Joint Global Change Research Institute (JGCRI) at the Pacific Northwest National Lab (PNNL). GCAM is an integrated assesment model (IAM), which means it pairs a climate model with a global economic model. GCAM is unique from many IAMs in that it has a fairly sophisticated representation of various sectors of the global economy modeled on a regional level…meaning the model is a beast compared to other IAMs (like DICE). I’ll be writing a later post on IAMs more generally. GCAM is coded in C++ and makes extensive use of xml databases for both input and output. The code is open source (available here), has a Wiki (here), and a community listserve where researches can pose questions to their peers. There are three flavors of the model available: one for Windows users, one for Mac users, and one for Unix users. The Windows version comes compiled, with a nice user-interface for post-processing the results. The Unix version comes as uncompiled code, and requires the installation of some third party C++ libraries. For those who’d like to sniff around GCAM the Windows or Mac versions are a good starting point, but if you’re going to be doing HPC heavy lifting, you’ll need to work with the Unix code. The GCAM Wiki and the detailed user manual provide excellent documentation about running GCAM the Model Interface tools, but are a bit limited when describing how to compile the Unix version of the code. One helpful page on the Wiki can be found here. GCAM comes with a version of the Boost C++ library in the standard Unix download (located in Main_User_Workspace/libs). GCAM also uses the Berkeley DBXML library, which can be downloaded here. You’ll want to download version 2.5.16 to be sure it runs well with GCAM. Once you’ve downloaded DBXML, you’ll need to build the library. As it turns out, this is fairly easy. Once you’ve ported the DBXML library onto your Unix cluster, simply navigate into the main directory and run the script buildall.sh. The buildall.sh script accepts flags that allow you to customize your build. We’re going to use the -b flag to build the 64-bit version of DBXML: sh buildall.sh -b 64  once you’ve build the DBXML library, you are nearly ready to build GCAM, but must first export the paths to the Boost library and to the DBXML include and lib directories. It seems that the full paths must be declared. In other words, they should read something like: export BOOST_INCLUDE=/home/fs02/pmr82_0001/jrl276/GCAM_Haewon_trans/libs/boost_1_43_0/ export DBXML_INCLUDE=/home/fs02/pmr82_0001/jrl276/dbxml-2.5.16/install/include/ export DBXML_LIB=/home/fs02/pmr82_0001/jrl276/dbxml-2.5.16/install/lib/  One tricky point here is that the full path must be typed out instead of using ‘~’ as a short-cut (I’m not sure why). Once this is done, you are ready to compile gcam. Navigate into the ‘Main_User_Workspace’ directory and type the command:  make gcam  Compiling GCAM takes several minutes, but at the end there should be a gcam.exe executable file located in the ‘/exe/’ folder of the ‘Main_User_Workspace’. It should be noted that there is a multi-thread version of gcam (more info here), which I’ll be compiling in the next few days. If there are any tricks to that process, I’ll post again with some tips. That’s it for now. # Setting up Eclipse for C/C++ IDEs are tools to make code development a lot easier, specially if your project has multiple files, classes, and functions. However, setting up the IDE can sometimes be as painful as developing complex codes without an IDE. This post will present a short tutorial about how to install and configure Eclipse for C/C++ on Windows 7 in a (hopefully) fairly painless manner. This tutorial is sequenced as follows: 1. Installation 1. Downloading the Java Runtime Environment. 2. Downloading the GCC compiler. 3. Downloading Eclipse. 2. First steps with Eclipse 1. Setting up a template (optional) 2. Creating a new project 3. Including libraries in your project INSTALLATION Downloading the Java Runtime Environment To check if you have the Java Runtime Environment installed, go to java.com with either Internet Explorer or Firefox (Chrome will block the plugin) and click on “Do I have Java?”. Accept running all the pluggins and, If the website tells you you do not have java, you will have to download and install it from the link displayed on the website. Downloading the GCC compiler After the check is done, you will have to download the GCC compiler, which can be done from http://www.equation.com. On the side menu, there will be a link to Programming Tools, which after expanded shows a link to Fortran, C, C++. Click on this link and download the right GCC version for your system (32/64 bit), as shown in the following screenshot. After downloading it, double click on the executable, accept the licence, and type “c:\MinGW” as the installation directory. This is important because this is the first folder where Eclipse will look for the compiler in your computer. Proceed with the installation. Downloading Eclipse Now it is time to download an install eclipse. Go to the Eclipse download website and download Eclipse IDE for C/C++ Developers. Be sure to select the right option for your computer (Windows, 32bit/64bit), otherwise eclipse may not install and even if it does it will not run after installed. If unsure about which version you should download, this information can be found at Control Panel -> System by looking at System type. After downloading it, extract the file contents to “C:\Program Files\eclipse” (“Program Files (x86) if installing the 32 bits version) so that everything is organized. Note that for this you will need to start WinRAR or any other file compression program with administrative privileges. This can be done by right clicking the name of the program on the start menu and clicking on Run as Administrator. Now, go to C:\Program Files\eclipse and double click on eclipse.exe to open eclipse. In case you get an error message saying, among other things:  Java was started but returned exit code=13 ... ... -os win32 -ws win32 ...  then delete the whole eclipse folder, go back to the eclipse download page, download eclipse 32 bit, and extract it as previously described. You should not see the same error again when trying to run eclipse.exe. Now that Eclipse is up and running, it is time to use it. FIRST STEPS WITH ECLIPSE The first thing eclipse will do is ask you to choose a workspace folder. This is the folder where all your code projects will be stored. It should not matter too much which folder you choose, so using the default is probably a good idea. Setting up templates (optional) It is helpful to create a code template in order to avoid retyping the same standard piece of code every time you create a new file or project. Many scientific codes have similar imports (such as math.h and stdio.h) and all of them must have a main method (as any C++ code). If we create a code template with a few common imports and the int main function, we can just tell Eclipse when creating a new project to add these to a new .cpp file. In order to create the mentioned template, go to Window -> Preferences. There, under C/C++ -> Code Style on the left panel, click on Code Templates. Under Configure generated code and comments, expand Files -> C++ Source File, and then click on New. Choose a meaningful name for your template (I chose “Cpp with main”) and type a short description. After that, copy and paste the template below under “Pattern”. /* File:${file_name}

Author: ${user} Date:${date}
*/

#include <iostream>
#include <string>
#include <math.h>
#include <stdio.h>
#include <string.h>

using namespace std;

int main()
{
// Your code here.

return 0;
}


Note ${file_name},${data}, and \${user} are variables, which means that they will be replaced by your file’s actual data. To see a list of the other variables that can be inserted in your template, click on Insert Variable…. Click Ok and Ok again and your template will be ready to be used!

Creating a new project

Click on File -> New -> C++ Project. Under Project type choose Empty Project, then under Toolchains choose MinGW GCC, and, finally, type “project1” as your project name an click on Finish.

After your project is created, click on File -> New -> Source File. Type “say_something.cpp” (no quotes and do not forget the .cpp after the file name) as the name of your source file and choose the template you created as the template. The window should then look like this:

Click on Finish. If you used the template, replace the comment “// Your code here.” by “cout << “Yay, it worked!” << endl;”. Your code should look like the snippet below. If you have not created the template, just type the following code to your file.

/*
File: say_something.cpp

Author: bct52
Date: Jun 26, 2015
*/

#include <iostream>
#include <string>
#include <math.h>
#include <stdio.>
#include <string.h>

using namespace std;

int main()
{
cout << "Yay, it worked!" << endl;

return 0;
}


Now, build the code by clicking on the small hammer above the code window and, after the project is built, click on the run button (green circle with white play sign in the center). If everything went well, your window should look like the screenshot below, which means your code compiled and is runs as expected.

Including libraries in your project

When developing code, often times other people have had to develop pieces of code to perform some of the intermediate steps we want our code to perform. These pieces of code are often publicly available in the form of libraries. Therefore, instead of reinventing the wheel, it may be better to simply use a library.

Some libraries are comprised of one or a few files only, and can be included in a project simply by dragging the file into the Eclipse project. Others, however, are more complex and should be installed in the computer and then called from the code. The procedure for the latter case will be described here, as it is the most general case . The process of installation and usage of the Boost library with MinGW (GCC) will be used here as a case study.

The first step is downloading the library. Download the Boost library from here and extract it anywhere in your computer, say in C:\Users\my_username\Downloads (it really doesn’t matter where because these files will not be used after installation is complete).

Now it is time to install it. For this:

1. Hold the Windows keyboard button and press R, type “cmd”, and press enter.
2. On the command prompt, type “cd C:\Users\bct52\Downloads\boost_1_58_0” (or the directory where you extracted boost to) and press enter.
3. There should be a file called bootstrap.bat in this folder. If that is the case, run the command:
bootstrap.bat mingw
4. In order to compile Boost to be used with MinGW, compile Boost with the gcc toolset. You will have to choose an installation directory for Boost, which WILL NOT be the same directory where you extracted the files earlier. In my case, I used C:\boost. For this, run the command:
b2 install --prefix=C:\boost toolset=gcc

Now go read a book or work on something else because this will take a while.

Now, if the installation worked with just warnings, it is time to run a code example from Boost’s website that, or course, uses the Boost library. Create a new project called “reveillon” and add a source file to it called “days_between_new_years.cpp” following the steps from the “Creating a new project” section. there is no need to use the template this time.

You should now have a blank source file in front of you. If not, delete any text/comments/codes in the file so that the file is blank. Now, copy and paste the following code, from Boost’s example, into your file.

 /* Provides a simple example of using a date_generator, and simple
* mathematical operatorations, to calculate the days since
* New Years day of this year, and days until next New Years day.
*
* Expected results:
* Adding together both durations will produce 366 (365 in a leap year).
*/
#include <iostream>
#include "boost/date_time/gregorian/gregorian.hpp"

int
main()
{

using namespace boost::gregorian;

date today = day_clock::local_day();
partial_date new_years_day(1,Jan);
//Subtract two dates to get a duration
days days_since_year_start = today - new_years_day.get_date(today.year());
std::cout << "Days since Jan 1: " << days_since_year_start.days()
<< std::endl;

days days_until_year_start = new_years_day.get_date(today.year()+1) - today;
std::cout << "Days until next Jan 1: " << days_until_year_start.days()
<< std::endl;
return 0;
};


Note that line 9 (“#include “boost/date_time/gregorian/gregorian.hpp””) is what tells your code what exactly is being used from Boost in your code. Line 15 (“using namespace boost::gregorian;”) saves you from having to type boost::gregorian every time you want to use one of its functions.

However, the project will still not compile in Eclipse because Eclipse still does not know where to look for the Boost library. This will require a couple of simple steps:

1. Right click on the project (reveillon), under the Project Explorer side window, then click on Properties. Under C/C++ Build->Settings, click on Includes under GCC C++ Compiler. On the right there should be two blank boxes, the top one called Include paths (-I) and the other called Include files (-include). Under Include paths (top one), add the path “C:\boost\include\boost-1_58” (note that this path must reflect the path where you installed Boost as well as which version of Boost you have). This is where the compiler will look for the header file specified in the code with the #include statement.
2. The compiled library files themselves must be included through the linker. This step is necessary only if you are using a compiled library. For this, on the same window, click on Libraries under MinGW C++ Linker. Add the path to the Boost libraries folder to the Library search path (-L) (bottom box). this path will be “C:\boost\lib” (again, if you installed Boost in a different folder your path will be slightly different). Now the actual compiled library must be added to the Libraries (-i) (top box). First, we need to figure out the name of the compiled library file used in the code. In this case, it is the file “libboost_date_time-mgw51-mt-d-1_58.a”. Therefore, add boost_date_time-mgw51-mt-d-1_58 (no lib prefix, no .a postfix, and be sure to match the name of your file) to Libraries (-i). Click Ok and Ok again.

Now compile the code by clicking on the hammer button and run the rode by clicking on the play button. Below is a screenshot reflecting both steps above as well as the expected output after running the program.

That’s it. After your model is in a good shape and it is time to run it with Borg (or other optimization algorithm), just change your “int main()” to a function with your model’s name and the right Borg’s arguments, add the standard Borg main, and change the makefile accordingly. Details on how to do all this for Borg will be explained in a future post.