Scripting

In a manner similar to using Matlab/Octave/R/Scilab, we may use SuanShu in an interactive or a scripting environment. This is most appropriate for research, data analysis and rapid prototyping. Often using a script language, the numerical problems may be expressed and solved in a reduced number of code lines as compared to using Java. Also, the users can get immediate feedback for each line execution, without having to first finish the whole program.

Here is a screenshot of using SuanShu in Groovy.

groovy_matrix_example

Tools

There are a number of ways to use SuanShu in a scripting environment. We recommend the following.

Advantages

There are a number of advantages to using SuanShu over other mathematics scripting languages such as Matlab/Octave/R/Scilab. In the following comparisons, we are illustrating the pros and cons using R, but the arguments probably carry over to the other mathematics scripting languages.

No Language Pollution

There is no need to learn and use a very specialized and proprietary language. In general, the proliferation (pollution) of languages causes confusion. There are very established programming languages that junior students and professional programmers already know. The most notable examples are Java, C# and C++. They are general enough to solve any computationally feasible problems. It makes no sense to invent a new language just for an application. The consequences of language pollution are:

  1. the users need to learn a new language just to use a new application, hence a barrier of entry; e.g., the R language to use R;
  2. the heterogeneity of languages introduces unnecessary challenges when combining them to build one integrated system; the more the worse;
  3. debugging a variety of languages is difficult;
  4. the hosting applications are almost never written in Matlab/Octave/R/Scilab; we would thus have to pass a large amount of data back and forth between the applications (e.g., in Java) and the data processing code (e.g., R); this is a severe performance issue;
  5. the specialized language has no use outside its environment; e.g., no one uses the R language outside the R computing environment.

Most programmers already know Java (or can convert from C# and C++ very easily). As SuanShu is written entirely in Java, there is no need to learn a new language to use it.

Better Code

The “real” programming languages (e.g., Java, Groovy, Scala) are better designed than the specialized mathematics languages. They are invented by scientists who really understand programming languages. The modern languages support (or shall I say, demand) object-oriented code, inheritance, interface, anonymous class, closure, functional programming, an advanced memory model and management, automatic garbage collection and etc. In contrast, the popular mathematics languages have few, if not none, of these features. For example, the R language is a by-and-large procedural imperative C-like language. (Its object-oriented feature is nothing more like a struct in C. It should not be confused with the real object-oriented support in e.g., Java.) It might be acceptable a few decades ago when computer science was at its dawn but certainly not after modern programming languages have come a long way.

In a nutshell, the users of modern programming languages can write better code, in terms of correctness, readability, easiness to debug (maintainability), modularity (object-oriented).

Let’s illustrate our opinion by numerically evaluating this integral:

$\int_0^1 \! x^2 \, \mathrm{d}x$

In Matlab/Octave/R/Scilab, we typically would write something like this.

n = 10000000
sum = 0
for (x in seq(0, 1, 1/n)) {
	sum = sum + x * x # (*)
}
value = sum / n

This would be a piece of perfectly acceptable code many decades ago when C was the most popular programming language. Judging from the modern software engineering perspective, we make one criticism: the mathematics concepts involved are not modular; they cannot be separated out. All this snippet does is to describe one particular procedure to integrate one particular function. There is no reusable code. For example,

  • if you want to integrate another function, you would need to copy and paste this snippet then modify (*);
  • the integrand cannot be manipulated (as a first class citizen); it cannot be evaluated without making a copy outside the integration for-loop; it cannot be passed to another operation, e.g., differentiation, or another way of doing integration, without copying and pasting the same code.

Here is a modern object-oriented way of doing integration using SuanShu in Groovy. Note that the integration algorithm is now independent of the integrand. The integrator can take different integrands, and the integrand is no longer tied to the integrator. This is not to mention that the code is a lot shorter. One statement is to define a Riemann integrator; another statement is to define the integrand. The two objects are independent.

I = new Riemann();//create an integrator
v = I.integrate(
	['evaluate' : {x -> return x * x}] as UnivariateRealFunction,//integrate the function y = x^2
	0, 1);//limit: [0, 1]

Software Engineering Process

Not only does the modern languages like Java allow the users to produce more elegant code, they also come with a suite of software engineering tools and paradigms to improve the coding experience as well as the quality of end products.

First, the computer science community has escaped from the procedural imperative programming mindset to object-oriented programming. We have designed many modern programming techniques, most notably design patterns for effective programming. However, it is not clear if the legacy mathematics languages support inheritance and interface to the extent that we can code patterns using them.

The modern languages come with tools to ensure code correctness. For example, for each mathematics concept, SuanShu has a corresponding class(es). The JUnit tool allows us to systematically and automatically do regression tests on all the source code regularly and when there are code changes. The same goes true for Groovy, Scala and etc. In contrast, suppose we have 100s of R scripts that depend on each other. When we make code changes in one of them, it is not clear how to systematically and automatically make sure that the changes are correct for this particular R file and that all other scripts will work after the changes. Specifically, suppose we make a fundamental change in the way the coefficients in an ARMA model are estimated, will all other functions that depend on this ARMA fitting still work? Theoretically maybe yes, but practically often no (unless you always write bug free code).

Below is a screenshot of the result of SuanShu unit testing.

junit_example

 

OPTIM

This is a collection of the advanced optimization algorithms.

  • linear programming (LP)
  • quadratic programming (QP)
  • second order conic programming (SOCP)
  • semi-definite programming (SDP)
  • sequential quadratic programming (SQP)
  • mixed linear integer programming (MIP)
  • genetic algorithms (GA)
  • differential evolution optimization (DE)

STATS

This is a collection of advanced statistics algorithms.

  • cointegration
  • covariance selection
  • factor analysis
  • hidden Markov model (Rabiner)
  • hidden Markov model (mixture models: Beta, Binomial, Exponential, Gamma, Log-Normal, Normal, Poisson)
  • hypothesis testing
  • Kalman filter
  • LASSO (Least Absolute Shrinkage and Selection Operator)
  • linear regression
  • principal component analysis
  • stochastic process simulation
  • time series analysis

Hypothesis Testing

  • Kolmogorov-Smirnov
  • D’Agostino
  • Jarque-Bera
  • Lilliefors
  • Shapiro-Wilk
  • One-way ANOVA
  • T
  • Kruskal-Wallis
  • Siegel-Tukey
  • Van der Waerden
  • Wilcoxon rank sum
  • Wilcoxon signed rank
  • Breusch-Pagan
  • Glejser
  • Harvey-Godfrey
  • White
  • ADF
  • Bartlett
  • Brown-Forsythe
  • F
  • Levene

Linear Regression

  • OLS regression
  • logistic regression
  • GLM regression
    • family: Binomial, Gamma, Gaussian, Inverse Gaussian, Poisson
    • quasi-family: Binomial, Gamma, Gaussian, Inverse Gaussian, Poisson
  • GLM model selection

Time Series Analysis

  • various moving average filters
  • univariate/multivariate time series analysis
    • ARIMA modeling, simulation, fitting, prediction
    • GARCH modeling, simulation, fitting, prediction
    • Auto-Covariance computation

Tutorials

The NM Dev Tutorials

The NM Dev tutorials are practical guides for those who would like to create mathematics applications using the NM Dev library. They include lessons with many complete and working examples. Groups of related lessons are organized into “trails”.

For the most up-to-date information on using the NM Dev library, please consult the Javadoc. Please check out the NM Forum, the place to discuss the NM Dev library, mathematics and programming. You are also encouraged to post your questions there to get help beyond documents.

Before you start these lessons, please make sure you have set up the development environment properly. You can find the details in the Setup Guide.

Basic Trails

Specialized Trails

 

HelloNMDev

Hello NM Dev

This tutorial walks you through a simple first example on using the NM Dev library. You will need to have a NM Dev project to follow the examples here. For details on how to create or download a NM Dev project, please read the Setup Guide.

We will create some matrices and vectors and do some simple linear algebra. Java is an object-oriented language. So, other than the primitive types, like int and double, all variables are objects, namely, instances of some object definitions called classes. For example, to have a matrix instance, a matrix variable, or simply a matrix, we create an instance of the Matrix class.

[pastacode lang=”java” manual=”Matrix%20A1%20%3D%20new%20DenseMatrix(%0A%20%20new%20double%5B%5D%5B%5D%7B%20%2F%2F%20create%20a%20matrix%0A%20%20%20%20%7B1%2C%202%2C%203%7D%2C%0A%20%20%20%20%7B4%2C%205%2C%206%7D%2C%0A%20%20%20%20%7B7%2C%208%2C%209%7D%0A%7D)%3B” message=”Java matrix” highlight=”” provider=”manual”/]

We can print out the Java matrix to check.

[pastacode lang=”java” manual=”System.out.println(A1)%3B” message=”” highlight=”” provider=”manual”/]

In NetBeans the short cut for System.out.println is to type “sout” (without the “”) and then hit the “Tab” key.

To fix the imports, hit “Ctrl+Shift+I”. There are multiple choices. Choose com.numericalmethod.suanshu.matrix.doubles.Matrix.

fix_all_imports_matrix

Hit “F6″ to run.

Vola! You will see this output.

[pastacode lang=”bash” manual=”%5B%2C1%5D%20%5B%2C2%5D%20%5B%2C3%5D%0A%5B1%2C%5D%201.000000%2C%202.000000%2C%203.000000%2C%0A%5B2%2C%5D%204.000000%2C%205.000000%2C%206.000000%2C%0A%5B3%2C%5D%207.000000%2C%208.000000%2C%209.000000%2C” message=”” highlight=”” provider=”manual”/]

Congratulations! You have just successfully created your first matrix using the NM Dev library. We are now ready to do more.

Let’s check whether this matrix has a full rank. We do

[pastacode lang=”java” manual=”int%20rank1%20%3D%20MatrixMeasure.rank(A1)%3B%20%2F%2F%20compute%20A1’s%20rank%0ASystem.out.println(rank1)%3B%20%2F%2F%20rank1%20is%202%20implies%20a%20deficient%20matrix” message=”” highlight=”” provider=”manual”/]

The MatrixMeasure.rank method computes the numerical rank. For more information, see Compute the numerical rank of a matrix.

Hit “F6” to run.

We see that the rank of A1 is 2 so A1 is a rank deficient matrix. We can make it a full rank matrix by changing one entry.

[pastacode lang=”java” manual=”A1.set(3%2C%203%2C%2010)%3B%20%2F%2F%20change%20A1%5B3%2C%203%5D%20to%2010″ message=”” highlight=”” provider=”manual”/]

Note that NM Dev counts from 1!

Let’s try again.

[pastacode lang=”java” manual=”rank1%20%3D%20MatrixMeasure.rank(A1)%3B%20%2F%2F%20compute%20modified%20A1’s%20rank%0ASystem.out.println(String.format(%22rank%20for%20the%20modified%20A1%20%3D%20%25d%22%2C%20rank1))%3B” message=”” highlight=”” provider=”manual”/]

OK. A1 is now a full rank matrix. We can now find its inverse.

[pastacode lang=”java” manual=”Matrix%20A2%20%3D%20new%20Inverse(A1)%3B%20%2F%2F%20compute%20A1%20inverse%0ASystem.out.println(A2)%3B” message=”” highlight=”” provider=”manual”/]

Hit “Ctrl+Shift+I” to fix the imports. Hit “F6″ to run. Here is the output.

[pastacode lang=”bash” manual=”3×3%0A%0A%5B%2C1%5D%20%5B%2C2%5D%20%5B%2C3%5D%0A%5B1%2C%5D%20%C2%A0-0.666667%2C%20-1.333333%2C%201.000000%2C%0A%5B2%2C%5D%20%C2%A0-0.666667%2C%203.666667%2C%20-2.000000%2C%0A%5B3%2C%5D%20%C2%A01.000000%2C%20-2.000000%2C%201.000000%2C” message=”” highlight=”” provider=”manual”/]

It looks good. To make sure, we multiply them together.

[pastacode lang=”java” manual=”Matrix%20A3%20%3D%20A1.multiply(A2)%3B%20%2F%2F%20A3%20%3D%20A1%20%25*%25%20A2%0A%2F**%20In%20numerical%20computing%2C%20we%20are%20usually%20satisfied%20to%20compare%20equality%20up%20to%20a%20precision.*%2F%0Afinal%20double%20precision%20%3D%201e-14%3B%20%2F%2F%20A3%20and%20the%203×3%20identity%20matrix%20are%20equal%20up%20to%2014%20decimal%20points%0Aif%20(AreMatrices.equal(A3%2C%20new%20DenseMatrix(3%2C%203).ONE()%2C%20precision))%20%7B%20%2F%2F%20compare%20A3%20to%20the%203×3%20identity%20matrix%0ASystem.out.println(%22A2%20is%20indeed%20A1’s%20inverse%22)%3B” message=”” highlight=”” provider=”manual”/]

Hit “Ctrl+Shift+I” to fix the imports. Hit “F6″ to run.

OK. It works.

Let’s try a vector. We do

[pastacode lang=”java” manual=”Vector%20y%20%3D%20new%20DenseVector(new%20double%5B%5D%7B%20%2F%2F%20create%20a%20vector%0A20%2C%2010%2C%204%0A%7D)%3B” message=”” highlight=”” provider=”manual”/]

Hit “Ctrl+Shift+I” to fix the imports. There are multiple choices. Choose “com.numericalmethod.suanshu.vector.doubles.Vector”.

fix_all_imports_vector2

We solve a linear equation of this form.

Ax = y

We need to set up a linear equation solver.

[pastacode lang=”java” manual=”LinearSystemSolver%20solver%20%3D%20new%20LinearSystemSolver(A1)%3B%20%2F%2F%20set%20up%20the%20linear%20equation%20solver%0AVector%20x%20%3D%20solver.solve(y)%3B%20%2F%2F%20solve%20A1%20%25*%25%20x%20%3D%20y%0ASystem.out.println(x)%3B” message=”” highlight=”” provider=”manual”/]

Hit “Ctrl+Shift+I” to fix the imports. There are multiple choices. Choose com.numericalmethod.suanshu.matrix.doubles.linearsystem.LinearSystemSolver.

x is a vector

[pastacode lang=”bash” manual=”%5B-22.666667%2C%2015.333333%2C%204.000000%5D” message=”” highlight=”” provider=”manual”/]

The source code can be browsed in our svn repository here.

NM Dev Setup Guide

Setup Guide

This tutorial is on how to set up the development environment for programming with the NM Dev library (the Java version). In order to use the NM Dev library, you need to know some Java programming. Sun (or Oracle) provides excellent online tutorials on learning Java. You can find a lot of easy-to-read-and-follow lessons in Sun’s Java Tutorials.

IDE

You will need an Integrated Development Environment (IDE) to program NM Dev. Although you may choose any IDE you like, we recommend Netbeans. You can download it from here.

 

To learn about how to use NetBeans, you will find a quick start tutorial here. The most important skill for you to learn is debugging. This is especially true for those who are converting from Matlab and R and are used to printing out values. NetBeans debugger allows you to easily inspect values, place conditional breakpoints, pause and resume your program and much more. A quick Google search gives these articles.

We cannot emphasize enough that it is imperative that you must master using a debugger to make your programming experience productive and enjoyable.

Hotkey is another tool that will increase your programming productivity. My favorite hotkey is Ctrl+Shift+F (or Alt+Shift+f if you are using windows). It makes my code look pretty and neat. A list of hotkeys is found here. You can customize them to your own taste by going to Netbeans -> Preferences -> Keymap (or Tools -> Options -> Keymap for Windows/Linux).

nmdev.jar

To program the NM Dev library, you will need a licensed copy of nmdev.jar (the actual file name depends on the version that you have). It is a compressed file that contains the NM Dev Java classes for numerical computing. To get a copy of the trial license, please email sales@nm.dev.

  nmdev-1.0.0.zip (13.1 MiB, 4,125 hits)

This section walks you through how to set up a NM Dev project. You may skip this section if you want to start programming NM Dev right away using the S2 IDE which comes in-built with NM Dev. Please read our Basic Trails to get started. If in the future you would like to update the NM Dev library for a newer version, you will need to come back to this section.

Create a NM Dev project

To create a NM Dev project in NetBeans, open the NetBeans IDE, click File -> New Project… (Ctrl+Shift+N). In ‘categories’, choose ‘Java with Ant’; in projects, choose ‘Java Application’. Click ‘Next’.

In ‘Project Name’, type whatever you like, e.g., HelloNMDev. In ‘Project Location’, type where you would like to save your project. Make sure the box ‘Use dedicated Folder for Storing Libraries’ is checked.Then click ‘Finish’.

We need to tell this Java project that we will be calling the NM Dev library. To do so, right click on the ‘Libraries’ Folder in the project tab and click ‘Add JAR/Folder’

Browse to where you saved nmdev.jar. Select nmdev.jar and then hit ‘Open’.

 

To install the NM Dev javadoc in NetBeans, we need to associate the nmdev.jar with its javadoc. Right click the nmdev.jar file you just added and select ‘Edit’.

In the popup window, under the ‘javadoc’ section, click ‘Browse’ and select the nmdev-javadoc.jar file. Then click ‘Open’.

Now, you have created an empty NM Dev project. Copy and paste these lines under “// TODO code application logic here“.

 

[box type=”shadow”]

System.out.println(“Hello NM Dev”);

 

Matrix A1 = new DenseMatrix(new double[][]{//create a matrix

{1, 2, 1},

{4, 5, 2},

{7, 8, 1}

});

System.out.println(A1);

 

Matrix B = new Inverse(A1);//compute the inverse of A1

Matrix I = A1.multiply(B);//this should be the identity matrix

System.out.println(String.format(“%s * %s = %s (the identity matrix)”,

A1,

B,

I));

[/box]

Hit Command+Shift+I (or Ctrl+Shift+I on Windows/Linux) to fix the necessary imports. Make sure these statements appear at the top of your program. Else you can input them manually.

[box type=”shadow”]

import dev.nm.algebra.linear.matrix.doubles.Matrix;
import dev.nm.algebra.linear.matrix.doubles.matrixtype.dense.DenseMatrix;
import dev.nm.algebra.linear.matrix.doubles.operation.Inverse;

[/box]

Hit Alt+Shift+F to beautify your code. Your code should look something like this.

 

Before you can run your code you will have to place your license file by doing

[box type=”shadow”]

import dev.nm.misc.license.License;

License.setLicenseFile(new java.io.File(“your license path”));

[/box]

Right click on the project name, e.g., HelloNMDev, and click ‘Run’ (or simply press F6 if you make this project your main project). Voila! You should see this in your output window (Ctrl+4).

Congratulations!

Javadoc

NM Dev’s javadoc is the complete reference to the NM Dev library. You can read it in a browser, e.g., Internet Explorer. You can bring them up dynamically during programming. Each time you press Ctrl+SPACE after the ‘.’ of any object you create, you will see a list of available methods for this object, and the javadoc for the method selected.

To get more information about programming NM Dev, please consult our tutorials.

 

 

OPDE

This is a collection of numerical algorithms to solve ordinary and partial differential equation problems.

  • Ordinary Differential Equation (ODE) solvers for initial value problem (IVP):
    • Euler’s method
    • Runge Kutta
      • 1st order Runge Kutta
      • 2nd order Runge Kutta
      • 3rd order Runge Kutta
      • 4th order Runge Kutta
      • 5th order Runge Kutta
      • 6th order Runge Kutta
      • 7th order Runge Kutta
      • 8th order Runge Kutta
      • 10th order Runge Kutta
    • Runge-Kutta-Fehlberg (RKF45) (adaptive step-size control)
    • Adams-Bashforth-Moulton predictor-corrector multi-step method
      • 1st order
      • 2nd order
      • 3rd order
      • 4th order
      • 5th order
    • solvers based on Richardson extrapolation
      • Burlisch-Stoer extrapolation
      • semi-implicit extrapolation (suitable for stiff systems)
    • first order system of ODEs
    • conversion from high order ODE to first order ODE system
  • Partial Differential Equation (PDE) solvers
    • finite difference methods:
      • elliptic problem:
        • iterative central difference method (for Poisson’s equation)
      • 1D hyperbolic problem:
        • explicit central difference method (for 1D wave equation)
      • 2D hyperbolic problem:
        • explicit central difference method (for 2D wave equation)
      • 1D parabolic problem:
        • Crank-Nicolson method (for 1D heat or diffusion equation)
      • 2D parabolic problem:
        • alternating direction implicit (ADI) method (for 2D heat or diffusion equation)