1 files changed, 56 insertions, 0 deletions
diff --git a/CakePenalties.R b/CakePenalties.R
new file mode 100644
index 0000000..3b2a0df
--- /dev/null
+++ b/CakePenalties.R
@@ -0,0 +1,56 @@
+# $Id: CakePenalties.R 216 2011-07-05 14:35:03Z nelr $
+# The CAKE R modules are based on mkin
+# CAKE (6245), by Tessella, for Syngenta: Copyright (C) 2011  Syngenta
+#
+#    This program is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation, either version 3 of the License, or
+#    (at your option) any later version.
+# 
+#    This program is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+# 
+#    You should have received a copy of the GNU General Public License
+#    along with this program.  If not, see <http://www.gnu.org/licenses/>.”
+
+# Penalty function for fits that should not be accepted
+CakePenalties<-function(params, modelled, obs, penalty.scale = 10, ...){
+	# Because the model cost is roughly proportional to points,
+	# make the penalties scale that way too for consistent behaviour
+	# with many or few data points
+	sum(CakePenaltiesLong(params, modelled, obs, penalty.scale, ...)$value) * dim(obs)[[1]]
+}
+CakePenaltiesLong<-function(params, modelled, obs, penalty.scale = 10, ...){	
+	r<-data.frame(name=character(0), value=numeric(0))
+	
+	# Flow fractions > 1
+	#-------------------------------
+	fails<-CakePenaltiesFF(params);
+	if(dim(fails)[[1]] > 0){
+		#print("Penalty failures:")
+		#print(fails)
+		# Add a failure record for each flow fraction that missed
+		names(fails)<-c('name', 'value')
+		fails$value = 100*(fails$value - 1)
+		fails$name = lapply(fails$name, function(x) paste("FF_", x, sep=""))
+		r <- rbind(r, fails)
+		#totalPenalties <- totalPenalties + 100*sum(fails$x - 1)	# Penalty for failure
+	}	
+
+	r$value <- r$value * penalty.scale
+	return(r)
+}
+
+CakePenaltiesFF<-function(params){
+	ff.values<-params["f_"==substr(names(params), 0, 2)]	# Flow fractions
+	if(length(ff.values) > 0){
+		ffs<-data.frame(t(sapply(strsplit(names(ff.values), '_'), FUN=function(x){ x[c(2,4)] })))	# Split into source/dest
+		names(ffs)<-c('source', 'dest')
+		ffs['value'] <- ff.values	# Add values
+		sums<-aggregate(ffs$value, list(s=ffs$source), sum)	# Sum of flows from each source
+		fails<-sums[sums$x>1.00001,]	# All compartments > 1
+	} else { fails<-data.frame(s=character(0), x=numeric(0)) }
+	return(fails)
+}
+\ No newline at end of file