r - How to calculate the volume under a surface defined by discrete data? -


i need determine volume beneath series of surfaces represented discrete data points. in data, each sample stored separate data frame within list of data frames. here (small) example data:

df1 <- data.frame(x=c(2,2,2,3,3,3,4,4,4,5,5,5,6,6,6),                   y=c(1,2,3,1,2,3,1,2,3,1,2,3,1,2,3),                   z=c(0,2,0,4,6,7,3,2,1,2,7,8,9,4,2))  df2 <- data.frame(x=c(2,2,2,3,3,3,4,4,4,5,5,5,6,6,6),                   y=c(1,2,3,1,2,3,1,2,3,1,2,3,1,2,3),                   z=c(1,1,2,3,5,6,2,1,3,3,8,9,8,3,1))  df <- list(df1,df2)

answers similar questions either in other languages (matlab, python), or answers not contain useable script address problem (as here). can think of 2 acceptable ways estimate volume beneath each surface: 1) write out discretized version of simpson's rule function in r applied across list of data frames (df); 2) calculate arbitrary relationship between x, y, , z , use multivariate numerical integration find volume under surface (with functions simpson2d / quad2d in package pracma or adaptintegrate in cubature).

regarding first approach, formula composite simpson's rule (that use) here, due complexity, have been unsuccessful in writing working double summation function. in expression, i(lambda(em) lambda(ex)) equal z in above datasets @ each x,y grid point, , delta(em) , delta(ex) represent interval between x , y points.

the second approach extend approach found here multivariate spline fits , pass predicted z values function integration. here's have tried far approach:

require(pracma)  df1.loess <- loess(z ~ x + y, data=df[[1]]) mod.fun <- function(x,y) predict(df1.loess, newdata=x,y)  simpson2d(mod.fun, x=c(2,6), y=c(1,3))

but not yield useful results.

in reality, have list of 100 data frames individual samples, need able express solution series of lapply functions automate these calculations across data frames in list. example looks this:

require(akima) df.splines <- lapply(df, function(x,y,z) interp(x = "x", y = "y", z = "z",                                                 linear=f, nx=4, ny=2))

unfortunately, produces exception missing values , infs. i'm extremely open suggestions how implement 1 of these strategies, or utilize different (simpler?) approach. kriging function (like km in dicekriging package) produce better fit passed on numerical integration?

i assuming volume surface mesh defined connecting points via straight lines. can find volume beneath surface via

  1. triangular tessellation of (x,y) grid triangles t_i area a_i
  2. finding corresponding z values z_i each of triangles t_i
  3. calculating volume v_i of truncated prisms (defined t_i , z_i) via v_i=a_i*sum(z_i)/3 (see https://en.wikipedia.org/wiki/prism_(geometry) , https://math.stackexchange.com/questions/2371139/volume-of-truncated-prism)
  4. summing truncated prism volumes v_i

keep in mind, however, volume depend on tessellation , tessellation not unique. problem not defined in sense not describe how 1 should interpolate between points. approach calculate volume have make additional assumptions.

going solution approach, points 1 , 2 can achieved via geometry package. here code

library(geometry)  getvolume=function(df) {   #find triangular tesselation of (x,y) grid   res=delaunayn(as.matrix(df[,-3]),full=true,options="qz")   #calulates sum of truncated prism volumes   sum(mapply(function(tripoints,a) a/3*sum(df[tripoints,"z"]),              split.data.frame(res$tri,seq_along(res$areas)),              res$areas)) }  sapply(df,getvolume) #[1] 32.50000 30.33333

since it's hard check whether results consistent, here simple example know right answer. it's cube side length 2 have cut out wedge along x axis. cut-out region 1/4 of total volume.

cutoutcube=expand.grid(c(0,1,2),c(0,1,2)) colnames(cutoutcube)=c("x","y") cutoutcube$z=ifelse(cutoutcube$x==1,1,2)  sapply(list(cutoutcube),getvolume) #[1] 6

that's correct since 2^3*(1-1/4)=6.

another sanity check can performed calculating "complement" of volume w.r.t. simple cuboid z values set max z value (in case max(z)=9 in both cases). simple cuboid volumes 72 both of cases. not let's define complement surfaces , sum volume , complement volume

df1c=df1 df1c$z=max(df1c$z)-df1c$z df2c=df2 df2c$z=max(df2c$z)-df2c$z dfc=list(df1c,df2c)  sapply(dfc,getvolume)+sapply(df,getvolume) #[1] 72 72

so volume , complement volume give right simple cuboid volume in both cases.


Comments

Post a Comment

Popular posts from this blog

Is there a better way to structure post methods in Class Based Views -

i'm using class based views in django , never specify url in ajax calls , omit action parameter in form because know these requests go through post method of class based view. notice view - profile page make many requests depending on user updating on page - post method become convoluted , if/else mess example: def post(self, request, *args, **kwargs): # handle user changing available date if request.post.get('availabledate') != none: self.updateavailabledate(request) return jsonresponse({'result':'success'}) if request.post.get('newprojectname') != none: creative_user = creativeuserprofile.objects.get(id = request.user.id) project = project.create_project(creative_user, request.post['newprojectname']) project.save() return jsonresponse({'projectid': project.id }) if request.post.get('projectimage') != none , request.files['file']: proj...
Read more

performance - Why is XCHG reg, reg a 3 micro-op instruction on modern Intel architectures? -

i'm doing micro-optimization on performance critical part of code , came across sequence of instructions (in at&t syntax): add %rax, %rbx mov %rdx, %rax mov %rbx, %rdx i thought had use case xchg allow me shave instruction , write: add %rbx, %rax xchg %rax, %rdx however, dimay found agner fog's instruction tables , xchg 3 micro-op instruction 2 cycle latency on sandy bridge, ivy bridge, broadwell, haswell , skylake. 3 whole micro-ops , 2 cycles of latency! 3 micro-ops throws off 4-1-1-1 cadence , 2 cycle latency makes worse original in best case since last 2 instructions in original might execute in parallel. now... cpu might breaking instruction micro-ops equivalent to: mov %rax, %tmp mov %rdx, %rax mov %tmp, %rdx where tmp anonymous internal register , suppose last 2 micro-ops run in parallel latency 2 cycles. given register renaming occurs on these micro-architectures, though, doesn't make sense me done way. why wouldn't register renamer ...
Read more

c# - Asp.net web api : redirect unauthorized requst to forbidden page -

im trying redirect unauthorized request forbidden page instead i'm getting forbidden page in response body , how can fix ? here's startup class : app.createperowincontext(storecontext.create); app.createperowincontext<applicationusermanager>(applicationusermanager.create); app.createperowincontext<applicationsigninmanager>(applicationsigninmanager.create); app.usecookieauthentication(new cookieauthenticationoptions { authenticationtype = defaultauthenticationtypes.applicationcookie, expiretimespan = timespan.fromdays(30), }); app.useexternalsignincookie(defaultauthenticationtypes.externalcookie); app.usetwofactorsignincookie(defaultauthenticationtypes.twofactorcookie, timespan.fromminutes(5)); app.usetwofactorrememberbrowsercookie(defaultauthenticationtypes.twofactorrememberbrowsercookie); the method im trying reach : [httpget] [authorize(roles = "admin")] public string getcurrentusername() { return userman...
Read more