Sunday, March 16, 2014

Align two vector, first try, not perfect results

See "align_curves.20140316.R" at github/project.H2O2.tolerance/Feb14,2011.H2O2.titration.bioscreen.star.strain



#20140316 align growth curves by minimizing sum of errors 

rm(list=ls())
debug = 9; 

#
pairwise_sum_of_errors = function( X1, X2, Start1,End1,Start2, End2 ) {
  if( (End1 - Start1) == (End2 - Start2) ) {
    return( sum( ( X1[Start1:End1]-X2[Start2:End2] )^2))
  } else {
    print("Error: X1 and X2 should have the same lengths")
    return( NA); 
  }
}

#testing 
x1 = 1:100;  x2 = 5:104
pairwise_sum_of_errors( x1, x2, 5, 100, 1,96)
sum((x1[5:100]-x2[1:96])^2)
pairwise_sum_of_errors( x1, x2, 5, 100, 1,100)

#align two vectors with equal lengths
#Let's assume start2 lags behind start1 
#start2 = start1 + delta # 2nd start lags behind 1st start by delta
pairwise_cost_of_shifted_vectors = function( delta, X1, X2,  low.threshold=0.01, debug=10 ){
  delta = floor( (delta+0.5) ) #delta must be an integer
  X1[X1<low.threshold]=0;
  X2[X2<low.threshold]=0;
  Start1 = 1;
  Start2 = Start1 + delta;  #left truncation at delta
  End2 = length(X2)
  End1 = length(X1) - delta; #right truncation at delta
  if(debug){ print(paste('Start1', Start1, "Start2", Start2, "End1", End1, "End2", End2))
  }
  if( (End1 - Start1) == (End2 - Start2) ) {
    return( sum( ( X1[Start1:End1]-X2[Start2:End2] )^2))
  } else {
    print("Error: X1 and X2 should have the same lengths")
    return( NA); 
  }
}

#testing 1
x1 = 5:104; x2 = 1:100;  
pairwise_sum_of_errors( x1, x2, 1,96, 5,100 )
pairwise_cost_of_shifted_vectors( 4, x1, x2)

delta = 4
res1 = optim(c(delta), fn=pairwise_cost_of_shifted_vectors, X1=x1, X2=x2, method="Brent", lower=c(1), upper=c(100))
#return 3 not 4? 

#testing 2
x1 = 50:149; x2 = 1:100;  
pairwise_sum_of_errors( x1, x2, 1,50, 50,99 )
pairwise_cost_of_shifted_vectors( 49, x1, x2)

delta = 50
res2 = optim(c(delta), fn=pairwise_cost_of_shifted_vectors, X1=x1, X2=x2, method="Brent", lower=c(1), upper=c(100))
#return 50 not 49?

No comments:

Post a Comment