Non-Linear Curve Fitting is Nature Publication Worthy?

Someone needs to explain to me why this is a Nature publication ...

Granted it is Nature Protocols which started its publication run when I was in grad school and is likely less high-impact as full blown Nature.  Seriously tough - using the Solver Add-in in Excel to do non-linear curve fitting gets published?  If a younger version of myself knew about this, I'd have a Nature publication as a junior in undergrad.

It should be noted that such Nature worthy protocols are accomplished in the following 6 lines of R code (in this case for fitting a Michaelis-Menton curve):

f = formula(v ~ vmax*s/(km + s))
data = data.frame(
   s = c(1250, 625, 313, 156, 78, 39, 20, 10),
   v = c(148.68, 123.66, 96.84, 69.76, 45.54, 27.84, 16.38, 9.24)
)

fit = nls(f, data=data, start=list(vmax=1, km=1),
          algorithm='port', trace=F, lower=0)
summary(fit)

(note: I'm including the blank line in my 6-line count) On a similar thread, does that mean I can get a publication for this?

PID Control-R

On a whim, I thought it might be fun to try to implement a PID control algorithm ... in R.
Yes, I know I have a strange idea of fun.

Background

PID control is a heuristic method of automatically controlling processes as wide ranging as water levels in tanks to the direction of ships running against a current.  What's important is that you don't need to know much about a process to control it with a PID controller.  Hence, why it is so broadly used.  I studied the theory behind it in college and even built a LabVIEW UI for use in a student teaching lab.

For more detailed information, the obligatory Wikipedia link should be a good starting point.  In fact, that's where I got all the information I needed since my control theory text is currently buried in a closet somewhere.

A PID controller produces an output (aka a "Manipulated Variable", MV(t) or u(t)) - usually from 0 to 100% - based on a weighted sum of three terms:

• Error between where you want the process to be (a setpoint) and where the process is
• How much error has accumulated over time
• How rapidly the error is changing

which naturally leads to the P, I, and D in the name and the three weighting parameters (gains) that need to be tuned:

• Proportional
• Integral
• Derivative

The math for this looks like:

Don't be scared, it isn't as gnarly as it looks.

For the moment, let's ignore the Kp on the outside while we sort out the individual parts being added together.

The main variable in the equation above is e(t).  This is the error between the process variable PV(t) - the thing you want to control - and the setpoint SP(t) - the value you want the process to be at.  So this makes e(t):

e(t) = PV(t) - SP(t)

This happens to also be the first term - the Proportional response.  This should make sense - move the control by about as much error that you currently see.

The second term - the Integral response - isn't as ugly as it looks.  Integrals can be thought of as the area under a curve - in this case e(t).  Numerically, this simplifies to the sum of e(t)*dt, where dt is the amount of time between readings of the process variable.  To complete the Integral response, this sum is divided by Ti, the integral time, which is an estimate of how long it will take to abolish all the accumulated error.  In terms of the overall sum, this adds more control output to account for error that was previously "missed".

Last is the third and final term - the Derivative response.  Here the derivative is of e(t) with respect to time - or how fast is e(t) changing at time t.  To compute this, take the difference between the current value of e(t) and its last value e(t-dt) and divide it by dt.  To complete the Derivative response, multiply the derivative value by Td, the derivative time, which is an estimate of how far in the future the error can be reliably predicted.  This term adds control output to account for error that might happen.

All of the above is multiplied by Kp, the proportional gain, which scales the whole control response.

Code

Process

First things first, we need a process variable.


The code above is a simple function that simulates a process that, uncontrolled, will grow exponentially over time.  For added realism, I've included some uniform noise to the process growth.  The value of the control variable (in this case `u`) decreases the process value.

Controller

To code the controller for the above process, we first need to initialize a few variables:

In the interest of book-keeping, I'm storing the values of the Integral and Derivative responses in vectors EI and ED, respectively.  This isn't necessary, and indeed if this were a continuously operating controller, I wouldn't want to for the sake of conserving memory space.

Next a setpoint.  Below is code to create a setpoint profile that has a few step changes in time.

From here on, it's a matter of running a simulation loop.


Folks savvy in calculus may note that my implementation of the integral isn't entirely textbook - meaning, I'm not using the Trapezoid rule or Simpson's rule to compute a more accurate area for each time step.  This is on purpose, making the algorithm easier to implement.  It also assumes that my time step, dt, is small enough that this doesn't matter much.

Running the simulation and plotting the necessary variables, shows that the controller works reasonably well. Note, I disabled the noise in the process variable for this first run.

Adding noise to the process doesn't affect the results much, just a noisier output in U (the manipulated variable):

Using a more complicated setpoint profile:

produces:

From here, one can have some fun (or do some work) playing with the parameters, process variable, and setpoint profile.

Summary

A PID controller is based on a reasonably simple algorithm that can be quickly implemented in R, or any other programming language for that matter.

Set operations on more than two sets in R

Problem

Set operations are a common place thing to do in R, and the enabling functions in the base stats package are:

• intersect(x, y)
• union(x, y)
• setdiff(x, y)
• setequal(x, y)

That said, you'll note that each ONLY takes two arguments - i.e. set X and set Y - which was a limitation I needed to overcome recently.

Googling around I found that you can apply set operations to multiple (>2) sets using the Reduce() function.  For example, an intersection of sets stored in the list x would be achieved with:
Reduce(intersect, x)

However, things get trickier if you want to do a setdiff() in this manner - that is find the elements that are only in set a, but not in sets b, c, d, ..., etc.  Since I'm not a master at set theory, I decided to write my own, brute force method to get intersections, unions, and setdiffs for an arbitrary number of sets.

Implementation

The function I wrote uses a truth table to do this where:

• rows are the union of all elements in all sets
• columns are sets

So each "cell" of the table is TRUE if the element represented by the row is in the set represented by the column.

To find an element that intersects all sets, the values across the row need to be all TRUE.
To find an element that is only in one set, only one row value is TRUE.  To determine which set that element is in, numeric values are applied to each column such that:

• col1 = 1
• col2 = 2
• col3 = 4
• col4 = 8
• ... and so on

The values above are multiplied by each row's TRUE/FALSE (ergo: 1/0) values and summed across to produce a row's inclusion value.

For example, if an element is only in set 2 (that is column 2) the corresponding row inclusion value is 2.  If an element is in both sets 2 and 3 (columns 2 and 3) the corresponding row inclusion value is 6.

Visually ...

	1	2	4
set:	x	y	z	inc.val
a	T	F	F	1
b	F	T	T	6
c	T	F	T	5
d	T	T	T	7

Thus, determining intersections or setdiffs becomes a simple matter of filtering by row sums.

Code

A copper toned publication!

At long last (1.5yrs since the first submission attempt to be exact), the research I worked on as a post-doctoral fellow has been published!

Click on the image above to head over to the article for some light reading.  A lot of work went into this paper, and there are many to thank who don't appear in the author list - basically everyone in my former lab. On a fun note, the model we developed and just about all the figures were produced in R.  Quite the conversion, since I started this work as a big proponent for Matlab.

Overly Honest Methods

A friend of mine shared this on Facebook - A curated list of tweets tagged #overlyhonestmethods.  Basically, how science is really done behind that sparkly-smarty-pants curtain.

http://storify.com/BeckiePort/overlyhonestmethods

Really filter by date in Picasa 3

I've been managing my photos in Google's Picasa since 2009.  I personally think it's a good application with all the social connectivity I need (Facebook, Google+).  I expected that Google, being a company built on information and search, would have made the search bar in Picasa as robust as their web search bar.  Unfortunately, not so.

If you use any of Google's web products, you should be familiar with searching by dates using keys like...

• date:
• between:
• before:
• after:

Well, these don't work in Picasa.  Worse, all that there is for filtering by date is a slider bar that limits the range from some time in the past to the current date. This is not at all useful if you are only looking for photos taken in the month of June, 2012.

Yes, Picasa naturally organizes photos in date named folders.  However, this convenience is lost when importing instant uploads from Google+, which just drops a pile of photos/videos into one "Instant Uploads" folder.

What the Picasa search bar does do is index all text fields/tags in your photos.  If EXIF data is saved with your photos, there is a field called "Camera Date" which holds the date (according to the camera) that the photo was taken.  This field is formatted as:

yyyy:mm:dd hh:mm:ss
(where hours are in 24-hr format)

So, if you want to search for June 2012 photos, type:

2012:06

into the search bar.  So, while this works to grab photos from a specific month, there still isn't a way to find photos from a specific date range.  Hopefully, this will be a new feature in an upcoming release.

My Job Just Got Shiny-er

Today I went to the RStudio site to post a feature request.  I noticed up in the top navigation a link called "Shiny".  I like things that are shiny.  It turns out Shiny is a new package the RStudio team has developed for easily making reactive web applications from R code.

One of the tools I built at my current job is a web application that uses R as the server side analysis engine with A LOT of HMTL/CSS and AJAX tomfoolery in the browser to make things responsive a pretty.

From the looks of it, the shiny new Shiny package will significantly reduce the UI development time.  In fact, it appears to make web UI development more like UI development in desktop python apps using [insert your favorite toolkit here] (note: so far I mostly have experience with wx).

Ever since I started using R (back in 2007) I often lamented the lack of a way to create easy to use GUIs that could encapsulate complex analysis scripts for non-R users.  In MATLAB, this is done using GUIDE or low level uicontrol() functions.  In 2009, RGG looked promising, and I also found SciViews-R, but I never really had/invested the time to look deeper into either of these tools.  Besides, doing things from the comfort of a web browser is all the rage these days (and certainly makes deploying app updates easier).

I'm looking forward to testing this out!