Back to BCIMA

For the moment I find myself back in the BCIMA code doing what I hoped I would be able to avoid: minor fixes that will facilitate or even permit the ABCIMA code to be written and work.

The big one was a redo of the proposalDensity function and its arguments. For some strange reason, I wrote the origional proposalDensity function to accept arbitrary inputs and then check to see if the right variable names came in. In hindsight, this is really bizarre as there are only four arguments to that function and this is really easy to set up. So I redid that function (cutting it down from like ten lines plus comments to five in the process) and had to redo all calls to it. However, the new structure is much more logical. Rather than having to pass genMean and genSigma (the proposal mean and proposal sigma respectively–I know, some of my `legacy’ naming conventions are horrid) as optional arguments which get checked and prompted for, now they are simply explicit arguments to the IMA, BCIMA, and eventually ABCIMA software. It’s not like these are odd arguments I’m not liable to use again or manipulate from higher levels (those I’ve tried to pick the `optimal’ way to use in the design), these are really big, important arguments. It was time.

All that said, the BCIMA code is still really slow. The IMA portion is already slow, so having to do a couple of these sorts of things takes even longer. Unless I figure out how to do some of it in C or Fortran, though, I’m not sure it’s going to get faster in any meaningful way. At this point I can only predict ABCIMA will run even slower because of the multiple BCIMA chains involved.

