I’m happy to announce that the first working release of spgl1++ is now available. It is based on the Matlab implementation of Michael Friedlander and Ewout van den Berg and is meant to be a standalone library which can be included in your C++ projects.

It’s still “early” in development, as such it only solves the regular basis pursuit problem with real measurement matrix :

subject to

Features I’m currently working on, which shouldn’t be far away, include solving the weighted version with denoising and complex measurement matrix :

subject to

I’ll be tweaking it and testing various matrix vector libraries. It currently works with Armadillo, a C++ linear algebra library whose syntax (API) is deliberately similar to Matlab.

Since the goal is to accomodate general matrix vector libraries, you willhave to provide spgl1++ with certain template specializations, such as how it does matrix vector products. For example, with Armadillo:

I say “should”, since I’ve only really tested it with Armadillo and I’m sure I’ve accidentally added some dependencies in the code and will have to work on making it as general as possible. In the coming days/weeks/months I’ll be doing testing and adding documentation.

One of my goals when I started my website just after I finished grad school was to start sharing little snippets of code I had used or developed over the course of my undergraduate and graduate degree. One thing I did use, and still use frequently, is the Fourier transform. The Fourier transform is used in countlessly many algorithms in seismic imaging because of it’s relationship to sampling theory, ability to identify frequency content in signals, and provide solutions to differential equations.

One day I stumbled upon a neat series of lectures by my undergraduate NSERC supervisor, Dr. Michael Lamoureux, for a Seismic Imaging Summer School in 2006. Hidden in a footnote of this set of lecture notes contains one of my favorite quotes about a mathematical theorem, namely the divergence theorem:

And, as I like to explain to my kids, the Divergence Theorem simply says you can measure how much methane gas a cow produces by either measuring how much is produced in each cubic centimeter of the cow, or simply by measuring how much leaks out its mouth and rear end, and other places on the surface of the cow. – Michael Lamoureux

However, I digress! In these lecture notes Michael describes a variety of theoretical considerations when investigating partial differential equations, and in particular the famous wave equation with initial conditions:

In the third lecture of these notes, the Fourier transform is used to find a solution of the equation. This is also given as an example in Elias Stein’s famous book on harmonic analysis (though I have it written here as Page 395, I’m not sure which volume. It’s likely volume 1). In this case, the Fourier transform of the solution can be written as the sum of sines and cosines:

In the future, I might include the derivation though for now this will be sufficient. We simply apply the inverse Fourier transform to obtain a solution. In 2009 or so, I wrote a little Matlab code to demonstrate this solution. There are some scaling factors to match that of the Fourier transform implementation of Matlab. Luckily I’m a bit of a hoarder when it comes to old school work so here it is:

To be honest, there are a lot of poor programming practices going on here. For example when I wrote this, prior to my introduction into software development, I probably gave myself a nice little pat on the back for putting in those comments. Those who employ clean coding practices see where I’m going with this. If you take a look through the code, every time there is a comment it could be naturally refactored into a function with a descriptive name. Moreover, who am I trying to kid with variable names like “W”! No one would quickly realize that I’m referring to a grid of frequencies. Though it is academic code in nature, it’s still embarrassing!

Moving forward, we can then use this function to plot the solution for a couple different initial conditions:

Notice that we see some funky stuff happening at the boundaries. The mathematical solution given is valid on the entire plane, so it should just disappear off our screen. The fact that there is reflections at the boundary is an artifact of the discrete nature of computational mathematics, and these are typically handled using things like absorbing boundary conditions or perfectly matched layers. The easiest, but most computational inefficient way of handling it, would be to make an even larger grid and only record for the amount of time it takes for the reflection to reach the area you are recording.

I have been using a library called libunittest for unit testing an open source project. A couple months ago I did a feature request for it, namely for an assertion for relative approximation. It’s an incredibly easy library to use and comes with a lot of documentation. Here is one of the examples which shows the “easy” way where you don’t have to register the test class:

C++

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

// easy.cpp

#include <libunittest/all.hpp>

usingnamespaceunittest::assertions;

TEST(test_value_is_true)

{

ASSERT_TRUE(true);

}

structfixture{

intvalue;

fixture():value(42){}

};

TEST_FIXTURE(fixture,test_with_fixture)

{

ASSERT_EQUAL(1337,value);

}

After installing, link against it and compile, then execute:

I personally like a different style that libunittest supports which is a bit more writing of code but I believe it’s a nice way of organizing a test suite:

C++

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

