Spline Interpolation using matrix math.

Started at 2 pm yesterday and finished up at at 1:50 am today.  Finished it after 3 hours today.  Total of 15 hours, with about 8 just implementing the matrix library.

https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/04_Spline.c
https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/matrix.c
https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/matrix.h

Compile with the following command:

gcc matrix.c 04_Spline.c

I did the matrix work earlier in support of this project.  That article is here: http://mystry-geek.blogspot.com/2014/02/a-small-c-matrix-library.html

The output and inline comments in the code step through the processing that it is doing.  I tend to keep the verbose output until I can get a good test harness in place so that I know I am not breaking things.

Output is:

time ./a.out 

h:     2.0000  2.0000  2.0000   1.0000   2.0000  

x: 1.0000   3.0000  5.0000  7.0000  8.0000  10.0000 
a: 6.0000   2.0000  5.0000 10.0000  4.0000   1.0000 
A = 
    1.000000000     0.000000000     0.000000000     0.000000000     0.000000000     0.000000000 
    2.000000000     8.000000000     2.000000000     0.000000000     0.000000000     0.000000000 
    0.000000000     2.000000000     8.000000000     2.000000000     0.000000000     0.000000000 
    0.000000000     0.000000000     2.000000000     6.000000000     1.000000000     0.000000000 
    0.000000000     0.000000000     0.000000000     1.000000000     6.000000000     2.000000000 
    0.000000000     0.000000000     0.000000000     0.000000000     0.000000000     1.000000000 

b_vector = 
    0.000000000 
   10.500000000 
    3.000000000 
  -25.500000000 
   13.500000000 
    0.000000000 

A^-1 = 
    1.000000000     0.000000000     0.000000000     0.000000000     0.000000000     0.000000000 
   -0.268343816     0.134171908    -0.036687631     0.012578616    -0.002096436     0.004192872 
    0.073375262    -0.036687631     0.146750524    -0.050314465     0.008385744    -0.016771488 
   -0.025157233     0.012578616    -0.050314465     0.188679245    -0.031446541     0.062893082 
    0.004192872    -0.002096436     0.008385744    -0.031446541     0.171907757    -0.343815514 
    0.000000000     0.000000000     0.000000000     0.000000000     0.000000000     1.000000000 

A^-1 * b_vector = c = 
    0.000000000 
    0.949685535 
    1.451257862 
   -5.254716981 
    3.125786164 
    0.000000000 

b = -2.6331  -0.7338   4.0681  -3.5388  -5.6677  
d =  0.1583   0.0836  -1.1177   2.7935  -0.5210  

Solution for x: -1.500000 with 6 points: nan 
First derivative at x: nan 

Solution for x: 1.500000 with 6 points: 4.703223 
First derivative at x: -2.514413 

Solution for x: 3.500000 with 6 points: 1.880994 
First derivative at x: 0.278629 

Solution for x: 4.500000 with 6 points: 3.318298 
First derivative at x: 2.679573 

Solution for x: 5.500000 with 6 points: 7.257174 
First derivative at x: 4.681145 

Solution for x: 6.500000 with 6 points: 10.595421 
First derivative at x: 0.877686 

Solution for x: 7.500000 with 6 points: 7.266116 
First derivative at x: -6.698375 

Solution for x: 8.500000 with 6 points: 1.882469 
First derivative at x: -2.932652 

Solution for x: 9.500000 with 6 points: 0.773192 
First derivative at x: 0.193134 

Solution for x: 15.000000 with 6 points: nan 
First derivative at x: nan 

real 0m0.009s
user 0m0.000s
sys 0m0.004s

As you can see from above the code does the following:

• Builds the (x, a) point list, with a= f(x).
• Builds the h list in points, this is the difference between 2 consecutive x values.
• Builds the A matrix from the h values.
• Builds the b vector from a formula using a and h
• Inverts the A matrix and multiplies it against the b vector to get the c_vector
• Copies the c_vector into points.
• Loads b into points using a formula involving h, j, and c.
• Loads d into points using a formula involving h and c.

Finished the following a few hours later, after sleeping:

• In the finished version we can release the matrices we generated and return the point value.
• We can reportedly generate interpolated values using the points table as is.
• Actually generate specific values.

To generate an f(x) for a given x, you have to:  

• figure out which interval between points you are in (j).
• if you are outside the bounds of the points, reject the request, return NaN.
• solve a third degree polynomial using the aj, bj, cj, dj, and xj and the provided x value.

Two ways to get the rate of change at any point

    Calculate a point on each side of the value you are looking at and find the slope of the line that connects the two points.

   You could easily program the computer to take the derivative of the simple formula for one of these lines, basically just pre-program a function that is the derivative.
  I went with calculating the derivative to get an absolutely exact estimated value. Yeah, exact estimate is a problem. :D 

---

We got a take home test with 21 points and the code I wrote calculated the result in .016 seconds right out of the box.  I was able to solve it before the class period was over.   It looks like a duck!  I also calculated the first derivative of the function so we can perform more analysis on the information. :D  Of course it took me about 20 hours to write the matrix code and this code to get to that point. :D

A small C matrix library.

I am taking both Linear Algebra and Numeric Analysis this quarter and the latest assignment in Numerics is crossing over into Linear territory.  We have to calculate Spline interpolation, which involves solving a system of equations equal to n+1 in   size, which is a linear matrix of n+1 by n+1 in size.  So first I needed a matrix library that I could use to handle that part of the problem.

The code is here:
https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/matrix.c
https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/matrix.h
https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/matrix_testharness.c


Most of it is very strait forward and only took about an hour or two to write. The one section that took the most time was the part that built the inverse.

5 hours worth of coding.

This bit was difficult because the first way I tried to solve the problem I had 4 layers of for loops and it was just too confusing to figure out. So I deleted all that code and  tried again.

I broke the problem down into elementary row operations and had two layers of for loops in each layer, which I could hold in my mind to solve this function.

This section may be worth putting into two
different functions so I can reuse the parts to solve other problems.

I also want to add in a check here to test the inverse against the original matrix to make sure it is works.  If it doesn't work, then I should free the used memory and return a NULL.

Tested the matrix inverse function against another system to make sure it worked

After the inverse was working I wrote the multiply function.  I tested mult first by multiplying a matrix against the identity matrix and corrected a few issues with that.  When you multiply by the identity you should get the original matrix.   Once that was working I multiplied the matrix by its inverse and should have gotten the identity matrix, but I did not.  It took me a bit to figure out what I had reversed and then it was working.  After that matrix addition, and scalar multiplication was easy.

I converted the one file into three files so that the matrix.c  can be used as a library.  If you look at the matrix_testharness.c file, you can see that main() is located there, and that out of all the files there is only ever a single main().  You can see that at the top of matrix_testharness.c file there is a #include "matrix.h" file which brings in all the function headers.   The command to compile this 3 file example is here:

gcc matrix.c matrix_testharness.c

If you look at the top of matrix.h you can see that the data type for matrix is mapped into void *.  This keeps the contents of all the matrix objects private so they can only be accessed in the matrix.c file. This is important because you want a single interface to each of your objects.  If you let developers peek into the data structures they will start writing into them all nimbly bimbly and soon half the logic for the object is fuzzed across 50 files, with code copy pasted almost randomly between them all.  A single interface keeps things clean and is testable.  If there is a bug that needs fixed, it can be fixed in a single place.

I went through and cleaned up every warning.  I used -Wall in gcc to get most of them, then cppcheck and splint to find the rest.

The determinant

I added in a recursive function to calculate the determinate.  First of all, only square matrices can be calculated.  If the matrix is 2x2 then it returns ad-cb.  If the matrix is larger it picks the top row, and iterates over the columns, multiplying the cell in the top row against a new matrix made by removing the row and column you are currently in.  The sign alternates each time, starting positive.  The cell from the top line you are on is multiplied by the sign and the determinate of the new smaller matrix.

There is an alternative method that doesn't use recursion documented here:  http://www.purplemath.com/modules/determs2.htm  This looks fairly strait forward. and can be done in just a few loops without recursion.  This might be important for a 1000 x 1000 matrix, that would have to have a stack of 1000 function calls the original way.

I implemented this second  way of calculating the determinate. Because it is not calling N! subroutines.  However there is a problem.  Apparently this only works for 3x3 matrices. So switching back to the other method for now.

So, a few days later...

