On 15 December 2021 at 11:22, Alex Ilich wrote: | I actually think that the issue may be related to conversions between | int/IntegerVectors and double/NumericVector/std::vector<double>. I was | talking with the terra package author to see if the focalCpp function can | take an IntegerVector instead of a NumericVector/ /std::vector<double>, as | it seemed initially to be working and he said that Rcpp automatically | coerces a std::vector<double> to an Rcpp::IntegerVector. However, when I
That is a misunderstanding. That would be a *major* bug. Brining Rob in just in case. Rcpp faithfully and reliably maps back and forth between std::vector<int32_t> and Rcpp::IntegerVector std::vector<double> and Rcpp::NumericVector as well as other numerics and of course std::string etc What you may be confusing is that writing something such as > Rcpp::cppFunction("int mysum(IntegerVector x) { return sum(x); }") > mysum(c(1.0, 2.0, 3.0)) [1] 6 > Rcpp::cppFunction("double myothersum(NumericVector x) { return sum(x); }") > myothersum(1:3) # the : creates an _int_ sequence [1] 6 > myothersum(c(1L,2L,3L)) [1] 6 > cast up/down from/to int and double _when access from R_ which is a _very_ different story. | replace the call to C_make_glcm which was causing the issue to just | out(i,0) = xw[0] (i.e. return the first value), the edges that should be NA | evaluate to -2147483648 so maybe this automatic conversion doesn't work if | there are NAs. I can reproduce this using only Rcpp by having a function | that takes an IntegerVector and returns the first value of that vector as a | double, in which case NA's are changed to -2147483648 (code below). Is this | expected behavior, and are there any ways to guard against this? It's a bug in your code as you don't cast right. And you also 'call it incorrectly' just as per my example above. NA is a _double_ as is 1. We then help you by making the _now int_ value a matching NA _as an int_. And here is a *very* important bit: Int NAs are an R-only thing. The assignment you make is a C/C++ one *which does not know this* and you loose NAness. No more, no less. Just a usage bug on your end. | library(Rcpp) | cppFunction('double C_fun2(IntegerVector x){ | size_t start = 0; | size_t end = x.length()-1; | IntegerVector xw = x[Rcpp::Range(start,end)]; //Current window of | elevation values | double out=xw[0]; | return(out); | }') | | C_fun2(c(NA,1)) | # -2147483648 Simpler fixed version: > Rcpp::cppFunction("double alex(IntegerVector x) { NumericVector y(x); return y[0];}") > alex(c(NA, 1)) [1] NA > alex(c(NA_integer_, 1L)) # "correct" call with actual int arguments [1] NA > This ensure we build a proper R-context aware vector of doubles and then return its first element. You can use the same logic in your code, we also do this (I think) when using std::vector<> of these base types as they are so common. Dirk | On Mon, Dec 13, 2021 at 3:43 PM Alex Ilich <alexilic...@gmail.com> wrote: | | > Thanks for digging into this for me. If focal_val is 16, then this is an | > out of bounds error as the matrix is 16x16 and the maximum index can be 15. | > What I'm confused about is that focal_val comes from the data which has | > been transformed to be restricted to integers in the range of 0-15 (0 - | > n_levels-1). This is explicitly checked for in the glcm_textures2 function | > with the line " if((terra::global(r, fun = max) > (n_levels-1)) | | > (terra::global(r, fun = min) < 0)){stop("Error: raster must have values | > between 0 and n_levels-1")}". Additionally, if I start with numbers in 0-15 | > and don't do any transformation (quantization="none"), I still run into | > this issue (see code below). | > | > Thanks, | > Alex | > | > library(terra) | > library(raster) | > library(GLCMTextures) | > | > r1b<- terra::rast(volcano) | > set.seed(5) | > values(r1b)<- sample(0:15, size = ncell(r1b), replace = TRUE) # Replace | > values with random integers from 0-15 | > r2b<- glcm_textures2(r1b, w=c(3,7), n_levels = 16, quantization = "none", | > shift=c(0,1)) #Often leads to Rsession Aborted | > | > On Mon, Dec 13, 2021 at 3:12 PM Bill Dunlap <williamwdun...@gmail.com> | > wrote: | > | >> In this example vgdb shows that valgrind's first complaint is at | >> (gdb) where 5 | >> #0 0x000000001c1ada01 in C_make_glcm (x=..., n_levels=16, shift=..., | >> na_opt=...) at glcm_cpp_functions.cpp:73 | >> #1 0x000000001c1b2238 in C_glcm_textures_helper2 (x=..., w2=..., | >> n_levels=16, shift=..., na_opt=..., ni=1, nw=21) | >> at glcm_cpp_functions.cpp:218 | >> #2 0x000000001c1a098e in _GLCMTextures_C_glcm_textures_helper2 | >> (xSEXP=0x19023168, w2SEXP=0x92cbde8, | >> n_levelsSEXP=0xa3138d8, shiftSEXP=0x92cc468, na_optSEXP=0xa313398, | >> niSEXP=0xb655560, nwSEXP=0x1aeec7d0) | >> at RcppExports.cpp:94 | >> #3 0x00000000001f6b02 in R_doDotCall ( | >> ofun=0x1c1a0732 <_GLCMTextures_C_glcm_textures_helper2(SEXP, SEXP, | >> SEXP, SEXP, SEXP, SEXP, SEXP)>, nargs=7, | >> cargs=0x1ffeff6830, call=0x15ed3c88) at | >> /mnt/c/R/R-svn/trunk/src/main/dotcode.c:622 | >> #4 0x00000000002637ca in bcEval (body=0x15ed4078, rho=0x19023210, | >> useCache=TRUE) | >> at /mnt/c/R/R-svn/trunk/src/main/eval.c:7695 | >> | >> Line 73 is | >> GLCM(neighbor_val,focal_val) = GLCM(neighbor_val,focal_val)+1; | >> where GLCM is a 16 x 16 IntegerMatrix and focal_val is 16: | >> (gdb) print neighbor_val | >> $1 = 9 | >> (gdb) print focal_val | >> $2 = 16 | >> (gdb) print i | >> $3 = 2 | >> (gdb) print j | >> $4 = 1 | >> | >> -Bill | >> | >> On Mon, Dec 13, 2021 at 11:37 AM Bill Dunlap <williamwdun...@gmail.com> | >> wrote: | >> | >>> Thanks, I hadn't noticed the --use-valgrind option to check. | >>> | >>> I often run R under valgrind with the command line: | >>> $ R --quiet --no-save --debugger=valgrind | >>> --debugger-args="--track-origins=yes --vgdb=full --vgdb-error=0" | >>> then in another window start a debugger process with | >>> $ gdb RHOME/bin/exec/R | >>> ... | >>> (gdb) target remote | vgdb | >>> (gdb) cont | >>> You may need to supply a complete path to vgdb, valgrind's interface to | >>> gdb. | >>> Then gdb will set breakpoints where valgrind reports memory misuse and | >>> you can inspect variables, etc., where an error occurs. | >>> You can check for memory leaks from gdb with | >>> (gdb) monitor leak_check full | >>> | >>> -Bill | >>> | >>> | >>> On Mon, Dec 13, 2021 at 10:29 AM Dirk Eddelbuettel <e...@debian.org> | >>> wrote: | >>> | >>>> | >>>> On 13 December 2021 at 09:15, Bill Dunlap wrote: | >>>> | I ran your example under valgrind on Linux (Ubuntu 20.04) | >>>> | >>>> Excellent call, as usual! Thanks for doing that. | >>>> | >>>> R actually makes is so easy to run with valgrind by calling | >>>> | >>>> R CMD check --use-valgrind somePkg_1.2-3.tar.gz | >>>> | >>>> that I started to add it to some CI setups. | >>>> | >>>> Dirk | >>>> | >>>> | and valgrind found some memory misuse: | >>>> | | >>>> | $ R --quiet --no-save --debugger=valgrind | >>>> | --debugger-args="--track-origins=yes" | >>>> | ==1533== Memcheck, a memory error detector | >>>> | ==1533== Copyright (C) 2002-2017, and GNU GPL'd, by Julian Seward et | >>>> al. | >>>> | ==1533== Using Valgrind-3.15.0 and LibVEX; rerun with -h for | >>>> copyright info | >>>> | ==1533== Command: /home/bill/R-devel/R-build/bin/exec/R --quiet | >>>> --no-save | >>>> | ==1533== | >>>> | > library(terra) | >>>> | terra version 1.5.2 | >>>> | > library(raster) | >>>> | Loading required package: sp | >>>> | > library(GLCMTextures) | >>>> | > | >>>> | > r1a<- raster::raster(volcano) | >>>> | > r2a<- glcm_textures(r1a, w=c(3,7), n_levels = 16, quantization = | >>>> "equal | >>>> | prob", shift=c(0,1)) | >>>> | > | >>>> | > r1b<- terra::rast(volcano) | >>>> | > r2b<- glcm_textures2(r1b, w=c(3,7), n_levels = 16, quantization = | >>>> "equal | >>>> | prob", shift=c(0,1)) #Often leads to Rsession Aborted | >>>> | ==1533== Invalid read of size 4 | >>>> | ==1533== at 0x1C1ADA01: C_make_glcm(Rcpp::Matrix<13, | >>>> | Rcpp::PreserveStorage>, int, Rcpp::Vector<13, Rcpp::PreserveStorage>, | >>>> | Rcpp::String) (glcm_cpp_functions.cpp:73) | >>>> | ==1533== by 0x1C1B2237: C_glcm_textures_helper2(Rcpp::Vector<13, | >>>> | Rcpp::PreserveStorage>, Rcpp::Vector<13, Rcpp::PreserveStorage>, int, | >>>> | Rcpp::Vector<13, Rcpp::PreserveStorage>, Rcpp::String, unsigned long, | >>>> | unsigned long) (glcm_cpp_functions.cpp:218) | >>>> | ==1533== by 0x1C1A098D: _GLCMTextures_C_glcm_textures_helper2 | >>>> | (RcppExports.cpp:94) | >>>> | ==1533== by 0x1F6B01: R_doDotCall (dotcode.c:622) | >>>> | ==1533== by 0x2637C9: bcEval (eval.c:7695) | >>>> | ==1533== by 0x240653: Rf_eval (eval.c:748) | >>>> | ==1533== by 0x243223: R_execClosure (eval.c:1918) | >>>> | ==1533== by 0x242EFD: Rf_applyClosure (eval.c:1844) | >>>> | ==1533== by 0x240E06: Rf_eval (eval.c:871) | >>>> | ==1533== by 0x246CD6: do_set (eval.c:2990) | >>>> | ==1533== by 0x240A8C: Rf_eval (eval.c:823) | >>>> | ==1533== by 0x2457DA: do_begin (eval.c:2538) | >>>> | ==1533== Address 0xa59c0c4 is 28 bytes before a block of size 64 in | >>>> arena | >>>> | "client" | >>>> | ==1533== | >>>> | ==1533== Invalid write of size 4 | >>>> | ==1533== at 0x1C1ADA44: C_make_glcm(Rcpp::Matrix<13, | >>>> | Rcpp::PreserveStorage>, int, Rcpp::Vector<13, Rcpp::PreserveStorage>, | >>>> | Rcpp::String) (glcm_cpp_functions.cpp:73) | >>>> | ==1533== by 0x1C1B2237: C_glcm_textures_helper2(Rcpp::Vector<13, | >>>> | Rcpp::PreserveStorage>, Rcpp::Vector<13, Rcpp::PreserveStorage>, int, | >>>> | Rcpp::Vector<13, Rcpp::PreserveStorage>, Rcpp::String, unsigned long, | >>>> | unsigned long) (glcm_cpp_functions.cpp:218) | >>>> | ==1533== by 0x1C1A098D: _GLCMTextures_C_glcm_textures_helper2 | >>>> | (RcppExports.cpp:94) | >>>> | ==1533== by 0x1F6B01: R_doDotCall (dotcode.c:622) | >>>> | ==1533== by 0x2637C9: bcEval (eval.c:7695) | >>>> | ==1533== by 0x240653: Rf_eval (eval.c:748) | >>>> | ==1533== by 0x243223: R_execClosure (eval.c:1918) | >>>> | ==1533== by 0x242EFD: Rf_applyClosure (eval.c:1844) | >>>> | ==1533== by 0x240E06: Rf_eval (eval.c:871) | >>>> | ==1533== by 0x246CD6: do_set (eval.c:2990) | >>>> | ==1533== by 0x240A8C: Rf_eval (eval.c:823) | >>>> | ==1533== by 0x2457DA: do_begin (eval.c:2538) | >>>> | ==1533== Address 0xa59c0c4 is 28 bytes before a block of size 64 in | >>>> arena | >>>> | "client" | >>>> | ==1533== | >>>> | --1533-- VALGRIND INTERNAL ERROR: Valgrind received a signal 11 | >>>> (SIGSEGV) - | >>>> | exiting | >>>> | --1533-- si_code=1; Faulting address: 0x10A59C138; sp: 0x10090e2e20 | >>>> | | >>>> | valgrind: the 'impossible' happened: | >>>> | Killed by fatal signal | >>>> | | >>>> | host stacktrace: | >>>> | ==1533== at 0x580505C4: ??? (in | >>>> | /usr/lib/x86_64-linux-gnu/valgrind/memcheck-amd64-linux) | >>>> | ==1533== by 0x58004EBB: ??? (in | >>>> | /usr/lib/x86_64-linux-gnu/valgrind/memcheck-amd64-linux) | >>>> | ==1533== by 0x58005DA7: ??? (in | >>>> | /usr/lib/x86_64-linux-gnu/valgrind/memcheck-amd64-linux) | >>>> | ==1533== by 0x580A7204: ??? (in | >>>> | /usr/lib/x86_64-linux-gnu/valgrind/memcheck-amd64-linux) | >>>> | ==1533== by 0x580F5FD4: ??? (in | >>>> | /usr/lib/x86_64-linux-gnu/valgrind/memcheck-amd64-linux) | >>>> | | >>>> | sched status: | >>>> | running_tid=1 | >>>> | | >>>> | Thread 1: status = VgTs_Runnable (lwpid 1533) | >>>> | ==1533== at 0x483BE63: operator new(unsigned long) (in | >>>> | /usr/lib/x86_64-linux-gnu/valgrind/vgpreload_memcheck-amd64-linux.so) | >>>> | ==1533== by 0x6E92D3C: void std::__cxx11::basic_string<char, | >>>> | std::char_traits<char>, std::allocator<char> >::_M_construct<char | >>>> | const*>(char const*, char const*, std::forward_iterator_tag) (in | >>>> | /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.28) | >>>> | ==1533== by 0x1C1B2377: C_glcm_textures_helper2(Rcpp::Vector<13, | >>>> | Rcpp::PreserveStorage>, Rcpp::Vector<13, Rcpp::PreserveStorage>, int, | >>>> | Rcpp::Vector<13, Rcpp::PreserveStorage>, Rcpp::String, unsigned long, | >>>> | unsigned long) (glcm_cpp_functions.cpp:221) | >>>> | ==1533== by 0x1C1A098D: _GLCMTextures_C_glcm_textures_helper2 | >>>> | (RcppExports.cpp:94) | >>>> | ==1533== by 0x1F6B01: R_doDotCall (dotcode.c:622) | >>>> | ==1533== by 0x2637C9: bcEval (eval.c:7695) | >>>> | ==1533== by 0x240653: Rf_eval (eval.c:748) | >>>> | ==1533== by 0x243223: R_execClosure (eval.c:1918) | >>>> | ==1533== by 0x242EFD: Rf_applyClosure (eval.c:1844) | >>>> | ==1533== by 0x240E06: Rf_eval (eval.c:871) | >>>> | ==1533== by 0x246CD6: do_set (eval.c:2990) | >>>> | ==1533== by 0x240A8C: Rf_eval (eval.c:823) | >>>> | ==1533== by 0x2457DA: do_begin (eval.c:2538) | >>>> | ==1533== by 0x240A8C: Rf_eval (eval.c:823) | >>>> | ==1533== by 0x244FC9: do_for (eval.c:2420) | >>>> | ==1533== by 0x240A8C: Rf_eval (eval.c:823) | >>>> | ==1533== by 0x2457DA: do_begin (eval.c:2538) | >>>> | ==1533== by 0x240A8C: Rf_eval (eval.c:823) | >>>> | ==1533== by 0x244FC9: do_for (eval.c:2420) | >>>> | ==1533== by 0x240A8C: Rf_eval (eval.c:823) | >>>> | ==1533== by 0x2457DA: do_begin (eval.c:2538) | >>>> | ==1533== by 0x240A8C: Rf_eval (eval.c:823) | >>>> | ==1533== by 0x243223: R_execClosure (eval.c:1918) | >>>> | ==1533== by 0x242EFD: Rf_applyClosure (eval.c:1844) | >>>> | ==1533== by 0x253F32: bcEval (eval.c:7107) | >>>> | ==1533== by 0x240653: Rf_eval (eval.c:748) | >>>> | ==1533== by 0x243223: R_execClosure (eval.c:1918) | >>>> | ==1533== by 0x243D9C: R_execMethod (eval.c:2094) | >>>> | ==1533== by 0x8800187: R_dispatchGeneric | >>>> (methods_list_dispatch.c:1145) | >>>> | ==1533== by 0x2B15F4: do_standardGeneric (objects.c:1285) | >>>> | ==1533== by 0x253DB3: bcEval (eval.c:7096) | >>>> | ==1533== by 0x240653: Rf_eval (eval.c:748) | >>>> | ==1533== by 0x243223: R_execClosure (eval.c:1918) | >>>> | ==1533== by 0x242EFD: Rf_applyClosure (eval.c:1844) | >>>> | ==1533== by 0x253F32: bcEval (eval.c:7107) | >>>> | ==1533== by 0x240653: Rf_eval (eval.c:748) | >>>> | ==1533== by 0x243223: R_execClosure (eval.c:1918) | >>>> | ==1533== by 0x242EFD: Rf_applyClosure (eval.c:1844) | >>>> | ==1533== by 0x240E06: Rf_eval (eval.c:871) | >>>> | ==1533== by 0x246CD6: do_set (eval.c:2990) | >>>> | ==1533== by 0x240A8C: Rf_eval (eval.c:823) | >>>> | ==1533== by 0x28EF5E: Rf_ReplIteration (main.c:264) | >>>> | ==1533== by 0x28F157: R_ReplConsole (main.c:316) | >>>> | ==1533== by 0x290BC0: run_Rmainloop (main.c:1130) | >>>> | ==1533== by 0x290BDA: Rf_mainloop (main.c:1137) | >>>> | ==1533== by 0x16F74B: main (Rmain.c:29) | >>>> | client stack range: [0x1FFEED3000 0x1FFF000FFF] client SP: | >>>> 0x1FFEFF35F0 | >>>> | valgrind stack range: [0x1008FE3000 0x10090E2FFF] top usage: 12336 of | >>>> | 1048576 | >>>> | | >>>> | | >>>> | Note: see also the FAQ in the source distribution. | >>>> | It contains workarounds to several common problems. | >>>> | In particular, if Valgrind aborted or crashed after | >>>> | identifying problems in your program, there's a good chance | >>>> | that fixing those problems will prevent Valgrind aborting or | >>>> | crashing, especially if it happened in m_mallocfree.c. | >>>> | | >>>> | If that doesn't help, please report this bug to: www.valgrind.org | >>>> | | >>>> | In the bug report, send all the above text, the valgrind | >>>> | version, and what OS and version you are using. Thanks. | >>>> | | >>>> | On Mon, Dec 13, 2021 at 8:15 AM Alexander Ilich <ail...@mail.usf.edu> | >>>> wrote: | >>>> | | >>>> | > Hi, I'm upgrading one of my R packages to rely on the terra package | >>>> | > instead of the raster package for the handling of spatial raster | >>>> data. | >>>> | > Previously when relying on the raster package I had to convert the | >>>> data to | >>>> | > a matrix and send it to C++ to loop through the cells manually, but | >>>> with | >>>> | > the new features in the terra package I can supply a C++ function | >>>> that | >>>> | > returns multiple values directly to terra::focalCpp to perform focal | >>>> | > operations | >>>> | > < | >>>> https://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-focal-statistics-works.htm | >>>> >. | >>>> | > I get the same values using both the old and new versions of the | >>>> functions, | >>>> | > but the new version often causes the R session to abort, especially | >>>> at | >>>> | > larger window sizes. For example, a 3x3 or a 3x5 window seems to | >>>> always | >>>> | > run, but a 3x7 window will often cause the R session to abort. | >>>> However, | >>>> | > when it doesn't, I get the correct values. Functions followed by 2 | >>>> are the | >>>> | > terra versions of the function. By commenting out things and | >>>> rebuilding I | >>>> | > was able to determine the line that causes the crash in the terra | >>>> version | >>>> | > is "NumericMatrix curr_GLCM = C_make_glcm(curr_window, n_levels, | >>>> shift, | >>>> | > na_opt); //Tabulate the GLCM"; however, this line is also included | >>>> in the | >>>> | > raster version of the function so I'm not sure why this would | >>>> happen. Any | >>>> | > help would be greatly appreciated. Here is the github repository | >>>> | > <https://github.com/ailich/GLCMTextures/tree/terra>, and I've | >>>> added some | >>>> | > sample code below to illustrate the issue. | >>>> | > | >>>> | > Thanks, | >>>> | > Alex | >>>> | > | >>>> | > install.packages('raster', repos='https://rspatial.r-universe.dev') | >>>> | > #install development version of raster | >>>> | > install.packages('terra', repos='https://rspatial.r-universe.dev') | >>>> | > #install development version of terra | >>>> | > | >>>> | > remotes::install_github("ailich/GLCMTextures", ref = "terra") | >>>> #Install | >>>> | > branch of my package testing terra versions of functions | >>>> | > | >>>> | > library(terra) | >>>> | > library(raster) | >>>> | > library(GLCMTextures) | >>>> | > | >>>> | > r1a<- raster::raster(volcano) | >>>> | > r2a<- glcm_textures(r1a, w=c(3,7), n_levels = 16, quantization = | >>>> "equal | >>>> | > prob", shift=c(0,1)) | >>>> | > | >>>> | > r1b<- terra::rast(volcano) | >>>> | > r2b<- glcm_textures2(r1b, w=c(3,7), n_levels = 16, quantization = | >>>> "equal | >>>> | > prob", shift=c(0,1)) #Often leads to Rsession Aborted | >>>> | > | >>>> | > all.equal(values(r2a),values(r2b)) #TRUE | >>>> | > | >>>> | > | >>>> | > | >>>> | > #System Information | >>>> | > #OS: Windows 10 | >>>> | > R.version | >>>> | > # platform x86_64-w64-mingw32 | >>>> | > # arch x86_64 | >>>> | > # os mingw32 | >>>> | > # system x86_64, mingw32 | >>>> | > # status | >>>> | > # major 4 | >>>> | > # minor 0.4 | >>>> | > # year 2021 | >>>> | > # month 02 | >>>> | > # day 15 | >>>> | > # svn rev 80002 | >>>> | > # language R | >>>> | > # version.string R version 4.0.4 (2021-02-15) | >>>> | > # nickname Lost Library Book | >>>> | > | >>>> | > packageVersion("terra") | >>>> | > # ‘1.5.2’ | >>>> | > packageVersion("raster") | >>>> | > # ‘3.5.10’ | >>>> | > packageVersion("Rcpp") | >>>> | > # ‘1.0.7’ | >>>> | > packageVersion("RcppArmadillo") | >>>> | > # ‘0.10.7.3.0’ | >>>> | > _______________________________________________ | >>>> | > Rcpp-devel mailing list | >>>> | > Rcpp-devel@lists.r-forge.r-project.org | >>>> | > | >>>> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel | >>>> | _______________________________________________ | >>>> | Rcpp-devel mailing list | >>>> | Rcpp-devel@lists.r-forge.r-project.org | >>>> | | >>>> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel | >>>> -- | >>>> https://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org | >>>> | >>> _______________________________________________ | >> Rcpp-devel mailing list | >> Rcpp-devel@lists.r-forge.r-project.org | >> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel | > | > -- https://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org _______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel