When I last taught our postgraduate statistical computing course a few years ago I gave an assignment on classes, methods, and group generics in R:
Your assignment is to write a class and methods to handle numbers with units in the MKS (\(m\), \(kg\),\(s\)) system and the alternative FFF (furlongs, firkins, fortnights) system (where 1 fur = 201.168m, 1 fir = 40.82333133kg, and 1 ftn = 1209600s). You may use either S3 or S4 objects.
I thought I’d put up the code as an example.
First, set up some boring utility constants and functions
.mksnames<-c("m","kg","s")
.fffnames<-c("fur","fir","ftn")
.fff2mks<-c(201.168, 40.82333133,1209600)
.mks2fff<-1/.fff2mks
isWholeNumber<-function(x, tolerance=1e-8){
max(abs(x-round(x)))<=tolerance
}
isDimensionless<-function(x) all(getUnits(x)==0)The isWholeNumber() function confused the students a bit: why isn’t it just is.integer()? Basically, is.integer() refers to how the number is stored, and isWholeNumber() refers to the value. The constant 69 gives TRUE for isWholeNumber() but not for is.integer() because it’s a numeric (real, double) value.
We obviously now need getUnits(), and presumably getValues() as well, and to define them we need to know something about the class representation
setClass("mks",representation(values="numeric",m="numeric",kg="numeric",s="numeric"))
setClass("fff",representation(values="numeric",fur="numeric",fir="numeric",ftn="numeric"))
setGeneric("getUnits",function(object) standardGeneric("getUnits")) ## [1] "getUnits"setMethod("getUnits",signature="mks",function(object) c(object@m,object@kg,object@s)) ## [1] "getUnits"setMethod("getUnits",signature="fff",function(object) c(object@fur,object@fir,object@ftn)) ## [1] "getUnits"setGeneric("getValues",function(object) standardGeneric("getValues"))## [1] "getValues"setMethod("getValues",signature="mks",function(object) object@values) ## [1] "getValues"setMethod("getValues",signature="fff",function(object) object@values) ## [1] "getValues"Now, a constructor would be nice, since calling new directly in user code leads to bugs, bad software karma, and perhaps even updog.
mks<-function(values,m=0,kg=0,s=0){
if (!is.numeric(values))
stop("values must be numeric")
if (length(m)!=1 || length(kg) != 1 || length(s)!=1)
stop("dimensions must be of length 1")
if (!isWholeNumber(c(m,kg,s)))
stop("dimensions must be whole numbers")
new("mks",values=values,m=round(m),kg=round(kg),s=round(s))
}
fff<-function(values,fur=0,fir=0,ftn=0){
if (!is.numeric(values))
stop("values must be numeric")
if (length(fur)!=1 || length(fir) != 1 || length(ftn)!=1)
stop("dimensions must be of length 1")
if (!isWholeNumber(c(fur,fir,ftn)))
stop("dimensions must be whole numbers")
new("fff",values=values,fur=round(fur),fir=round(fir),ftn=round(ftn))
}The next basic step is to be able to print the values out. One component of that is printing the units, where we want to use something like \(ms^{-2}\) for acceleration. We can start with m^1.kg^0.s^-2, then drop exponents of 1 and drop dimensions with exponent 0.
formatUnits<-function(unitvalues, unitnames){
units<-data.frame(values=unitvalues,names=unitnames,stringsAsFactors=FALSE)
units<-subset(units, unitvalues!=0)
needspower<-units$values!=1
units$names[needspower]<-paste(units$names[needspower],units$values[needspower],sep="^")
paste(units$names[order(units$values,decreasing=TRUE)],collapse=".")
}
setMethod("show", signature="mks",
function(object) {
units<-formatUnits(getUnits(object),.mksnames)
values<-format(getValues(object),digits=getOption("digits"))
print(paste(values,units))
})## [1] "show"setMethod("show", signature="fff",
function(object) {
units<-formatUnits(getUnits(object),.fffnames)
values<-format(getValues(object),digits=getOption("digits"))
print(paste(values,units))
})## [1] "show"If we did this for more unit systems, life would become more pleasant with some sort of template system. But for two it isn’t bad.
And we can now have data: the densities and surface gravities of the nine planets
planet.g <-mks(c(3.7,8.9,9.8,3.7,23.1,9.0,8.7,11.0,0.6), m=1, kg=0, s=-2)
planet.d <- mks(c(5427,5243,5515,3933,1326,687,1270,1638,2390), m=-3,kg=1,s=0)
planet.d## [1] "5427 kg.m^-3" "5243 kg.m^-3" "5515 kg.m^-3" "3933 kg.m^-3"
## [5] "1326 kg.m^-3" " 687 kg.m^-3" "1270 kg.m^-3" "1638 kg.m^-3"
## [9] "2390 kg.m^-3"planet.g## [1] " 3.7 m.s^-2" " 8.9 m.s^-2" " 9.8 m.s^-2" " 3.7 m.s^-2" "23.1 m.s^-2"
## [6] " 9.0 m.s^-2" " 8.7 m.s^-2" "11.0 m.s^-2" " 0.6 m.s^-2"Yes, ok, pedants. The eight planets and the most important dwarf planet.
Having two unit systems isn’t much use if you can’t convert
as.mks<-function(x){
if (is(x,"mks")){
x
} else if (is(x,"fff")){
units<-getUnits(x)
conversion<-prod(.fff2mks^units)
mks(getValues(x)*conversion, units[1],units[2],units[3])
} else if (is.numeric(x)){
mks(x,0,0,0)
} else stop("Huh?")
}
as.fff<-function(x){
if (is(x,"fff")){
x
} else if (is(x,"mks")){
units<-getUnits(x)
conversion<-prod(.mks2fff^units)
fff(getValues(x)*conversion, units[1],units[2],units[3])
} else if (is.numeric(x)){
fff(x,0,0,0)
} else stop("Huh?")
} You could also use setAs() to make these into methods for as()
setAs("mks","numeric",
function(from) {if (isDimensionless(from)) getValues(from) else stop("can't drop units")})
setAs("fff","numeric",
function(from) {if (isDimensionless(from)) getValues(from) else stop("can't drop units")})The data, again
as.fff(planet.g)## [1] " 26910785970 fur.ftn^-2" " 64731350036 fur.ftn^-2"
## [3] " 71277216893 fur.ftn^-2" " 26910785970 fur.ftn^-2"
## [5] "168010582677 fur.ftn^-2" " 65458668576 fur.ftn^-2"
## [7] " 63276712956 fur.ftn^-2" " 80005039370 fur.ftn^-2"
## [9] " 4363911238 fur.ftn^-2"as.fff(planet.d)## [1] "1082251196 fir.fur^-3" "1045557955 fir.fur^-3" "1099800137 fir.fur^-3"
## [4] " 784318031 fir.fur^-3" " 264430640 fir.fur^-3" " 137001395 fir.fur^-3"
## [7] " 253263132 fir.fur^-3" " 326649615 fir.fur^-3" " 476613296 fir.fur^-3"Now, maths. You can multiply or divide any units by any other units, though you might need to convert first. Multplication adds the exponents; division subtracts them. If the two operands use different unit systems, the first operand wins.
There are lots of cases, but we make some use of ANY to simplify
setMethod("*",signature=c("mks","ANY"),
function(e1,e2){
v1<-getValues(e1)
u1<-getUnits(e1)
v2<-getValues(as.mks(e2))
u2<-getUnits(as.mks(e2))
u12<-u1+u2
mks(v1*v2,u12[1],u12[2],u12[3])
}) ## [1] "*"setMethod("*",signature=c("fff","ANY"),
function(e1,e2){
v1<-getValues(e1)
u1<-getUnits(e1)
v2<-getValues(as.fff(e2))
u2<-getUnits(as.fff(e2))
u12<-u1+u2
fff(v1*v2,u12[1],u12[2],u12[3])
} ) ## [1] "*"setMethod("*", signature=c("numeric","mks"),
function(e1,e2){
v2<-getValues(e2)
u2<-getUnits(e2)
mks(e1*v2,u2[1],u2[2],u2[3])
} )## [1] "*"setMethod("*", signature=c("numeric","fff"),
function(e1,e2){
v2<-getValues(e2)
u2<-getUnits(e2)
fff(e1*v2,u2[1],u2[2],u2[3])
} )## [1] "*"setMethod("/",signature=c("mks","ANY"),
function(e1,e2){
v1<-getValues(e1)
u1<-getUnits(e1)
v2<-getValues(as.mks(e2))
u2<-getUnits(as.mks(e2))
u12<-u1-u2
mks(v1/v2,u12[1],u12[2],u12[3])
}) ## [1] "/"setMethod("/",signature=c("fff","ANY"),
function(e1,e2){
v1<-getValues(e1)
u1<-getUnits(e1)
v2<-getValues(as.fff(e2))
u2<-getUnits(as.fff(e2))
u12<-u1-u2
fff(v1/v2,u12[1],u12[2],u12[3])
} ) ## [1] "/"setMethod("/", signature=c("numeric","mks"),
function(e1,e2){
v2<-getValues(e2)
u2<-getUnits(e2)
mks(e1/v2,-u2[1],-u2[2],-u2[3])
} )## [1] "/"setMethod("/", signature=c("numeric","fff"),
function(e1,e2){
v2<-getValues(e2)
u2<-getUnits(e2)
fff(e1/v2,-u2[1],-u2[2],-u2[3])
} )## [1] "/"Some more data: weight of a gallon of water on each planet (the firkin unit is the mass of a firkin of water, which is nine gallons).
mass.gallon.water <- fff(1/9,fir=1)
planet.g*mass.gallon.water## [1] " 16.782925 m.kg.s^-2" " 40.369739 m.kg.s^-2" " 44.452072 m.kg.s^-2"
## [4] " 16.782925 m.kg.s^-2" "104.779884 m.kg.s^-2" " 40.823331 m.kg.s^-2"
## [7] " 39.462554 m.kg.s^-2" " 49.895183 m.kg.s^-2" " 2.721555 m.kg.s^-2"Newton<-mks(1,kg=1)*mks(1,m=1,s=-2)
planet.g*mass.gallon.water/Newton## [1] " 16.782925 " " 40.369739 " " 44.452072 " " 16.782925 " "104.779884 "
## [6] " 40.823331 " " 39.462554 " " 49.895183 " " 2.721555 "Addition and subtraction only make sense if the dimensions are the same. We use the Ops group generic, whichs defines methods for all of +,-,*,/, but the ones for * and / won’t be used because we have more specific methods defined
setMethod("Ops", c("mks","ANY"), function(e1,e2){
v1<-getValues(e1)
u1<-getUnits(e1)
v2<-getValues(as.mks(e2))
u2<-getUnits(as.mks(e2))
if(!all(u1==u2)) stop("incompatible units")
mks(callGeneric(v1,v2),u1[1],u1[2],u1[3])
} )## [1] "Ops"setMethod("Ops", c("fff","ANY"), function(e1,e2){
v1<-getValues(e1)
u1<-getUnits(e1)
v2<-getValues(as.fff(e2))
u2<-getUnits(as.fff(e2))
if(!all(u1==u2)) stop("incompatible units")
fff(callGeneric(v1,v2),u1[1],u1[2],u1[3])
} )## [1] "Ops"setMethod("Ops", c("numeric","mks"), function(e1,e2){
u2<-getUnits(e2)
if(!all(u2==0)) stop("incompatible units")
mks(callGeneric(e1,getValues(e2)),0,0,0)
} )## [1] "Ops"setMethod("Ops", c("numeric","fff"), function(e1,e2){
u2<-getUnits(e2)
if(!all(u2==0)) stop("incompatible units")
fff(callGeneric(e1,getValues(e2)),0,0,0)
} ) ## [1] "Ops"Comparison operations are the same: we need the same dimensions for both operands, and we can use the Compare group generic. Once we know the units are the same, we just keep the units of the first operand.
setMethod("Compare", c("mks","ANY"), function(e1,e2){
v1<-getValues(e1)
u1<-getUnits(e1)
v2<-getValues(as.mks(e2))
u2<-getUnits(as.mks(e2))
if(!all(u1==u2)) stop("incompatible units")
mks(callGeneric(v1,v2),u1[1],u1[2],u1[3])
} )## [1] "Compare"setMethod("Compare", c("fff","ANY"), function(e1,e2){
v1<-getValues(e1)
u1<-getUnits(e1)
v2<-getValues(as.fff(e2))
u2<-getUnits(as.fff(e2))
if(!all(u1==u2)) stop("incompatible units")
mks(callGeneric(v1,v2),u1[1],u1[2],u1[3])
} )## [1] "Compare"setMethod("Compare", c("numeric","mks"), function(e1,e2){
u2<-getUnits(e2)
if(!all(u2==0)) stop("incompatible units")
mks(callGeneric(e1,getValues(e2)),0,0,0)
} )## [1] "Compare"setMethod("Compare", c("numeric","fff"), function(e1,e2){
u2<-getUnits(e2)
if(!all(u2==0)) stop("incompatible units")
fff(callGeneric(e1,getValues(e2)),0,0,0)
} ) ## [1] "Compare"Powers are interesting: we want to allow integer powers, and to allow fractional powers if the resulting dimensions are integers
setMethod("^",c("mks","numeric"), function(e1,e2){
if(isWholeNumber(e2)){
values<-getValues(e1)
units<-getUnits(e1)
mks(values^e2,units[1]*e2,units[2]*e2,units[3]*e2)
} else stop("Only whole-number powers allowed")
}) ## [1] "^"setMethod("^",c("fff","numeric"), function(e1,e2){
if(isWholeNumber(e2)){
values<-getValues(e1)
units<-getUnits(e1)
fff(values^e2,units[1]*e2,units[2]*e2,units[3]*e2)
} else stop("Only whole-number powers allowed")
})## [1] "^"setGeneric("root",function(e1,e2) standardGeneric("root")) ## [1] "root"setMethod("root",c("mks","numeric"),
function(e1,e2){
if (isWholeNumber(e2)){
units<-getUnits(e1)/e2
if (!isWholeNumber(units))
stop(paste("fractional dimensions:",paste(units,collapse=",")))
mks(getValues(e1)^(1/e2),units[1],units[2],units[3])
} else stop("only whole-number roots supported")
})## [1] "root"setMethod("root",c("fff","numeric"),
function(e1,e2){
if (isWholeNumber(e2)){
units<-getUnits(e1)/e2
if (!isWholeNumber(units))
stop(paste("fractional dimensions:",paste(units,collapse=",")))
fff(getValues(e1)^(1/e2),units[1],units[2],units[3])
} else stop("only whole-number roots supported")
})## [1] "root"setMethod("sqrt", "mks", function(x) root(x,2)) ## [1] "sqrt"setMethod("sqrt", "fff", function(x) root(x,2)) ## [1] "sqrt"And other maths stuff makes no sense unless the number is actually dimensionless
setMethod("Math","mks", function(x) {
if (isDimensionless(x))
callGeneric(getValues(x))
else stop("not defined for numbers with units")
})## [1] "Math"setMethod("Math","fff", function(x) {
if (isDimensionless(x))
callGeneric(getValues(x))
else stop("not defined for numbers with units")
})## [1] "Math"setMethod("^",c("mks","mks"), function(e1,e2) stop("You cannot be serious"))## [1] "^"setMethod("^",c("fff","fff"), function(e1,e2) stop("You cannot be serious"))## [1] "^"setMethod("^",c("mks","fff"), function(e1,e2) stop("You cannot be serious"))## [1] "^"setMethod("^",c("fff","mks"), function(e1,e2) stop("You cannot be serious"))## [1] "^"And some more examples
speed.of.light<-fff(1.8026175e12,fur=1,fir=0,ftn=-1)
one.year<-fff(26,ftn=1)
distance.to.star <- 4.3*speed.of.light*one.year
as.mks(distance.to.star)## [1] "4.054192e+16 m"speed.limit<-mks(100*1000,m=1)/mks(60*60,s=1)
time.to.star<-distance.to.star/speed.limit
one.century<-100*one.year
time.to.star/one.century## [1] "464078.7 "New Zealand was having a drought that year, and there was discussion of water prices for farmers. I asked how much it would cost at city water prices to make up the deficit for a \(700m^2\) urban section.
mm<-mks(0.001,m=1)
litre<-mks(1/1000,m=3)
sectionsize<-mks(700,m=2)
deficitvolume<- 100*mm*sectionsize
waterprice<-1.343/(1000*litre)
deficitvolume*waterprice## [1] "94.01 "An example of addition: mass of a litre of water plus a gallon of water
water.density<-mks(1e3,kg=1,m=-3)
litre*water.density+mass.gallon.water## [1] "5.535926 kg"And cube roots make sense for volumes: the side length of a cube that would hold the rainfall deficit
root(deficitvolume,3)## [1] "4.121285 m"