I ran this library against a 21 point spline interpolation and the determinate calculation took 2/10ths of a second.  Without the calculation the code ran in .007 seconds.  So I realized that the fractions I was pulling out of each line were the determinate, accumulated it and saved it to a new det field I created in the matrix object. This cached value is NAN until it is set.  The inverse matrix det is set to the original matrices 1/det.  If you multiply by a scalar the det increases by x^r, where x is the scalar and r is the number of rows, but only for a square matrix.   I rewrote Mat_Determinate() to use a more light weight variety of the same algorithm I used to calculate the inverse, stopping when I am in upper triangular form, with the pivots all one and the row multiplier options used to perform this accumulated into the det variable.

To do:

Write an automated test harness to test all the functions.  Worked toward that by creating a .h file and the start of a test harness file.

I noticed that the way I am trying to find the inverse requires that the pivot cannot contain a zero.  So I need to add in a recursive search to swap rows around to make this happen. So I also need the elementary row operation of swapping rows.  This needs to look down at the rows under the current row.  Because we already checked the determinate, we know there is a solution. 

Check for memory leaks.

Hermite Interpolation.

My implementation is here:  https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/04_Hermite.c

For the Hermite interpolation to work you need to not just measure the x and f(x) for your plot points, but also your f'(x).  So if you are measuring time and distance, you should also grab the speed at the same points.  Often this extra data is overlooked, which is a shame because it can often be collected at very little additional cost.

Then you double the points, and use the f'(x) at depth 1 when otherwise using the two points would result in division by zero.

This procedure only gives you the constant prefixes for each polynomial in the solution.  This starts with f(x) at the zeroth level, then the first value calculated at each level after that.  I push them all into a stack and then build the polynomial as I return from the recursive calculation of the results.

The depth you reach is stored in the depth of pushed values.  The degree of polynomial you generate is given by depth minus one.  The polynomial is calculated as shown by the printout.

Official Makerbot filament feed guide replacement is now installed.

As previously stated, after a couple of weeks after buying my Replicator 2 I had problems with "air prints."  The prints would seem to start and then in the middle they would mysteriously just stop. This was caused by a defective design combined with natural variation in filament diameters.  The original design had a hard plastic nub with only a small o-ring to act as a spring, so that as it wore down and it hit a narrow spot on the filament, the printer would stop feeding plastic.  You could adjust the printer to work again, but this involved taking off a fan and fan housing from the side and messing around with a small hex set screw.  

I saw where many other people were having similar issues, so I downloaded the fix file that someone had put onto thingiverse, ordered parts from various sites, adjusted the printer one more time the old way, and then printed the new part in pla plastic.  

Over the summer Makerbot offered a replacement injection molded part with the springs and bearings and screws to finally fix the issue.  So I ordered it and waited about 5 months to receive the part.  I threw it on a shelf because the fix I had in place was working just fine.  

I had not printed in a while, and started to print something about 3 weeks ago.  I noticed that the filament would not feed and the replacement part was broken, after about a year of use.   So I installed the factory replacement part and began printing again. :D

So far the factory replacement part has been as good as the part that I printed from the internet.  Here is hoping that the fix will last the life of the printer. :D

Creating a customizer file for Thingiverse.

A friend of mine wanted me to duplicate a movement tray for his war gaming figures.  These trays keep a bunch of figs together so that you can move the tray, instead of having to move a bunch of individual figs.

So I created an openscad file that looked like this:

difference () {

    cube ( [107, 53, 3], center=true ) ;

    translate ([0,2,1]) {
        cube ( [107-2, 53-1, 3], center=true ) ;
    }
}

This created a 3D tray by creating a box shape and chopping out a slightly smaller block that is moved to one side and up a bit.

You work with vectors of [x, y, z] numbers because this is 3 dimensions or 3D.  There is a basic tutorial on openscad here:  http://peak5390.wordpress.com/2013/01/08/an-introduction-to-3d-modeling-with-openscad-openscad-basics/  There are many tutorials, but this one seemed to work for me.

And this is all fine and good, if all you wanted was that single movement tray.  But he wants to be able to build any number of different movement trays, for many different figures.  So I took the numbers and made them into variables and put the variables at the top of the file, and that allowed me to easily change the numbers to anything I wanted and the math is handled.

I had been seeing the customizer things on http://www.thingiverse.com and I thought that this would be a great first try for this techniques.  All I had to do was put in some markup in comments around the variables like the following:

