How to make IPF code faster and more concise in R Iterative proportional fitting is a way of adjusting 'internal cells' in a multidimensional matrix to optimise fit. It is also known as 'raking' and can be seen as a subset of 'entropy maximisation'.
The purpose for which I use IPF is to allocated individuals to zones. My code iteratively adjusts a weight for each individual in each zone, depending on how representative the individual is of each zone. A practical introduction to the technique, with many links to other resources, can be found here.
Now, my problem is that the code seems longer and slower than it needs to be. I have thought this for a while, but am only asking the question now because I've recently discovered 4 alternative R implementations of IPF:


*

*This very recently published ipfp package implements IPF in C, so will be fast.

*There is a longer-standing version of IPF within the long-standing categorical data analysis package, cat.

*There is also a custom ipf() function that seems to be used by some, described on CRAN email list servers, dating back to 2008.

*An implementation for the Alaska Department of Labor and Workforce Development 
Options 2-4 are mentioned in a previous Cross Validated post that is very useful and which triggered this question. The top answer to this question proposed yet another method: glm().
Now, I can see that there are large potential benefits of using one or more of these options in my code to speed things up and to make things run faster. However, I can't see how to implement these new functions in the framework of my code. The problem, I suspect, is that I'm re-weighting individual-level data, whereas the other implementations are re-weighting aggregate data.
So, to put the question more concisely: which option is best suited to my application? Any guidance (preferably based on the reproducible example) on implementation greatly appreciated.
Reproducible example
This reproducible example can be found here. Cloning the repo https://github.com/Robinlovelace/IPF-performance-testing enables running additional iterations.
# Minimal IPF example
start_time <- Sys.time()
num.its <- 10 # iterations

# Read-in data (manually to start, will use scripts in future)
c.names <- c("id", "age", "sex")
ind <- c(       1, "50+", "m",
                2, "50+", "m", 
                3, "16-49", "m", 
                4, "50+", "f", 
                5, "16-49", "f")
ind <- matrix(ind, nrow = 5, byrow = T) # Convert long data into matrix, by row
ind <- data.frame(ind) # Convert this into a dataframe
names(ind) <- c.names # Add correct column names

# Read in the data in constraints
category.labels <- c("16-49", "50+",  "m", "f") # Age, and sex constraints
all.msim <- c(        8,        4,     6,    6  
) # each row represents an area, each column a category
all.msim <- matrix(all.msim, ncol = 4, byrow = T) # Convert long data into matrix, by row
all.msim <- data.frame(all.msim) # Convert this into a dataframe
names(all.msim) <- category.labels # Add labels
con1 <- all.msim[,1:2] ; con2 <- all.msim[,3:4]

ind.cat <- data.frame(cbind(model.matrix(~ind$sex - 1)[,c(2,1)], model.matrix(~ind$age - 1)))
names(ind.cat) <- category.labels

# create weights in 3D matrix (individuals, areas, iteration)
weights <- array(dim=c(nrow(ind),nrow(all.msim),num.cons+1)) 
weights[,,num.cons+1][] <- 1 # sets initial weights to 1
ini.ws <- weights[,,num.cons+1]

# convert survey data into aggregates to compare with census (3D matix)
ind.agg <- array(dim=c(nrow(all.msim),ncol(all.msim),num.cons+1))
for (i in 1:nrow(all.msim)){
  ind.agg[i,,1]   <- colSums(ind.cat) * weights[1,i,num.cons+1]}

# re-weighting for constraint 1 via IPF 
for (j in 1:nrow(all.msim)){
  for(i in 1:ncol(con1)){
 weights[which(ind.cat[,i] == 1),j,1] <- con1[j,i] /ind.agg[j,i,1]}}
for (i in 1:nrow(all.msim)){ # convert con1 weights back into aggregates
  ind.agg[i,,2]   <- colSums(ind.cat * weights[,i,num.cons+1] * weights[,i,1])}

# test results for first row
ind.agg[1,1:2,2] - all.msim[1,1:2] # should be zero

# second constraint
for (j in 1:nrow(all.msim)){
  for(i in 1:ncol(con2) + ncol(con1)){
  weights[which(ind.cat[,i] == 1),j,2] <- all.msim[j,i] /ind.agg[j,i,2]}}  
for (i in 1:nrow(all.msim)){ # convert con2 back into aggregate
ind.agg[i,,3]   <- colSums(ind.cat * weights[,i,num.cons+1] * weights[,i,1] * weights[,i,2])}

# for multiple iterations
wf <- array(dim=c(dim(weights), num.its, 1)) # array to store weights its, wei
indf <- array(dim=c(dim(ind.agg), num.its, 1))
wf[,,,1,1] <- weights 
indf[,,,1,1] <- ind.agg

# loop for multiple iterations (run e2.R repeatedly, saving each time)
for(it in 2:num.its){
source(file="e2.R")
wf[,,,it,1] <- weights
indf[,,,it,1] <- ind.agg
}

# Analysis - in general, see analyis files
a.v <- as.vector(as.matrix(all.msim)) # constraints in long form, for cor
g.v <- as.vector(as.matrix(indf[,,1,1,1]))
cor(a.v,g.v)
(total_time <- Sys.time() - start_time)

 A: You can also have a look at the mipfp package (http://cran.r-project.org/web/packages/mipfp/index.html), which implements  an iterative proportional procedure.
This package is able to update an initial N-dimensional table to fit to a set of multi-dimensional target margins (the dimension of the margins can vary). The margins can be incomplete (i.e. containing NA values) and not consistent. In the later case, the set of margins are scaled to probabilities so that they all sum to 1.
A: Found a way to make the fast ipfp() function from the ipfp package work. This is both fast and concise. Tests suggests it is 50 times faster than the original R code, and uses ~1/2 the number of lines. Reproducible example: https://github.com/Robinlovelace/IPF-performance-testing/blob/master/models/simple/etsim-ipfp.R
The only question that remains: how to view the source of ipfp()? But this certainly is not the forum for such a question. Will continue to follow advice from the excellent Advanced R wiki book for that: http://adv-r.had.co.nz/
A: You might get some value from taking a look at the mmap package in R. 
http://cran.r-project.org/web/packages/mmap/mmap.pdf
