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
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?
Thanks,
Alex
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
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
>
>
_______________________________________________
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