// Width of the tray as it faces the enemy?
tray_width = 107; // [50:1000]
// Depth of the tray away from the enemy?
tray_depth = 53; // [25:500]
// Edge height of tray?
edge_height = 3;  // [2:5]
// How wide is the edge?
edge_width = 2; // [1:5]

difference () {

    cube ( [tray_width, tray_depth, edge_height], center=true ) ;

    translate ([0,edge_width,1]) {
        cube ( [tray_width-edge_width*2, tray_depth-edge_width+1, edge_height], center=true ) ;
    }
}

The comment just before each variable becomes the text in the customer for that control, and the range in the comment following the variable becomes the control type and range.  I just used sliders, but there are a variety of other controls you can use in your own customizer tags.  This is discussed here:  http://customizer.makerbot.com/docs  It  looks like you can even have controls in various named tabs, or create a hidden section so you can use variables that don't appear.

Then I uploaded the openscad file, added the customizer tag (this is important or you will not see the customer option for your thing), and a couple of other tags and was able to see the new object running in the customizer on thingiverse here:  http://www.thingiverse.com/thing:251785  I also created an stl file on my local system and uploaded that so that there was something for people to see.

There are some best practices to use here:  http://customizer.makerbot.com/best_practices
These tell you how to design for regular people, which is what these files are aimed at.  Anyone that is already advanced can already download your openscad file and configure the output anyway they like.

--

I took a second shot at this design, calling it the Easy Movement Tray.  http://www.thingiverse.com/thing:253049

This was after using the first tray designer to lay out a couple of trays.  I realized that for most people there was a much simpler way to do the design.

Instead of having to do all the measuring yourself, I allowed people to select a base size from a drop down list of sizes, then tell me how many figures wide and deep the tray should be.  I also put in a few little divots in the middle of the tray, between all the corners, for an easy way to see the sections where the bases should be.  This also helps to break up the large featureless flat plane.

I actually add one more mm to the base size, for base variation, so that 28mm becomes 29mm. After painting and flocking bases can be a little bigger than they started, so I wanted to compensate for this.   After that I calculate everything by the tray width and tray depth, just like the first tray.  I used a //[hidden]  tab for these variables, so that they wouldn't show up in the UI.

The main update I need to make to this file is to have a better drop down list of base sizes that reflect more sizes that people use in the real world.

La Grange and Neville's method implementation to fit an nth degree polynomial to a set of points

In numeric analysis class I had to create a solution for some curve fitting algorithms.  I choose to use C as well as implementing these functions in excel, because when I would like to automate data processing I would rather use a small program to perform the task as quickly as possible, rather than a manual process using Excel.

These functions can fit a curve to hundreds of points in hundredths of a second on a relatively slow net book.  The reason you would do this is that you may not know the function that governs a system, but you need to have a smooth continuous function that gives an accurate y value when you input an x value.  These curve fitting methods allow you to model an unknown function, inside the range of your data points, calculating the nth degree polynomial for your n points you know.  

The first is the Lagrange method, which builds a polynomial one term at a time, calling a second function to construct the parts of the polynomial for a specific x.

https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/04_Lagrange.c

The second one was Neville's method and I didn't understand how this was working until I manually did a small problem using excel to calculate the result.  I finally worked it all out and came up with this very simple solution:

https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/04_Neville.c

Unfortunately, as easy to understand as the above program is, it has a serious problem where it recalculates the sub results many times.  This results in taking many minutes for relatively small point sets.

So I re-implemented the solution in the below code to build the solution from the "base of the triagle" rather than from the apex, so that each part of the solution is only calculated once.  This allows things that took minutes to run in sub second time intervals.

https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/04_Neville-mod.c

The only tricky thing about this function is that I had to pass in the original set of points, unchanged, then use the depth that we are at in order to index the correct Point[i].x.  I tried just pushing a single value up, but the x values are selected diagonally so a single calculated value could need one of two returns for f(x).

The depth is really only an indication of the spread of the x points we need for the current calculation.

A benefit of this solution is that you can print out all the results of the sub functions and the data you used to generate the result, but uncommenting the printfs in this section.

This function also only calls itself n-1 times, and calculates just n-depth sub results in each call.
This value is stored in a net set of points and passed into the next lower depth of the function as it calls itself.

The function returns and unwinds from all the calls when it calls itself with just a single point.  All we do is return the value we last calculated and we are done. :D