Optimizing .*: Details of Vectorization and Metaprogramming


For many of us mathematicians, we were taught to use MATLAB, and we were taught to vectorize everything. I mean obviously if we have matrices $$A$$, $$B$$, and $$C$$ and want to multiply element-wise (say to solve a reaction-equation at each point in space), then the optimal code is

A.*B.*C

No questions to ask, right? Actually, this code isn’t as optimized as you’d think. Lets dig deeper.

BLAS, Linpack, and MKL

The reason you are always told by “the lords of numerical math” to vectorize your code is because very smart programmers worked really hard on making basic things work well. Most of the “standardized” vectorized computations are calling subroutines from packages known as BLAS and LINPACK. To see what version your MATLAB is using, you can call

READ MORE

What Exactly is Brownian Motion and Why Does it Matter?


When you talk about randomness, stochastics, and noise, everything always seems to be going back to “Brownian motion” and “white noise”. But what exactly are these things?

The problem is that if you ask someone who researches the subject, you most likely get many different definitions at the same time. The reason is because, although there are many equivalent definitions (and terms), every single one only seems to be a small part of the intuition of what this thing actually is. So I am going to state without proof equivalent formulations of $$B(t)$$, which is the same as $$B_{t}$$, which is the same as $$W_{t}$$… READ MORE

Quick Optimizations in Julia for Performance: A Practical Example


January 19 2016 in Julia, MATLAB | Tags: , , , , | Author: Christopher Rackauckas

Let’s take a program which plots the standard logistic map:

r = 2.9:.00005:4;
numAttract = 100;
steady = ones(length(r),1)*.25;
for i=1:300 ## Get to steady state
  steady = r.*steady.*(1-steady);
end
x = zeros(length(steady),numAttract);
x[:,1] = steady;
for i=2:numAttract ## Now grab some values at the attractor
  x[:,i] = r.*x[:,i-1].*(1-x[:,i-1]);
end
using PyPlot;
fig = figure(figsize=(20,10));
plot(collect(r),x,"b.",markersize=.06)
savefig("plot.png",dpi=300);

This plots the logistic map. If you take the same code and change the array … READ MORE

Julia iFEM1: Porting Mesh Generation


My first project on the quest for a Julia finite element method is a simple homework problem. Just some background, this is for UC Irvine’s graduate Computational PDEs 226B course where in the first quarter we did all sorts of finite difference methods and now is our first foray into finite element methods. The purpose of the project is to grasp the data structure enough to use simple tools (i.e. mesh creation and plotting) to create a finite element solver for Poisson’s equation in 2D and check the performance differences. For those who haven’t programmed much this is a great learning exercise, but being a pretty standard exercise the MATLAB code took no time to writeup and so I started thinking: how does Julia compare to MATLAB for solving simple PDEs with the finite element method?

Testing this would take … READ MORE

WordPress Julia Syntax Highlighting


January 19 2016 in Julia | Tags: , , | Author: Christopher Rackauckas

using Julia
print("Test for Julia")
A = [0 1 ; 0 1];

Syntax highlighting in Julia exists? Yes, you can make it work even though no plugins currently support it. The way I was able to get it to work was to use the WP-GeSHi-Highlight WordPress plugin which uses GeSHi to perform the syntax highlighting. Sure enough, an unofficial Julia syntax file for GeSHi is hosted on GitHub. I manually plopped this file into the GeSHi folder in my WordPress install, used the command lang=”julia” and you can see at the top of this page that it worked!

The New Project


January 19 2016 in Uncategorized | Tags: | Author: Christopher Rackauckas

Hello everyone! I have just finished setting up a WordPress theme which matches my personal website. Now I am ready to start blogging. This will soon be updated with all sorts of fun insights into stochastic dynamical systems, their numerical approximations, and what they illuminate about biological organisms. Yes, I will be making liberal use of $$\LaTeX$$ and code. Sounds like fun to me.

Why am I doing this? There are many reasons. Honestly, for the most part it’s personal. Terance Tao among many others repeatedly writes/says that the only way to learn math/science is to simply do it. I have quickly learned that “reading” isn’t “doing”. It’s easy to procrastinate by simply reading a whole bunch of lovely math and science papers, but you really aren’t doing much until you’re digging deep into it, whether you’re explaining what you read or … READ MORE