diff --git a/NAMESPACE b/NAMESPACE
index 091250b1c8f4cbffceafa762627eb5cd6cf44505..9dca9b7de7bf205588d1f37293a6d74a6f7c2803 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -24,6 +24,12 @@ S3method(UWEDGE,default)
S3method(UWEDGE,rcomp)
S3method(accuracy,DataFrameStack)
S3method(accuracy,data.frame)
+S3method(as.AnisotropyRangeMatrix,AnisotropyRangeMatrix)
+S3method(as.AnisotropyRangeMatrix,AnisotropyScaling)
+S3method(as.AnisotropyRangeMatrix,default)
+S3method(as.AnisotropyScaling,AnisotropyRangeMatrix)
+S3method(as.AnisotropyScaling,AnisotropyScaling)
+S3method(as.AnisotropyScaling,numeric)
S3method(as.CompLinModCoReg,CompLinModCoReg)
S3method(as.CompLinModCoReg,LMCAnisCompo)
S3method(as.DataFrameStack,array)
diff --git a/R/gmAnisotropy.R b/R/gmAnisotropy.R
index 3e2ad18cde33c6be7cf371f82928866931562231..98202c706416ffeb65a87be612a9382423e5a09d 100644
--- a/R/gmAnisotropy.R
+++ b/R/gmAnisotropy.R
@@ -3,6 +3,30 @@
#### Anisotropy ------------------
+#' Convert to anisotropy scaling matrix
+#'
+#' Convert an anisotropy specification to a scaling matrix
+#'
+#' @param x an matrix to be tagged as anisotropy scaling matrix
+#'
+#' @return An anisotropy scaling matrix \eqn{A} is such that for any
+#' lag vector \eqn{h}, the variogram model turns isotropic in terms
+#' of \eqn{u'=h'\cdot A}. This function does not check any special
+#' property for this matrix! You should probably be using `anis_GSLIBpar2A()`
+#' isntead, and leave `AnisotropyScaling()` for internal uses.
+#'
+#'
+#' @export
+#' @family anisotropy
+#' @examples
+#' ( l = anis_GSLIBpar2A(angles=30, ratios=0.5))
+#' ( ll = unclass(l) )
+#' AnisotropyScaling(l)
+AnisotropyScaling = function(x){
+ class(x) = "AnisotropyScaling"
+ return(x)
+}
+
#' Convert to anisotropy scaling matrix
#'
#' Convert an anisotropy specification to a scaling matrix
@@ -21,31 +45,27 @@
#' are given in degrees, all ratios must be smaller or equal to 1. This follows gstat (and hence
#' GSlib) conventions; see gstat::vgm() for details.
#'
+#' @family anisotropy
#' @export
+#' @aliases as.AnisotropyScaling.numeric
#' @examples
#' ( l = anis_GSLIBpar2A(angles=30, ratios=0.5))
#' ( ll = unclass(l) )
-#' AnisotropyScaling(l)
-AnisotropyScaling = function(x){
- class(x) = "AnisotropyScaling"
- return(x)
-}
-
-#' @describeIn AnisotropyScaling coerce function to anisotropy scaling matrix
+#' as.AnisotropyScaling(ll)
#' @export
as.AnisotropyScaling <- function(x){ UseMethod("as.AnisotropyScaling", x) }
-# @describeIn AnisotropyScaling identity
-# @method as.AnisotropyScaling AnisotropyScaling
-# @export
+#' @describeIn as.AnisotropyScaling identity method
+#' @method as.AnisotropyScaling AnisotropyScaling
+#' @export
as.AnisotropyScaling.AnisotropyScaling = function(x){
x
}
-# @describeIn AnisotropyScaling from a vector of numbers
-# @method as.AnisotropyScaling numeric
-# @export
+#' @describeIn as.AnisotropyScaling from a vector of numbers
+#' @method as.AnisotropyScaling numeric
+#' @export
as.AnisotropyScaling.double <-as.AnisotropyScaling.numeric <- function(x){
if(length(dim(x))==0){
if(length(x)==2) return(anis2D_par2A(ratio=x[2], angle=x[1], inv=TRUE))
@@ -61,9 +81,9 @@ as.AnisotropyScaling.double <-as.AnisotropyScaling.numeric <- function(x){
}
}
-# @describeIn AnisotropyScaling from an AnisotropicRangeMatrix
-# @method as.AnisotropyScaling AnisotropyRangeMatrix
-# @export
+#' @describeIn as.AnisotropyScaling from an AnisotropicRangeMatrix
+#' @method as.AnisotropyScaling AnisotropyRangeMatrix
+#' @export
as.AnisotropyScaling.AnisotropyRangeMatrix = function(x){
AnisotropyScaling(solve(chol(x)))
}
@@ -81,8 +101,8 @@ as.AnisotropyScaling.AnisotropyRangeMatrix = function(x){
#' @param checkValidity boolean, should validity be checked?
#'
#' @return the same matrix with a class attribute
+#' @family anisotropy
#' @export
-#'
AnisotropyRangeMatrix = function(x, checkValidity=TRUE){
if(checkValidity){
if(length(dim(x))==3){
@@ -99,24 +119,35 @@ AnisotropyRangeMatrix = function(x, checkValidity=TRUE){
return(x)
}
-#' @describeIn AnisotropyRangeMatrix Convert to anisotropy range matrix
+#' Force a matrix to be anisotropy range matrix,
+#'
+#' Force a matrix M to be considered an anisotropy range matrix, i.e
+#' with ranges and orientations,
+#' such that \eqn{u = sqrt(h' * M^{-1} * h)} allows to use an isotropic
+#' variogram.
+#'
+#' @param x matrix simmetric positive definite (i.e. M above)
+#'
+#' @return the same anisotropy, specified as `M`
+#' @family anisotropy
#' @export
as.AnisotropyRangeMatrix <- function(x){ UseMethod("as.AnisotropyRangeMatrix", x) }
-# @describeIn AnisotropyRangeMatrix Convert to anisotropy range matrix
-# @method as.AnisotropyRangeMatrix default
-# @export
+
+#' @describeIn as.AnisotropyRangeMatrix Default conversion to anisotropy range matrix
+#' @method as.AnisotropyRangeMatrix default
+#' @export
as.AnisotropyRangeMatrix.default <- function(x) as.AnisotropyRangeMatrix(as.AnisotropyScaling(x))
-# @describeIn AnisotropyRangeMatrix Convert to anisotropy range matrix
-# @method as.AnisotropyRangeMatrix AnisotropyRangeMatrix
-# @export
+#' @describeIn as.AnisotropyRangeMatrix identity conversion
+#' @method as.AnisotropyRangeMatrix AnisotropyRangeMatrix
+#' @export
as.AnisotropyRangeMatrix.AnisotropyRangeMatrix <- function(x) x
-# @describeIn AnisotropyRangeMatrix Convert to anisotropy range matrix
-# @method as.AnisotropyRangeMatrix AnisotropyScaling
-# @export
+#' @describeIn AnisotropyRangeMatrix Convert from AnisotropyScaling
+#' @method as.AnisotropyRangeMatrix AnisotropyScaling
+#' @export
as.AnisotropyRangeMatrix.AnisotropyScaling <- function(x) AnisotropyRangeMatrix(solve(tcrossprod(x)), checkValidity=FALSE)
@@ -145,8 +176,9 @@ as.AnisotropyRangeMatrix.AnisotropyScaling <- function(x) AnisotropyRangeMatrix(
#'
#' Other values are meaningless.
#'
+#' @family anisotropy
#' @export
-anis_GSLIBpar2A = function(ratios, angles, inv=FALSE){
+ anis_GSLIBpar2A = function(ratios, angles, inv=FALSE){
if( (length(ratios)==1) & length(angles)==1){
anis2D_par2A(ratio=ratios, angle=angles, inv=inv)
}else if( (length(ratios)==1) & length(angles)==1){
@@ -196,7 +228,7 @@ anis2D_par2A = function(ratio, angle, inv=FALSE){
# rotation to new coordinates and scaling:
A = R %*% S # t(Rinv) %*% S # because t(Rinv) = R in orthogonal matrices
A = rbind(cbind(A,0),c(0,0,1))
- class(A) = c("AnisotropyScaling","Anisotropy")
+ class(A) = "AnisotropyScaling"
return(A)
}
@@ -240,6 +272,6 @@ anis3D_par2A = function(ratios, angles, inv=FALSE){
S = diag(c(1, ratios^((1-2*inv)/2) ))
# rotation to new coordinates and scaling:
A = R %*% S # t(Rinv) %*% S # because t(Rinv) = R in orthogonal matrices
- class(A) = c("AnisotropyScaling","Anisotropy")
+ class(A) = "AnisotropyScaling"
return(A)
}
diff --git a/R/variograms.R b/R/variograms.R
index 01502ace7235d15f473bb3dad2037ca77ec88155..d39a3f70a6c35d5de6648e36f818ef2e67cdd5c7 100644
--- a/R/variograms.R
+++ b/R/variograms.R
@@ -512,13 +512,14 @@ is.isotropic.LMCAnisCompo = function(x, tol=1e-10, ...){
#'
#' @return a logical, TRUE if the object is an anisotropy specification; FALSE otherwise
#' @export
+#' @family anisotropy
#'
#' @examples
#' a = anis2D_par2A(0.5, 30)
#' a
#' is.anisotropySpecification(a)
is.anisotropySpecification = function(x){
- "Anisotropy" %in% class(x) | inherits(x, "Anisotropy")
+ length(grep("Anisotropy", class(x))) | inherits(x, "Anisotropy")
}
diff --git a/vignettes/SimulatingMicrostructures.Rmd b/vignettes/SimulatingMicrostructures.Rmd
index a25d0dd9af691b727aa2da7f9e29628fd71ac82a..cadb809954f2b9b0467f4f946506791e88a29285 100644
--- a/vignettes/SimulatingMicrostructures.Rmd
+++ b/vignettes/SimulatingMicrostructures.Rmd
@@ -91,18 +91,18 @@ E.g. `angle = 45` with `ratio = 0.5` would result into NE-SW elongated crystals.
For Quartz, let's assume we have mostly roundish/quadratic grains and only a minor preferred orientation of the crystals, so that they are slightly oriented along NE-SW:
```{r p-qu}
-p1_qu = anis2D_par2A(ratio = .9, angle = 45)
+p1_qu = anis_GSLIBpar2A(ratios = .9, angles = 45)
```
The Feldspars have a more pronounced elongated shape, also with preferred orientation along the same direction as the Quartz grains:
```{r p-fsp}
-p2_fsp = anis2D_par2A(ratio = .6, angle = 45)
+p2_fsp = anis_GSLIBpar2A(ratios = .6, angles = 45)
```
-Function `anis2D_par2A` returns the Cholesky decomposition of the anisotropy matrix. So, if the anistropy of the mineral phases had been generated by this function, then we need to calculate the cross product of the matrices to obtain the anistropy matrices `A`:
+Function `anis_GSLIBpar2A` returns the inverse Cholesky decomposition of the anisotropy matrix. So, if the anistropy of the mineral phases had been generated by this function, then we need to re-calculate the anistropy range matrices `A`:
```{r}
-A1_qu = tcrossprod(p1_qu)
-A2_fsp = tcrossprod(p2_fsp)
+A1_qu = as.AnisotropyRangeMatrix(p1_qu)
+A2_fsp = as.AnisotropyRangeMatrix(p2_fsp)
```
This results into
@@ -121,20 +121,26 @@ There is also the option to directly generate the anisotropy matrix.
In this case, the matrix can be written manually:
The Biotite occurs as thin plates and has a slightly different preferred orientation then the two main minerals:
```{r}
-A3_bi = matrix(data = c(1, .6, 0, .6, 1, 0, 0, 0, 1), ncol = 3)
+A3_bi = AnisotropyRangeMatrix(
+ matrix(data = c(1, .6, 0, .6, 1, 0, 0, 0, 1), ncol = 3)
+ )
A3_bi
```
For the fracture and joint systems the anisotropy is very pronounced, because they are very thin but passing through longer distances through the rock.
Accordingly, the anisotropy matrix should be chosen similar to this:
```{r}
-A4_frac = matrix(data = c(1, -.999, 0, -.999, 1, 0, 0, 0, 1), ncol = 3)
+A4_frac = AnisotropyRangeMatrix(
+ matrix(data = c(1, -.999, 0, -.999, 1, 0, 0, 0, 1), ncol = 3)
+ )
A4_frac
```
And let's assume, that the joint system is orthogonal to the fracture system:
```{r}
-A5_joints = matrix(data = c(1, .999, 0, .999, 1, 0, 0, 0, 1), ncol = 3)
+A5_joints = AnisotropyRangeMatrix(
+ matrix(data = c(1, .999, 0, .999, 1, 0, 0, 0, 1), ncol = 3)
+ )
A5_joints
```
@@ -287,7 +293,7 @@ For example, in the simulation above, the geologist would expect, that the fract
Changing the orientation of the fractures:
```{r}
-p4_frac = anis2D_par2A(ratio = .01, angle = 95)
+p4_frac = anis_GSLIBpar2A(ratio = .01, angle = 95)
A4_frac = tcrossprod(p4_frac)
A4_frac
```