summaryrefslogtreecommitdiff
path: root/CakePenalties.R
blob: 5506b7f8fa52f307f29aabbe7595f96c3ae7f45b (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
# $Id$
# Some of the CAKE R modules are based on mkin
# CAKE (6245, 7247, 8361, 7414), by Tessella, for Syngenta: Copyright (C) 2011-2016 Syngenta
#
#    The CAKE R modules are 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)
}

Contact - Imprint