CHNOSZ: R/objective.R

  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
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
# CHNOSZ/objective.R
# objective functions for revisit(), findit()
# all objective functions in one file, added attributes  20121009 jmd 

# input: a1 or loga1 is matrix with a column for each species and a row for each condition
# output: vector with length equal to number of rows of a1 or loga1
# attributes:
# optimum - the target for optimization (minimal or maximal)
## absolute - is the absolute value of the function used in the calculation? (TRUE or FALSE)

SD <- structure(
  function(a1) {
    # calculate standard deviation
    SD <- apply(a1, 1, sd)
    return(SD)
  },
  optimum="minimal"
)

CV <- structure(
  function(a1) {
    # get the coefficient of variation
    SD <- SD(a1)
    CV <- SD / rowMeans(a1)
    return(CV)
  },
  optimum="minimal"
)

shannon <- structure(
  function(a1) {
    # shannon entropy
    p <- a1/rowSums(a1)
    H <- rowSums(-p*log(p))
    return(H)
  },
  optimum="maximal"
)

DGmix <- structure(
  function(loga1) {
    # Gibbs energy/2.303RT of ideal mixing
    a1 <- 10^loga1
    DGmix <- rowSums(a1*loga1)
    return(DGmix)
  },
  optimum="minimal"
)

qqr <- structure(
  function(loga1) {
    # correlation coefficient of a q-q plot (for log-normality comparisons) 20100901
    # first, a function to calculate qqr
    qqrfun <- function(x) {
      # this is to catch errors from qqnorm
      qqn <- try(qqnorm(x, plot.it=FALSE), silent=TRUE)
      if(class(qqn)=="try-error") qqr <- NA
      else qqr <- cor(qqn[[1]], qqn[[2]])
    }
    # apply the function to the rows of loga1
    qqr <- apply(loga1, 1, qqrfun)
    return(qqr)
  },
  optimum="maximal"
)

logact <- structure(
  function(loga1, loga2) {
    # the logarithm of activity of the species identified in loga2 (species number)
    return(loga1[, loga2[1]])
  },
  optimum="maximal"
)

spearman <- structure(
  function(loga1, loga2) {
    # Spearman's rho (rank correlation coefficient) 20100912
    spearfun <- function(a, b) {
      # based on help(dSpearman) in package SuppDists
      if(length(a)!=length(b)) stop("a and b must be same length")
      if(any(is.na(a)) | any(is.na(b))) return(NA)
      ra <- rank(a)
      rb <- rank(b)
      dr <- rb - ra
      d <- sum(dr^2)
      r <- length(a)
      x <- 1-6*d/(r*(r^2-1))
      return(x)
    }
    rho <- apply(loga1, 1, spearfun, b=loga2)
    return(rho)
  },
  optimum="maximal"
)

pearson <- structure(
  function(loga1, loga2) {
    # pearson correlation coefficient
    H <- apply(loga1, 1, cor, y=loga2)
    return(H)
  },
  optimum="maximal"
)

RMSD <- structure(
  function(loga1, loga2) {
    rmsd <- function(a, b) {
      # to calculate the root mean square deviation
      d <- b - a
      sd <- d^2
      msd <- mean(sd)
      rmsd <- sqrt(msd)
      return(rmsd)
    }
    RMSD <- apply(loga1, 1, rmsd, b=loga2)
    return(RMSD)
  },
  optimum="minimal"
)

CVRMSD <- structure(
  function(loga1, loga2) {
    # coefficient of variation of the root mean squared deviation
    RMSD <- RMSD(loga1, loga2)
    CVRMSD <- RMSD/rowMeans(loga1)
    return(CVRMSD)
  },
  optimum="minimal"
)

DDGmix <- structure(
  function(loga1, loga2) {
    # difference in Gibbs energy/2.303RT of ideal mixing
    a2 <- 10^loga2
    DGmix2 <- sum(a2*loga2)
    DGmix1 <- DGmix(loga1)
    DDGmix <- DGmix2 - DGmix1
    return(DDGmix)
  },
  optimum="minimal"
)

DGinf <- structure(
  function(a1, a2) {
    dginf <- function(a1, a2) {
      # informatic Gibbs energy/2.303RT difference between assemblages
      p1 <- a1/sum(a1)
      p2 <- a2/sum(a2)
      sum(p2 * log10(p2/p1))
    }
    DGinf <- apply(a1, 1, dginf, a2=a2)
    return(DGinf)
  },
  optimum="minimal"
)

DGtr <- structure(
  function(loga1, loga2, Astar) {
    dgtr <- function(loga1, loga2, Astar) {
      # calculate the Gibbs energy/2.303RT of transformation
      # from an initial to a final assemblage at constant T, P and 
      # chemical potentials of basis species 20120917
      # loga1 - logarithms of activity of species in the initial assemblage
      # loga2 - logarithms of activity of species in the final assemblage
      # Astar - starved (of activity of the species of interest) values of chemical affinity/2.303RT
      # remove the logarithms
      a1 <- 10^loga1
      a2 <- 10^loga2
      # the molal Gibbs energy in the initial and final states
      # (derived from integrating A = Astar - RTln(a) from a1 to a2)
      G1 <- -a1*Astar + a1*loga1 - a1/log(10)
      G2 <- -a2*Astar + a2*loga2 - a2/log(10)
      # calculate the change in molal Gibbs energy for each species
      DG <- G2 - G1
      # return the sum
      return(sum(DG))
    }
    # we need to index both loga1 and Astar
    DGtr <- unlist(lapply(seq(nrow(loga1)), function(i) {
      dgtr(loga1[i, ], loga2, Astar[i, ])
    }))
    return(DGtr)
  },
  optimum="minimal"
)

### unexported functions ###

# return the objective function named in objective,
#  or produce an error if the function has no optimum attribute
get.objfun <- function(objective) {
  # get the function with the name given in 'objective'
  objfun <- get(objective)
  # perform a check on its usability
  if(!"optimum" %in% names(attributes(objfun)))
    stop(paste(objective, "is not a function with an attribute named 'optimum'"))
  # return the function
  return(objfun)
}

Questions? Problems? Suggestions? or email at ian@mutexlabs.com.

Please suggest features or report bugs with the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.