Introduction
I’ve implemented dynamic simulation models for chemical engineering applications for more than 50 years now, and am still doing it. You would think it should have gotten pretty old by now, but truth is, it is still fascinating to me. The model principles have not changed but the tools to get them into a computer have.
In this post I will take you on a 50 year journey through a language landscape I personally visited. I’m not covering every language I’ve used along the way but focussing on two old ones; FORTRAN and C, one intermediate, Java, and three for me new ones; Swift, Python and Kotlin. The examples will give me an opportunity to highlight some problems with the old languages later solved by the newer ones.
The example I’ll use for this post is how to implement an interpolation method based on Cubic Splines. You can easily find explanations of what Cubic Splines are but this brief intro is good for starters: https://mathworld.wolfram.com/CubicSpline.html
Admittedly, I don’t use splines very often in my dynamic modeling work but I chose these because both FORTRAN and C codes are well documented (see Forsythe et. al. Computer Methods for Mathematical Computations Prentice-Hall, Inc., 1977 and Press et. al. Numerical Recipes in C 2nd ed., Cambridge University Press 1992). For both the FORTRAN and C codes I have typed in the code from the two books. Since I don’t have a FORTRAN compiler anymore, I have not been able to check that code directly. I do have C compilers and verified that the C code works as expected. The other codes are my own translations from the original FORTRAN code. They also work as they should.
FORTRAN
Back in the 70’s when I started computer modeling in all earnest, FORTRAN was pretty much the only game in town for scientific and engineering calculations. Code was entered on punch cards (I had boxes full of them) and ran on multi-million dollar mainframe computers that were slower than my current iPhone or iPad!
As you can tell from the code fragments below, the formatting was very strict; only one line per card, short variable names (often just a single letter), code had to start at the 7th column (as I recall), comment lines started with a capital C in the left most column, etc.
Awkward as it was to work with punch cards and slow computers, I have to admit that FORTRAN as a language was pretty easy to learn and work with. It had a very limited syntax and no fancy constructs or external libraries to remember. And certainly no pointers or memory management to worry about (more on that later). And you have to agree with me, the code is very readable. All caps, each row starting at the same place, no squiggly brackets or semicolons, just plain text. As you see the code snippet of C below, you will appreciate the clean looking nature of FORTRAN.
But there is one nasty feature that I just have to mention and that is the GO TO statement and the line numbers. For example, in the first subroutine, a few code lines down, there is a statement “IF ( N .LT. 3 ) GO TO 50“. It says that if there are less than 3 data points available, go to the line in the code with the label “50”. That is towards the end of the subroutine. I used to joke with my friends and say that the GO TO is not the problem, it’s the “WHERE FROM” that causes all the headache. To understand what I’m talking about you have to image sifting through paper printouts of hundreds of lines of code and seeing some lines with a number next to them. If the number was not followed by the word CONTINUE you had to suspect that it was part of a potentially distant GO TO statement. Debugging was hard as you basically had to manually inspect the code and try to figure out what was wrong. The best medicine was to get it right the first time around. Also, breaking the code into smaller functions and subroutines was also a good way of isolating problems and avoid GO TO’s.
SUBROUTINE SPLINE(N, X, Y, B, C, D)
INTEGER N
REAL X(N), Y(N), B(N), C(N), D(N)
C THE COEFFICIENTS B(I), C(I), AND D(I), I=1,2,...,N ARE COMPUTED
C FOR A CUBIC INTERPOLATING SPLINE
C
C S(X) = Y(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2 + D(I)*(X-X(I))**3
C FOR X(I) .LE. X .LE. X(I+1)
C
C INPUT..
C
C N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2)
C X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER
C Y = THE ORDINATES OF THE KNOTS
C
C OUTPUT..
C
C B, C, D = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE
C
C USING P TO DENOTE LEVEL OF DIFFERENTIATION
C
C Y(I) = S(X(I))
C B(I) = SP(X(I))
C C(I) = SPP(X(I))/2
C D(I) = SPPP(X(I))/6 (DERIVATIVE FROM THE RIGHT)
C
C THE ACCOMPANYING FUNCTION SUBPROGRAM SEVAL CAN BE USED
C TO EVALUATE THE SPLINE
C
INTEGER NM1, IB, I
REAL T
C
NM1 = N-1
IF ( N .LT. 2 ) RETURN
IF ( N .LT. 3 ) GO TO 50
C
D(1) = X(2) - X(1)
C(2) = (Y(2) - Y(1))/D(1)
DO 10 I = 2, NM1
D(I) = X(I+1) - X(I)
B(I) = 2.*(D(I-1) + D(I))
C(I+1) = (Y(I+1) - Y(I))/D(I)
C(I) = C(I+1) - C(I)
10 CONTINUE
C
B(1) = -D(1)
B(N) = -D(N-1)
C(1) = 0.
C(N) = 0.
IF ( N .EQ. 3 ) GO TO 15
C(1) = C(3)/X(4)-X(2)) - C(2)/(X(3)-X(1))
C(N) = C(N-1)/(X(N)-X(N-2)) - C(N-2)/(X(N-1)-X(N-3))
C(1) = C(1)*D(1)**2/(X(4)-X(1))
C(N) = -C(N)*D(N-1)**2/(X(N)-X(N-3))
C
C FORWARD ELIMINATION
C
15 DO 20 I = 2, N
T = D(I-1)/B(I-1)
B(I) = B(I) - T*D(I-1)
C(I) = C(I) - T*C(I-1)
20 CONTINUE
C
C BACK SUBSTITUTION
C
C(N) = C(N)/B(N)
DO 30 IB = 1, NM1
I = N-IB
C(I) = (C(I) - D(I)*C(I+1))/B(I)
30 CONTINUE
C
C COMPUTE POLYNOMIAL COEFFICIENTS
C
B(N) = (Y(N) - Y(NM1))/D(NM1) + D(NM1)*(C(NM1) + 2.0*C(N))
DO 40 I = 1, NM1
B(I) = (Y(I+1) - Y(I))/D(I) - D(I)*(C(I+1) + 2.0*C(I))
D(I) = (C(I+1) - C(I))/D(I)
C(I) = 3.*C(I)
40 CONTINUE
C(N) = 3.*C(N)
D(N) = D(N-1)
RETURN
C
50 B(1) = (Y(2)-Y(1))/X(2)-X(1))
C(1) = 0.
D(1) = 0.
B(2) = B(1)
C(2) = 0.
D(2) 0.
RETURN
END
REAL FUNCTION SEVAL(N, U, X, Y, B, C, D)
INTEGER N
REAL U, X(N), Y(N), B(N), C(N), D(N)
C
C THIS SUBROUTINE EVALUATES THE CUBIC SPLINE FUNCTION
C
C SEVAL = Y(I) + B(I)*(U-X(I)) + C(I)*(U-X(I))**2 + D(I)*(U-X(I))**3
C
INTEGER I, J, K
REAL DX
DATA I/1/
IF ( I .GE. N) I = 1
IF ( U .LT. X(1) ) GO TO 10
IF ( U .LE. X(I+1) ) GO TO 30
C
C BINARY SEARCH
C
10 I = 1
J = N + 1
20 K = (I + J)/2
IF ( U .LT. X(K) ) J = K
IF ( U .GE. X(K) ) I = K
IF ( J .GT. I+1 ) GO TO 20
C
C EVALUATE SPLINE
C
30 DX = U - X(I)
SEVAL = Y(I) + DX*(B(I) + DX*(C(I) + DX*D(I)))
RETURN
ENDC
While FORTRAN was the best choice for scientific calculations on a mainframe, it was not useful for writing new operating systems or drivers for peripheral equipment. Those tasks were done in ASSEMBLY language requiring great programming skills. A better language was needed as the computer industry grew and different architectures became available. During the course of developing the operating system, UNIX, the engineers at Bell Labs devised a new, higher level language they named “C”. Just as FORTRAN was very restrictive as to how you could manipulate a computer, C is quite the opposite. We used to joke that when programming in C the language gives you a lot freedom and enough rope “to hang yourself”. The two biggest problems with C are the liberal use of pointers and the need for manual memory management. Both of those issues are illustrated in the two code snippets below.
The first code snippet defines the various functions used to compute the interpolating splines. The second snippet shows how to use these functions. Again, these codes are typed in from the Numerical Recipes in C book and are not using exactly the same algorithm as in the FORTRAN code. That’s ok because I’m just illustrating some significant differences between the two languages.
#include <stdio.h>
#include <stdlib.h>
#define NR_END 1
#define FREE_ARG char*
void nrerror(char error_text[])
{
fprintf(stderr,"%s\n", error_text);
exit(1);
}
float *vector(long nl, long nh)
{
float *v;
v=(float *)malloc((size_t) ((nh-nl+NR_END)*sizeof(float)));
if (!v) nrerror("allocation failure in vector()");
return v-nl+NR_END;
}
void free_vector(float *v, long nl, long nh)
{
free((FREE_ARG) (v+nl-NR_END));
}
void spline(const float x[], const float y[], int n, float yp1, float ypn, float y2[])
{
int i,k;
float p,qn,sig,un,*u;
u=vector(1,n-1);
if (yp1 > 0.99e30)
y2[1]=u[1]=0.0f;
else {
y2[1] = -0.5f;
u[1]=(3.0f/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
}
for (i=2;i<=n-1;i++) {
sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
p=sig*y2[i-1]+2.0f;
y2[i]=(sig-1.0f)/p;
u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
u[i]=(6.0f*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
}
if (ypn > 0.99e30)
qn=un=0.0f;
else {
qn=0.5f;
un=(3.0f/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
}
y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0f);
for (k=n-1;k>=1;k--)
y2[k]=y2[k]*y2[k+1]+u[k];
free_vector(u,1,n-1);
}
void splint(const float xa[], const float ya[], const float y2a[], int n, float x, float *y)
{
int klo,khi,k;
float h,b,a;
klo=1;
khi=n;
while (khi-klo > 1) {
k=(khi+klo) >> 1;
if (xa[k] > x) khi=k;
else klo=k;
}
h=xa[khi]-xa[klo];
if (h == 0.0f) nrerror("bad xa input to routine splint");
a=(xa[khi]-x)/h;
b=(x-xa[klo])/h;
*y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0f;
}Let’s first focus on the function “float *vector(long nl, long nh)” above. It claims that, given a starting index (nl) and an ending index (nh), the function will return a pointer to the beginning of a stretch of memory containing a string of real numbers of type float. The vector is allocated on the computer’s memory called the heap. This memory has to be reclaimed later or else we can run into serious issues during run time. The freeing of memory is done with the function “void free_vector(float *v, long nl, long nh)“. This function takes the pointer to the start of the memory that needs to be cleared along with the size of the vector.
In the application example below we see how we first allocate space for the vectors xa, ya, and y2 depending on how many data points, NP, we have. Inside the spline function we locally have another variable, u, that also needs memory before it can be used. Failing to allocate memory to a vector before it is used leads to a crash and failing to remove memory after it has been used can also eventually lead to a crash. Both of these bugs are hard to find.
But don’t get me wrong, despite its flaws, C was a great advance in computer science. It is incredibly flexible, terse, compact, and very fast. It is still used today for CPU and memory critical applications and is often the foundation for other languages (e.g. Python).
Before we leave the C code, I’d like you to pay attention to the reduced readability of the code compared to FORTRAN. The lower case is easier to type on a keyboard but harder to read, I think. The formatting is loose and optional since individual statements end with a semicolon and blocks of code are surrounded by squiggly brackets. Those features make it clear to the computer what to do but do not make it particularly clear to a human reading the code. There are good reasons for this. Back 40-50 years ago when computers were slow, CPU time was the limitation and code needed to be entered as efficiently as possible. Today computers are fast and “people time” is the limitation. There is much more emphasis on code clarity and readability as you will see in the modern language examples below.
#include <stdio.h>
#include <math.h>
#include "spline.h"
#define NP 10
#define PI 3.1415926
int main() {
int i,nfunc;
float f,x,y,yp1,ypn,*xa,*ya,*y2;
xa=vector(1,NP);
ya=vector(1,NP);
y2=vector(1,NP);
for (nfunc=1;nfunc<=2;nfunc++) {
if (nfunc == 1) {
printf("\nsine function from 0 to pi\n");
for (i = 1; i < NP; i++) {
xa[i] = (float) (i * PI / NP);
ya[i] = (float) (sin(xa[i]));
}
yp1 = (float) cos(xa[1]);
ypn = (float) cos(xa[NP]);
} else if (nfunc == 2) {
printf("\nexponential function from 0 to 1\n");
for (i = 1; i < NP; i++) {
xa[i] = (float) (i*1.0 / NP);
ya[i] = (float) (exp(xa[i]));
}
yp1 = (float) exp(xa[1]);
ypn = (float) exp(xa[NP]);
} else {
break;
}
spline(xa, ya, NP, yp1, ypn, y2);
for (i = 1; i <= 10; i++) {
if (nfunc == 1) {
x = (float) ((-0.05 + i / 10.0) * PI);
f = (float) sin(x);
} else if (nfunc == 2) {
x = (float) (-0.05 + i / 10.0);
f = (float) exp(x);
}
splint(xa, ya, y2, NP, x, &y);
printf("%12.6f %12.6f %12.6f\n", x, f, y);
}
}
free_vector(y2,1,NP);
free_vector(ya,1,NP);
free_vector(xa,1,NP);
return 0;
}Object-Oriented Programming – Java
Alongside advances in hardware and the advent of the personal computer in the early 80’s, there was also a quest to organize programs around chunks of data rather than just using collections of subroutines and functions operating on banks of common data. The idea is this: instead of having a library of public functions that use public data as parameters (and also modify the same common data), how about focussing on separate chunks of data (objects) that could have special functions (methods) associated with them to modify their data? Data in one object are not shared with other objects unless specifically requested. This level of encapsulation proved to be very beneficial for productivity and error reduction. Logical errors, that is, since a good editor and development environment could point out syntax errors for you.
FORTRAN and C are both procedural languages and have no built-in mechanisms for object-orientation. A solution to this problem came from Bjarne Stroustrup at Bell Labs when he invented C++. The new language has often been called a “better C”, but it is more than that. It has natural ways of creating classes and creating objects from them. While pointers and memory management are still present in the language they are greatly simplified for the programmer. Best of all, C++ is compiled and has the same great speed as C.
Other object-oriented languages came on the scene roughly at the same time or shortly after C++ was released. These were Smalltalk, Objective-C and Java. They completely eliminated the need for the programmer to worry about memory management and they also removed direct use of pointers. In addition, they were strictly object-oriented (everything became an object) and the notion of stand alone functions and functional programming was minimized or removed.
I personally switched from C and C++ to Java in the mid- to late 90’s and productively used that language for over 15 years. Before Python became firmly established as a language of choice, Java was likely the most popular language of its day. Java is incredibly versatile, while also rigorous, and lends itself to large stand-alone applications. I have built numerous simulation models that were used in my own research but could also easily be shared with my customers. The saying “Write Once, Run anywhere” rings true for Java in my own experience.
To complete the discussion of Java I show the code snippet that implements the spline functions, or the “SplineGenerator” class as it were in case of Java.
This code illustrates many of the points of object-oriented programming I made earlier. First, we see that the class contains all the vectors it needs and it protects them (keyword private) from being accessed from the outside other than through specific getter methods (e.g. public double[] getUX()). Memory for these vectors is allocated by the “new double[n]” statement, which hides the cumbersome pointer code used in the C example. Furthermore, there is no need to explicitly deallocate this memory, it’s done for you behind the scenes by the so-called garbage collector. Without pointers in the program the syntax looks more like FORTRAN than C. But it is much more expressive than FORTRAN. Notice for example how elegantly we get around the GO TO’s by using a switch statement in Java. All three cases of n = 1, 2 or >=3 are easily handled in one place without a need to jump around in the code.
A final word around code readability. You notice that the semicolons and the squiggly brackets are still in the Java code. From that standpoint formatting is optional but a good code editor helps keep things formatted and in order. I don’t know for sure but I could suspect that at the time when Java was invented there was still a need to make the code clear to the compiler rather than making the code easy to read for humans. But Java code is easier to read than C code, that’s for sure.
public class SplineGenerator {
private double[] b;
private double[] c;
private double[] d;
private double[] ux;
private double[] uxx;
public void spline(double[] x, double[] y) {
int n = x.length;
if (n != y.length) {
System.out.println("x and y are not the same length in spline!!!");
return;
}
b = new double[n];
c = new double[n];
d = new double[n];
switch (n) {
case 1:
System.out.println("Arrays must have more than one data points in spline!!!");
break;
case 2:
double b0 = (y[1] - y[0]) / (x[1] - x[0]);
b[0] = b0;
b[1] = b0;
break;
default:
int nm1 = n - 1;
d[0] = x[1] - x[0];
c[1] = (y[1] - y[0]) / d[0];
for (int i = 1; i < nm1; i++) {
d[i] = x[i + 1] - x[i];
b[i] = 2 * (d[i - 1] + d[i]);
c[i + 1] = (y[i + 1] - y[i]) / d[i];
c[i] = c[i + 1] - c[i];
}
b[0] = -d[0];
b[nm1] = -d[nm1 - 1];
c[0] = 0;
c[nm1] = 0;
if (n > 3) {
c[0] = c[2] / (x[3] - x[1]) - c[1] / (x[2] - x[0]);
c[nm1] = c[nm1 - 1] / (x[nm1] - x[nm1 - 2]) - c[nm1 - 2] / (x[nm1 - 1] - x[nm1 - 3]);
c[0] = c[0] * d[0] * d[0] / (x[3] - x[0]);
c[nm1] = -c[nm1] * d[nm1 - 1] * d[nm1 - 1] / (x[nm1] - x[nm1 - 3]);
}
//forward elimination
for (int i = 1; i < n; i++) {
double t = d[i - 1] / b[i - 1];
b[i] -= t * d[i - 1];
c[i] -= t * c[i - 1];
}
//back substitution
c[nm1] /= b[nm1];
for (int ib = 0; ib < nm1; ib++) {
int i = n - 2 - ib;
c[i] = (c[i] - d[i] * c[i + 1]) / b[i];
}
//compute polynomial coefficients
b[nm1] = (y[nm1] - y[nm1 - 1]) / d[nm1 - 1] + d[nm1 - 1] * (c[nm1 - 1] + 2 * c[nm1]);
for (int i = 0; i < nm1; i++) {
b[i] = (y[i + 1] - y[i]) / d[i] - d[i] * (c[i + 1] + 2 * c[i]);
d[i] = (c[i + 1] - c[i]) / d[i];
c[i] = 3 * c[i];
}
c[nm1] = 3 * c[nm1];
d[nm1] = d[nm1 - 1];
}
}
public double evaluateSplineAt(double u, double[] x, double[] y) {
int i = 0;
int j = x.length;
do {
int k = (i + j) / 2;
if (u < x[k]) {
j = k;
}
if (u >= x[k]) {
i = k;
}
} while (j > i + 1);
double dx = u - x[i];
return y[i] + dx * (b[i] + dx * (c[i] + dx * d[i]));
}
public void derivativesAtGridpoints(double xl, double xu, int n, double[] u, double percent) {
//compute unequal grid spacing with percent increase
double[] x = new double[n];
ux = new double[n];
uxx = new double[n];
double sum = 0.0;
for (int i = 0; i < (n - 1); i++) {
sum += Math.pow((1.0 + percent / 100), i);
}
double dx = (xu - xl) / sum;
for (int i = 0; i < n; i++) {
if (i == 0) {
x[i] = xl;
} else {
x[i] = x[i - 1] + dx * Math.pow((1.0 + percent / 100), (i - 1));
}
}
spline(x, u);
for (int i = 0; i < n; i++) {
ux[i] = b[i];
uxx[i] = 2.0 * c[i];
}
}
public double[] getUX() {
return ux;
}
public double[] getUXX() {
return uxx;
}
}Swift
When I started SySim Consulting in 2014 I transitioned from Java to Objective-C and later to Swift. The reason was partly that I wanted to try something new and partly that I became interested in exploring the iPad as a tool for interactive simulations.
I found Objective-C to be unintuitive and quite difficult to learn and use. I must not have been alone since with Apple’s new language, Swift, they took a whole new approach to language syntax, much more like C++ and Java.
Swift is a very modern language drawing from the experiences of programmers who had struggled with the limitations and quirks found in older paradigms. It took me a couple of years of using the language both for Mac OS and iOS applications before I got completely comfortable with it. The learning curve for me was extra steep because the language kept evolving as I was learning. New features were introduced with frequent updates making it hard for me to keep up.
There are too many novel features in Swift to cover here but I will mention three concepts that were new to me and I found very powerful: (1) Protocol extensions, (2) Optionals, and (3) Tuples. The first one has to do with objects and is not illustrated in the code snippet below. But the idea of Optionals and Tuples are.
The use of Optionals is shown in the return type of the spline function as -> coefficientVectors? The question mark after the struct expresses the notion that the spline function optionally may not return anything. In fact, the code shows that if the x and y input vectors are of different lengths we cannot proceed and return with nothing (nil). Also if n = 1 we must exit without the struct and also return nil.
How do we then deal with the spline function when we call it? Built into the Swift language are a couple of different ways to recognize an Optional. If you look towards the end of the code snippet you find a line that says: “if let coefVectors = spline(xArray: x, yArray: u)“. In words this says that if the spline function returns a struct then assign it to the variable coefVectors, otherwise skip the code that follows the if. That way we are assured not to crash from pointing at empty space.
A Tuple is a collection of variables, not necessarily of the same type as in an array. There are often times when you like to return more than one value from a function or a class method. FORTRAN, C or Java do not permit this. In C you have to go through the function’s argument list and pass pointers to the variables you want to consume or produce. In Java you need to form a new object that contains the variables and return this object or provide getter functions to the variables of interest. The use of Tuples in Swift is far more elegant and transparent. Both the second and third functions below return Tuples. No new objects required, just put what you want to return in parentheses and return it!
I find Swift to be the most elegant programming language I’ve used to date. The strongly typed syntax, the many modern constructs and the fast execution speed makes it ideal for building major applications. I’ve do so for several different processes and implemented the code both on Mac’s (MacOS) and iPads (iOS). It works. And I have to mentioned the readability of the code. Notice that the semicolons are gone. Gone are also parentheses around for loops and if-statements. The squiggly brackets are still there but you’ll see in the Python code even those are gone. In short, a modern language emphasizes human readability.
Well, nothing is perfect, at least not yet. While Swift is now open source https://www.swift.org it is still very much tied to Apple devices and their many support libraries. And I discovered in talks and presentations I’ve given that the engineering community is more Windows based than MacOS oriented. So unlike with Java, where I could easily share my models, the Swift models are still confined to my Mac’s and iPads. Not ideal.
import Foundation
public struct coefficientVectors {
var b: [Double]
var c: [Double]
var d: [Double]
}
public func spline(xArray x: [Double], yArray y: [Double]) -> coefficientVectors? {
let n = x.count
if n != y.count {
return nil
}
switch n {
case 1:
return nil
case 2:
let b0 = (y[1] - y[0]) / (x[1] - x[0])
return coefficientVectors(b: [b0, b0], c: [0.0, 0.0], d: [0.0, 0.0])
default:
var b = Array(repeating: 0.0, count: n)
var c = b
var d = b
let nm1 = n - 1
d[0] = x[1] - x[0]
c[1] = (y[1] - y[0]) / d[0]
for i in 1..<nm1 {
d[i] = x[i+1] - x[i]
b[i] = 2 * (d[i-1] + d[i])
c[i+1] = (y[i+1] - y[i]) / d[i]
c[i] = c[i+1] - c[i]
}
b[0] = -d[0]
b[nm1] = -d[nm1-1]
c[0] = 0
c[nm1] = 0
if n > 3 {
c[0] = c[2] / (x[3] - x[1]) - c[1] / (x[2] - x[0])
c[nm1] = c[nm1-1] / (x[nm1] - x[nm1-2]) - c[nm1-2] / (x[nm1-1] - x[nm1-3])
c[0] = c[0] * d[0] * d[0] / (x[3] - x[0])
c[nm1] = -c[nm1] * d[nm1-1] * d[nm1-1] / (x[nm1] - x[nm1-3])
}
//forward elimination
for i in 1..<n {
let t = d[i-1] / b[i-1]
b[i] -= t * d[i-1]
c[i] -= t * c[i-1]
}
//back substitution
c[nm1] /= b[nm1]
for ib in 0..<nm1 {
let i = n - 2 - ib
c[i] = (c[i] - d[i] * c[i+1]) / b[i]
}
//compute polynomial coefficients
b[nm1] = (y[nm1] - y[nm1-1]) / d[nm1-1] + d[nm1-1] * (c[nm1-1] + 2 * c[nm1])
for i in 0..<nm1 {
b[i] = (y[i+1] - y[i]) / d[i] - d[i] * (c[i+1] + 2 * c[i])
d[i] = (c[i+1] - c[i]) / d[i]
c[i] = 3 * c[i]
}
c[nm1] = 3 * c[nm1]
d[nm1] = d[nm1-1]
return coefficientVectors(b: b, c: c, d: d)
}
}
public func evaluateSpline(at u: Double, xArray x: [Double], yArray y: [Double], coefficientVectors cv:coefficientVectors) -> (s: Double, sx: Double, sxx: Double) {
var i = 0
var j = x.count
repeat {
let k = (i + j) / 2
if u < x[k] { j = k }
if u >= x[k] { i = k }
} while j > i + 1
let dx = u - x[i]
let s = y[i] + dx * (cv.b[i] + dx * (cv.c[i] + dx * cv.d[i]))
let sx = cv.b[i] + dx * (2 * cv.c[i] + dx * 3 * cv.d[i])
let sxx = 2 * cv.c[i] + dx * 6 * cv.d[i]
return (s, sx, sxx)
}
public func derivativesAtGridpoints(lowerXBoundary xL: Double, upperXBoundary xU: Double, gridPoints n: Int, valuesAtGridPoints u: [Double], percentIncrease percent: Double = 0) -> (ux: [Double], uxx: [Double]) {
var ux = u
var uxx = u
var x = u
//compute unequal grid spacing with percent increase
var sum = 0.0
for i in 0..<(n-1) {
sum += pow((1.0 + percent / 100.0), Double(i))
}
let dx = (xU - xL) / sum
for i in 0..<x.count {
if i == 0 {
x[i] = xL
} else {
x[i] = x[i-1] + dx * pow((1.0 + percent / 100.0), Double(i-1))
}
}
if let coefVectors = spline(xArray: x, yArray: u) {
ux = coefVectors.b
for i in 0..<n {
uxx[i] = 2.0 * coefVectors.c[i]
}
}
return (ux, uxx)
}Python
Prior to giving it a serious look I had heard a lot about Python and its use in developing websites, task automation, data analysis, data visualization and machine learning but was not sure it was suitable for dynamic modeling and interactive simulations. I was skeptical of its performance since I knew it was a scripting language and not directly compiled into machine code. But if you have read some of my other posts you know that I now have embraced Python completely. In my blog post comparing Python and Swift you find many of the reasons why I like Python. Here I will merely add a few observations around the code itself specifically related to what I said around some of previously mentioned languages.
If you look at the Python code below for implementing the cubic splines we have practically come full circle around to FORTRAN clarity! Just like in the FORTRAN code there are some useful comments about how use of the code. However, the Python comments are much more than a service to the programmer, they also help the user. For example, when I have the environment open where the spline function lives and I type help(spline), I get the output shown right below. It gives me all the information for using the function without having to look at the code.
Help on function spline in module symods.integrators.utils.spline:
spline(x, y)
Calculate the cubic spline for a set of x- and y-values
Parameters
----------
x : np.array of x-values
y : np.array of y-values
Returns
-------
tuple of np.arrays of polynomial coefficients b, c, dThe next feature of code clarity is the avoidance of semicolons, squiggly brackets, and redundant parentheses. Instead, Python relies on indentation to tell the interpreter where one block ends and another starts. In addition to being a language requirement, the indentations help structure the code and make it very readable.
But as I have said before, nothing is perfect. I miss the Optionals and Switch statements from Swift along with Interfaces and Protocols from Java and Swift. Python has Tuples though and they are great. Python is slower than both Java and Swift but shares the advantage with Java that it can run pretty much anywhere. And that is important for the work I do.
import numpy as np
def spline(x, y):
"""Calculate the cubic spline for a set of x- and y-values
Parameters
----------
x : np.array of x-values
y : np.array of y-values
Returns
-------
tuple of np.arrays of polynomial coefficients b, c, d
"""
n = len(x)
if n == 1:
print('Dimension must be greater than 1')
return [0], [0], [0]
if n == 2:
b0 = (y[1] - y[0]) / (x[1] - x[0])
return [b0, b0], [0.0, 0.0], [0.0, 0.0]
b = np.zeros(n)
c = np.zeros(n)
d = np.zeros(n)
nm1 = n - 1
d[0] = x[1] - x[0]
c[1] = (y[1] - y[0]) / d[0]
for i in range(1, nm1):
d[i] = x[i+1] - x[i]
b[i] = 2 * (d[i-1] + d[i])
c[i+1] = (y[i+1] - y[i]) / d[i]
c[i] = c[i+1] - c[i]
b[0] = -d[0]
b[nm1] = -d[nm1-1]
c[0] = 0
c[nm1] = 0
if n > 3:
c[0] = c[2] / (x[3] - x[1]) - c[1] / (x[2] - x[0])
c[nm1] = c[nm1-1] / (x[nm1] - x[nm1-2]) - c[nm1-2] / (x[nm1-1] - x[nm1-3])
c[0] = c[0] * d[0] * d[0] / (x[3] - x[0])
c[nm1] = -c[nm1] * d[nm1-1] * d[nm1-1] / (x[nm1] - x[nm1-3])
# forward elimination
for i in range(1, n):
t = d[i-1] / b[i-1]
b[i] -= t * d[i-1]
c[i] -= t * c[i-1]
# back substitution
c[nm1] /= b[nm1]
for ib in range(0, nm1):
i = n - 2 - ib
c[i] = (c[i] - d[i] * c[i+1]) / b[i]
# compute polynomial coefficients
b[nm1] = (y[nm1] - y[nm1-1]) / d[nm1-1] + d[nm1-1] * (c[nm1-1] + 2 * c[nm1])
for i in range(0, nm1):
b[i] = (y[i+1] - y[i]) / d[i] - d[i] * (c[i+1] + 2 * c[i])
d[i] = (c[i+1] - c[i]) / d[i]
c[i] = 3 * c[i]
c[nm1] = 3 * c[nm1]
d[nm1] = d[nm1-1]
return b, c, d
def evaluate_spline_at(u, x, y, cv):
"""Evaluate cubic spline at an arbitrary input point u.
Parameters
----------
u : float
x : np.array of x-values
y : np.array of y-values
cv : tuple of polynomial coefficient arrays, (b, c, d) received from a call to spline.
Returns
-------
tuple of three floats, (s, sx, sxx)
where s = the y-value at u
sx = the first derivative at u
sxx = the second derivative at u
"""
i = 0
j = len(x)
b, c, d = cv
while j > i + 1:
k = (i + j) // 2
if u < x[k]:
j = k
if u >= x[k]:
i = k
dx = u - x[i]
s = y[i] + dx * (b[i] + dx * (c[i] + dx * d[i]))
sx = b[i] + dx * (2 * c[i] + dx * 3 * d[i])
sxx = 2 * c[i] + dx * 6 * d[i]
return s, sx, sxx
def derivatives_at_gridpoints(x_lower, x_upper, n, y, percent=0.0):
"""Calculate the first and second derivatives at the n gridpoints that may be unequally spaced
Parameters
----------
x_lower : float
lower bound of the x-range
x_upper : float
upper bound of the x-range
n : int
number of gridpoints
y : np.array of y-values at the gridpoints
percent : float
(Default value = 0.0)
The spacing of the n gridpoints will increase by percent (%) amount for each subsequent gridpoint
Returns
-------
tuple of two np.arrays (ux, uxx),
containing the first and second derivative vectors at the gridpoints
"""
# compute unequal grid spacing with percent increase
x = np.zeros(n)
uxx = np.zeros(n)
a_sum = 0.0
for i in range(0, n-1):
a_sum += pow((1.0 + percent / 100.0), i)
dx = (x_upper - x_lower) / a_sum
for i in range(0, len(x)):
if i == 0:
x[i] = x_lower
else:
x[i] = x[i-1] + dx * pow((1.0 + percent / 100.0), i)
b, c, d = spline(x, y)
ux = b
for i in range(0, n):
uxx[i] = 2.0 * c[i]
return ux, uxxKotlin
By now I’ve been working with Python for well over a year and found it to be a great tool. It shines when it comes to data analysis, visualization and machine learning. However, due to its interpreted nature it is a bit slow for substantial dynamic simulations.
So, I started to look back at Java to see how that language had progressed since I left it in 2014. In my search I came across a really good blog post https://www.marcobehler.com/guides/a-guide-to-java-versions-and-features that made me realize that the Java I knew well from years ago is many versions behind the current one. Today’s Java has a Swift-like Switch block, a form of Optionals, improved Interfaces, etc. But the basic structure of the language is still Java. I was looking for something beyond.
Then I found Kotlin: https://kotlinlang.org. On its home page it is advertised as “a modern programming language that makes developers happier”. That was enough for me to take a closer look.
I quickly realized that this language has some real potential for me. First, it is multi-platform since it compiles into Byte code that runs on Java’s virtual machine (JVM). Second, it is fast since it is essentially like Java when executed by the JVM. Third, and perhaps best of all, it looks like Swift! https://www.raywenderlich.com/6754-a-comparison-of-swift-and-kotlin-languages.
Below I have translated the spline code to Kotlin. If you compare it to the Swift version you will find that they are nearly identical with a few minor differences and only one major: Kotlin does not have Tuples. Instead, like Java, you have to create a data class to hold the variables you want to return to the caller. It is not as elegant as a Tuple but not a show stopper either.
At this time I don’t have more than a few weeks experience with Kotlin but they have been productive. By using Jetbrain’s IntelliJ IDEA https://www.jetbrains.com/help/idea/get-started-with-kotlin.html I have been able to quickly translate Java code to Kotlin. Furthermore, I’ve discovered that if I paste in Swift code into the editor I have a very good base for manually turning it into Kotlin. In future posts I will share my use of Kotlin and also make some comparisons with Python.
data class TupleOfDoubleArrays(var b: DoubleArray, var c: DoubleArray, var d: DoubleArray)
public fun spline(x: DoubleArray, y: DoubleArray): TupleOfDoubleArrays? {
val n = x.size
if (n != y.size) {
return null
}
val b = DoubleArray(n)
val c = DoubleArray(n)
val d = DoubleArray(n)
when (n) {
1 -> {
return null
}
2 -> {
val b0 = (y[1] - y[0]) / (x[1] - x[0])
b[0] = b0
b[1] = b0
return TupleOfDoubleArrays(b, c, d)
}
else -> {
val nm1 = n - 1
d[0] = x[1] - x[0]
c[1] = (y[1] - y[0]) / d[0]
for (i in 1 until nm1) {
d[i] = x[i+1] - x[i]
b[i] = 2 * (d[i-1] + d[i])
c[i+1] = (y[i+1] - y[i]) / d[i]
c[i] = c[i+1] - c[i]
}
b[0] = -d[0]
b[nm1] = -d[nm1-1]
c[0] = 0.0
c[nm1] = 0.0
if (n > 3) {
c[0] = c[2] / (x[3] - x[1]) - c[1] / (x[2] - x[0])
c[nm1] = c[nm1-1] / (x[nm1] - x[nm1-2]) - c[nm1-2] / (x[nm1-1] - x[nm1-3])
c[0] = c[0] * d[0] * d[0] / (x[3] - x[0])
c[nm1] = -c[nm1] * d[nm1-1] * d[nm1-1] / (x[nm1] - x[nm1-3])
}
//forward elimination
for (i in 1 until n) {
val t = d[i-1] / b[i-1]
b[i] -= t * d[i-1]
c[i] -= t * c[i-1]
}
//back substitution
c[nm1] /= b[nm1]
for (ib in 0 until nm1) {
val i = n - 2 - ib
c[i] = (c[i] - d[i] * c[i+1]) / b[i]
}
//compute polynomial coefficients
b[nm1] = (y[nm1] - y[nm1-1]) / d[nm1-1] + d[nm1-1] * (c[nm1-1] + 2 * c[nm1])
for (i in 0 until nm1) {
b[i] = (y[i+1] - y[i]) / d[i] - d[i] * (c[i+1] + 2 * c[i])
d[i] = (c[i+1] - c[i]) / d[i]
c[i] = 3 * c[i]
}
c[nm1] = 3 * c[nm1]
d[nm1] = d[nm1-1]
return TupleOfDoubleArrays(b, c, d)
}
}
}
data class TupleOfSplineResults(var s: Double, var sx: Double, var sxx: Double)
public fun evaluateSplineAt(u: Double, x: DoubleArray, y: DoubleArray, cv:TupleOfDoubleArrays): TupleOfSplineResults {
var i = 0
var j = x.size
do {
val k = (i + j) / 2
if (u < x[k]) { j = k }
if (u >= x[k]) { i = k }
} while (j > i + 1)
val dx = u - x[i]
val s = y[i] + dx * (cv.b[i] + dx * (cv.c[i] + dx * cv.d[i]))
val sx = cv.b[i] + dx * (2 * cv.c[i] + dx * 3 * cv.d[i])
val sxx = 2 * cv.c[i] + dx * 6 * cv.d[i]
return TupleOfSplineResults(s, sx, sxx)
}
public fun derivativesAtGridpoints(xL: Double, xU: Double, n: Int, u: DoubleArray, percent: Double = 0.0): TupleOfDoubleArrays {
val uxx = u.copyOf()
val x = u.copyOf()
//compute unequal grid spacing with percent increase
var sum = 0.0
for (i in 0 until n-1) {
sum += Math.pow((1.0 + percent / 100.0), i.toDouble())
}
val dx = (xU - xL) / sum
for (i in x.indices) {
if (i == 0) {
x[i] = xL
} else {
x[i] = x[i-1] + dx * Math.pow((1.0 + percent / 100.0), (i-1).toDouble())
}
}
val coefTuple = spline(x, u) ?: throw IllegalArgumentException("Spline Function Returning Null")
val ux = coefTuple.b
for (i in 0 until n) {
uxx[i] = 2.0 * coefTuple.c[i]
}
return TupleOfDoubleArrays(ux, uxx, DoubleArray(n))
}