On a final note, I still need to impliment a flag to signal failure when the backchain hits the max backstep limit and make it do a restart on the backchain with warning to user. This will keep the software from failing out when the backchain fails to converge in a `timely’ fashion, but means instead that if it is really impossible there is a potential infinite loop, so I’ll need to impliment a counter and not let it try more than, say five times.

BCIMA

Well, it should be no surprise that BCIMA works and is in the package. It was working before and I modified it. A couple notes:

1. It really is much easier to do business with the new objects vice the older method. It is much more logical too.
2. I whipped up a sample output image. It looks good, but it also looks like I can’t get it to upload to mediafire right now. As soon as I do I’ll post the link here.
3. Speaking of what the image revealed, a 250 point chain gave a very nice estimate with a pretty off-center source point.

Once I get the image uploaded, I’m going to be pretty happy.

EDIT: Here’s the image. The red circle is the real center, the green the mean of the chain.

As for next steps, starting to work on the ABCIMA software is a big goal. I also need to run the check package utility on my most recent version in order to see if I can figure out what if anything makes it unhappy so I can get to work fixing those issues for eventual CRAN upload. After ABCIMA it’ll be time to get cracking on a couple hundred runs of IMA, BCIMA, and ABCIMA complete with time and error statistics for the requisite comparative analysis of the different approaches to this little problem.

Weekend Progress

This weekend’s progress was decent. I have the revised infimum function working again and have learned the hard way that the best way to handle arguments in a do.call method is to make a list outside of the method and then pass it to it vice trying to construct the list of arguments inside the do.call. So this means the next step is to go through all the bcima code (.backstep and bcMetropolis) and redo every one of my do.call statements to reflect this method. Truth is, it’s probably better this way because it makes the code easier to read.

One thing I’m going to have to read more about and figure out how to use is environments. They seems like a handy way to encapsulate what’s going on and then kill it all off, but I’m not sure what it’s going to get me in the short run. As such, I’m going to keep on ignoring it for the moment until I find I need to use it (and I won’t I suspect) but it will be a good tool to get used to for the sake of building more complicated structures over time. We’ll see what we see.

Today’s Progress

Today I knocked out most of my list from yesterday and incorporated the bcima code into the package, built the package, and am going to start testing tomorrow. Most of the changes were cosmetic in many ways but in the grand scheme are also a (I think) more intelligent way to handle the myriad data an analysis like this produces.

There are still some lingering variables which make me nervous because they are conditional on having other things be or not be true. For instance, if we are using a uniform proposal, then there is no need to supply a mean nor sigma for the bivariate normal distribution we would have needed when using a gaussian proposal. But right now the way this is handled is to make the mean and sigma optional arguments and have the code check if the user supplied it when needed. Part of me doesn’t like trusting the user this much, but since the commands aren’t really interactive there’s not much chance of my doing it differently.

Once bcima is debugged and working properly I’m going to take a hiatus from writing code and write documentation. There’s 20 something functions I need to write or rewrite documentation for and there’s no better time to start than soon.

Starting to Rework the Backward Coupling Code

Yes, I need to in many ways redo the backward coupling code and that means redoing some of the transition code and some other odds and ends (maybe).

So I spent my time today making a list of things to do to the code. On the one hand, it’s going to mean a lot of work. On the other hand, much of what I’ve identified means passing the already defined objects and just cherry picking some variables to use and I’m done!

One example is the tsargs argument to bcMetropolis. It turns out all I need to actually pass is the spectArcPolygon object previously created. This is, yes, easier. Another is the transition function. Several of the variables can be wrapped up and replaced by passing spectIC, spectArcPolygon, and spectData objects.

At some point, and after this round of code revisions gets done it should happen, I need to get the ABCIMA code written. Then I need to write wrappers which will automate the multiple cases I’ll need to run to conduct the comparative analysis of the three methods. Plus there’s the wrappers to write which will take a spectIC object and a spectData object and run an entire analysis from there. But this is in the future–and by future I mean months out.

IMA has been Modified

The IMA software has been modified to use the new, slightly object oriented version of the other scripts. Better still, it’s been tested and works great.

The reason I keep going over and over and over the code for the stuff up to and including IMA is mainly because it’s the base of all the other work and really if all of that stuff works well it’s much easier to make the stuff down the road work. As such, putting in the legwork to make the older stuff easier to use and maintain pays off down the road when I can’t come back to it for six months and want to make forward progress (like now).

The next thing to do is look at the backward coupling code and do whatever it is I think it will need to make things work, which really means redoing the inputs to the IMA script contained in the BCIMA code. Once that gets done I can start looking at what it will take to make adaptation work.

Perhaps the most gratifying part about working on this code when it is as advanced as it has become is that I can make good forward progress without needing to much around with the underlying software much. The pieces are already, in many ways, there for the using.

On the other hand, the documentation is almost non-existant and without that there’s not much to be said for the software. In addition, I would like to get the kinks worked out so that the package can go up on CRAN (which means figuring out how to do the licensing . . . grrr) and writing a paper which uses it to compare and contrast performance of the three different algorithms. This is all plenty of work to be done before I can call this project anywhere near finished.

Truncated Samples Now with (Some) Object Orientation

I’ve added a little object orientation to the truncated samples bit. This particular function really doesn’t need anything as far as an output object goes as it just returns a list of points. On the other hand it is pretty handy to have some object orientation going in.

First, having a well-built arcPolygon object saves lots of headaches as it keeps a bunch of relevant data together in one place for the using.

Second, I’m finding in general that it helps to force the use of an object of a particular class because it means I can rely on certain datapoints being there and having certain kinds of attributes. This makes it easy for me to write code because I don’t have to worry about what crazy names people might give to variables since the object is a standard sort of a thing and less liable to be influenced by the whimsy of people.

Now I just need to move forward and make an required modifications to the IMA and BCIMA software and then after that make the BCIMA software work and beyond that the ABCIMA software. At some point writing some good documentation of all this is going to be a requirement, but for the moment it feels good to just be programming again.

Some Progress (Finally)

A quick progress report from my first day working on the code base in a couple months.

• I’ve made an initial conditions object and it works (very nicely too . . . no need to build the initial conditions from scratch anymore)
• There’s a new spectData object which contains the info on a raw SPECT image apart from the initial condition information
• There’s the new arcPolygonVertices object which has the info for the Arc Polygon which can be fed into such gems as the Truncated Distribution Sampler.

Right now I’m working on modifying the truncated distribution sampler to work with the new data object. Unlike the earlier bits where no massive rewrites were required, this one might take just that in order to make all the pieces work together nicely.

Progress Tonight

Took a couple hours to work through some edits to the source files.  The big one is that I broke the huge block of arguments which are normally fed into the data generator into two smaller bits, one of which is now the output of an initial conditions function and the rest are still in the data generator. The upshot is that lots of the stuff which we might reuse is in the initial conditions and most of it doesn’t get used again. The big exception is the heights of the planes/the radius of the projection. I might incorporate that information into the dataset what’s in the dataset is the core information which will get used over and over again by different things.

Some Reworking Of The Codebase

On a whim I looked up how to do some basic OOP in R. After doing that I thought about what it would take to add some to my software and realized while OOP would only get minimal play in what I am doing, the awkwardness of having to manipulate bizare sets of variables in order to make different functions work was getting old, and consequently I’ve started to rework how some of the variables are passed to functions by in a sense forcing the user to employ the same trick I had been using in my R environment: create named lists for each of the generic sets or arguments I needed to pass and cheat and use those.

It’s a rather unsatisfactory state of affairs I was in, having to worry about these lists and then needing to write lots of checks into the functions to ensure all the right arguments are there. Instead I will be moving a lot of that out into seperate functions or otherwise handling it.