// test_template.cpp

#include <libunittest/all.hpp>

usingnamespaceunittest::assertions;

structtest_template:unittest::testcase<>{

staticvoidrun()

{

UNITTEST_CLASS(test_template)

UNITTEST_RUN(test1)

UNITTEST_RUN(test2)

}

voidtest1()

{

ASSERT_EQUAL(true,true);

}

voidtest2()

{

ASSERT_EQUAL(false,false);

}

};

REGISTER(test_template)

For learning about the other types of assertions you can do, you may want to peruse the testing code:

I had a problem this weekend where the following snippet would compile within Eclipse, however the IDE would still complain with things like “Symbol ‘array’ could not be resolved”:

C++

1

2

3

4

5

6

7

8

9

10

11

12

13

14

#include <array>

#include <algorithm>

#include <iostream>

intmain()

{

constintarraysize=10;

std::array<int,arraysize>a;

for(inti=0;i<arraysize;i++){

a[i]=((i*32315^393)&15);

}

std::for_each(a.begin(),a.end(),[](intv){

std::cout<<v<<'\n';});

}

There were lots of people who were able to fix this problem by simply adding the appropriate paths, but since it compiles in my case this wasn’t the solution. Anyway, some genius on Stackoverflow figured it out: http://stackoverflow.com/questions/17131744/eclipse-cdt-indexer-does-not-know-c11-containers

Add the symbol with name “__cplusplus” (which in this case is apparently an override):

__cplusplus

with the following value:
201103L

Now use “Run C/C++ Code Analysis and the red underlining from the unresolved imports goes away:

There is still that red on the left hand side which is just highlighting of the scope which bothers me. This can be removed by right clicking on the red portion, going to preferences, text editors, and removing the “Show range indicator.” check mark:

Finally something nice looking:

Note that for this example I still needed to ensure I was using the C++11 Language Standard:

This caused me a bit of grief but the solution is somewhat anticlimactic for how long it took me! When I tried to mex compile some C++11 code with the “out of the box” mex compiler I ran into the following errors of the following type:

warning: rvalue references are a C++11 extension [-Wc++11-extensions]

warning: deleted function definitions are a C++11 extension [-Wc++11-extensions]

no type named ‘shared_ptr’ in namespace ‘std’

And all sorts of weird stuff even though it compiled fine with clang++ on the command line. Essentially it came down to make the following adjustments to the mexopts file:

and then place pumpkin.png and theme.css into userContent folder in the home directory.

Once the plugin is installed, you will have access to some extra configuration settings in “Configure System” under “Manage Jenkins” that you will need to update:

So I’ve decided to take the plunge and have begun making repositories for my projects (old and new). I can’t believe I used to code without version control, it’s insanity! During grad school I had a neat project from my course in sparse matrix computation on Google Pagerank that I’m sharing on Github now. In my quest for making cool project names, I’ve dubbed this project LRPR(Low Rank Page Rank) and it can be found here: https://github.com/bee-rock/LRPR. Though it contains only one type of low rank approximation at the moment, perhaps that will change in the future.

The project starts out with a short introduction to Pagerank, formulated as en eigenvalue problem:

.

Abstract:

We investigate what the effect of a low rank approximation for the transition matrix has on the power method and an inner-outer iteration for solving the Pagerank problem. The purpose of the low rank approximation is two fold: (1) to reduce memory requirements (2) to decrease computational time. We show that we see an improvement in storage requirements and a decrease in computational time if we discard the time it takes to perform the low rank approximation, however at the sacrifice of accuracy.

I used to have some code that would extract a sample view count from the site: http://circle.ubc.ca/handle/2429/46448. The UBC Thesis collection, namely cIRcle: http://circle.ubc.ca/ is a large collection of dissertations, thesis, etc that students need to submit in order to satisfy a requirement for their degree. On each submission, there is a collection of statistics that monitors how many people have viewed and downloaded your thesis and from which country this was done from. UBC has since changed the way they show the statistics so it stopped working! So I changed the code around a bit so I could still provide a neat little script to do some datamining from a webpage.

I’ll be using two extremely useful Python libraries. One is called Requests, which essentially makes communicating with the internet websites a breeze. The second one is called Beautiful Soup, which makes the data extraction from web pages extremely simple.

I was always a bit confused about when and where try/catch blocks should go. Well, as a rule of thumb:

You should include in your try block a sequence of statements such that, if one of those statements fails, the rest of the sequence should be abandoned. No more and no fewer statements than that.