Skip to content

Instantly share code, notes, and snippets.

@tslumley

tslumley/unit-examples.R

Created August 17, 2016 23:14
Show Gist options
  • Save tslumley/b2eaac120b45a41cfdc3e9a67aa1cb5f to your computer and use it in GitHub Desktop.
Save tslumley/b2eaac120b45a41cfdc3e9a67aa1cb5f to your computer and use it in GitHub Desktop.
Download ZIP
Arithmetic and logical operators for mks and furlong/firkin/fortnight units
Raw
unit-examples.R
gravity<-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)
density<-mks(c(5427,5243,5515,3933,1326,687,1270,1638,2390), m=-3,kg=1,s=0)
speed.of.light<-fff(1.8026175e12,fur=1,fir=0,ftn=-1)
setMethod("mean","mks",function(x,...){
u<-getUnits(x)
mks(mean(getValues(x)),u[1],u[2],u[3])
})
sqrt(mean(gravity^2)-mean(gravity)^2)
sqrt(mean(density^2)-mean(density)^2)
## how much does a litre of milk weigh on the other planets?
litre<-mks(1/1000,m=3)
densityofmilk<-mks(1030,m=-3,kg=1)
weightonplanets<-litre*densityofmilk*gravity
## simple tests
gravity*density
2*gravity
gravity/gravity
1-density/density
1+density
gravity+density
gravity-gravity
gravity*0
gravity^0
gravity/0
## mixed unit systems
gravity*2-as.fff(gravity)
## mass of litre of milk plus firkin of wate
litre*densityofmilk
mass.gallon.of.beer<-fff(1/9,fir=1)
mass.gallon.of.beer+litre*densityofmilk
sin(density)
log(speed.of.light)
atanh(gravity)
## ignoring relativity, how long to the speed of light at one g acceleration?
one.g<-mks(9.8,m=1,s=-2)
speed.of.light/one.g
## how long to get to Alpha Centauri at the NZ speed limit of 100kph.?
one.year<-fff(26,ftn=1)
distance.to.star <- 4.3*speed.of.light*one.year
as.mks(distance.to.star)
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
## NZ drought
## rainfall.deficit was 100 mm
## How much would this cost in Auckland tap water
## for a 700m^2 Auckland section
mm<-mks(0.001,m=1)
sectionsize<-mks(700,m=2)
deficitvolume<- 100*mm*sectionsize
waterprice<-1.343/(1000*litre)
deficitvolume*waterprice
## picky,picky
a<-mks(1,m=1+1e-14)
b<-mks(1,m=1)
a+b
2^one.year
mm^2
mm+one.year
mks(1,m=1.1)
sqrt(speed.of.light)
sqrt(speed.of.light^2)
root(deficitvolume,3)
root(deficitvolume,2)
## math operations can work if units are m^0kg^0s^0
sectionsize/(mm*mm)
log10(sectionsize/(mm*mm))
Raw
units4.R
## S4 solution
## this is easier, but more verbose
## multiple dispatch handles binary operators better
## internal/private. these wouldn't be exported from NAMESPACE in a package
.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
}
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=".")
}
isDimensionless<-function(x) all(getUnits(x)==0)
## end of private functions/variables
setClass("mks",representation(values="numeric",m="numeric",kg="numeric",s="numeric"))
setClass("fff",representation(values="numeric",fur="numeric",fir="numeric",ftn="numeric"))
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))
}
setGeneric("getUnits",function(object) standardGeneric("getUnits"))
setMethod("getUnits",signature="mks",function(object) c(object@m,object@kg,object@s))
setMethod("getUnits",signature="fff",function(object) c(object@fur,object@fir,object@ftn))
setGeneric("getValues",function(object) standardGeneric("getValues"))
setMethod("getValues",signature="mks",function(object) object@values)
setMethod("getValues",signature="fff",function(object) object@values)
setMethod("show", signature="mks",
function(object) {
units<-formatUnits(getUnits(object),.mksnames)
values<-format(getValues(object),digits=getOption("digits"))
print(paste(values,units))
})
setMethod("show", signature="fff",
function(object) {
units<-formatUnits(getUnits(object),.fffnames)
values<-format(getValues(object),digits=getOption("digits"))
print(paste(values,units))
})
## These could be generic, with methods for the other unit type and for numeric
## Or they could use setAs()
## This is simple but not terribly elegant
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 write out all eight cases explicitly
##
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])
})
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])
} )
setMethod("*", signature=c("numeric","mks"),
function(e1,e2){
v2<-getValues(e2)
u2<-getUnits(e2)
mks(e1*v2,u2[1],u2[2],u2[3])
} )
setMethod("*", signature=c("numeric","fff"),
function(e1,e2){
v2<-getValues(e2)
u2<-getUnits(e2)
fff(e1*v2,u2[1],u2[2],u2[3])
} )
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])
})
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])
} )
setMethod("/", signature=c("numeric","mks"),
function(e1,e2){
v2<-getValues(e2)
u2<-getUnits(e2)
mks(e1/v2,-u2[1],-u2[2],-u2[3])
} )
setMethod("/", signature=c("numeric","fff"),
function(e1,e2){
v2<-getValues(e2)
u2<-getUnits(e2)
fff(e1/v2,-u2[1],-u2[2],-u2[3])
} )
##
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])
} )
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])
} )
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)
} )
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)
} )
##
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")
})
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")
})
setMethod("^",c("mks","mks"), function(e1,e2) stop("You cannot be serious"))
setMethod("^",c("fff","fff"), function(e1,e2) stop("You cannot be serious"))
setMethod("^",c("mks","fff"), function(e1,e2) stop("You cannot be serious"))
setMethod("^",c("fff","mks"), function(e1,e2) stop("You cannot be serious"))
##
setGeneric("root",function(e1,e2) standardGeneric("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")
})
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")
})
setMethod("sqrt", "mks", function(x) root(x,2))
setMethod("sqrt", "fff", function(x) root(x,2))
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])
} )
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])
} )
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)
} )
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)
} )
setMethod("Math","mks", function(x) {
if (isDimensionless(x))
callGeneric(getValues(x))
else stop("not defined for numbers with units")
})
setMethod("Math","fff", function(x) {
if (isDimensionless(x))
callGeneric(getValues(x))
else stop("not defined for numbers with units")
})
## not part of assignment, but could be useful
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")})
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment