'''
Module which contains the required functions to apply the statistical model for
a sample of :math:`n` second-order tensors
`Jelínek (1978) <https://doi.org/10.1007/BF01613632>`_ in order to obtain the
mean tensor :math:`\\boldsymbol{k}` of the sample, the
:math:`\\boldsymbol{k}`'s principal values :math:`k_1, k_2\ \\&\ k_3`, with
their confidence intervals, and the :math:`\\boldsymbol{k}`'s principal
directions :math:`\\boldsymbol{p}_1, \\boldsymbol{p}_2\ \&\ \\boldsymbol{p}_3`
with their confidence regions.
This application program is able to plot the summary of the statistical model
described above in a stereographic projection for a better understanding of the
outcomes.
Note:
* The packages `numpy <http://www.numpy.org/>`_,\
`scipy <https://www.scipy.org/>`_,\
`matplotlib <https://matplotlib.org/>`_\
and `mplstereonet <https://pypi.python.org/pypi/mplstereonet>`_ are\
required for using the ``jelinekstat.py`` module. All of them are\
downloadable from the PyPI repository.
* The mathematical notation in this documentation is taken from the\
original reference\
`Jelínek (1978) <https://doi.org/10.1007/BF01613632>`_.
* Copyright (c) 2018, Universidad Nacional de Colombia, Medellín.\
Copyright (c) 2018, Exneyder A. Monotoya-Araque and Ludger O.\
Suarez-Burgoa.\
`BSD-2-Clause <https://opensource.org/licenses/BSD-2-Clause>`_ or\
higher.
'''
[docs]def normalizeTensors(sample):
'''Divides all the tensor's elements by the mean susceptibility
:math:`{k}`, *i.e.* gets :math:`\\boldsymbol{k}_\\mathrm{norm}` using the
equations (8) of `Jelínek (1978) <https://doi.org/10.1007/BF01613632>`_.
Parameters:
sample (`numpy.ndarray`): :math:`(n \\times 6)` array that
contains the values obtained from the ``extractdata`` fuction.
Returns:
(`numpy.ndarray`): :math:`(n \\times 6)` array that contains the\
tensors :math:`\\boldsymbol{k}_\\mathrm{norm}` with the same\
format and structure of the ``extractdata`` function's output.
Examples:
>>> from jelinekstat.tools import dataFromFile
>>> sample, numTensors = dataFromFile('inputDataExample.txt')
>>> normalizeTensors(sample)
array(
[[ 1.02327, 1.02946, 0.94727, -0.01495, -0.03599, -0.05574],
[1.02315, 1.01803, 0.95882, -0.00924, -0.02058, -0.03151],
[1.02801, 1.03572, 0.93627, -0.03029, -0.03491, -0.06088],
[1.02775343, 1.00633335, 0.96591322, -0.01635005, -0.04148014,
-0.02006007],
[1.02143, 1.01775, 0.96082, -0.02798, -0.04727, -0.02384],
[1.01822661, 1.01202663, 0.96974677, -0.01125996, -0.02832991,
-0.03648988],
[1.01486338, 1.0206734 , 0.96446321, -0.01046003, -0.01913006,
-0.03864013],
[1.04596, 1.01133, 0.94271, -0.0166 , -0.04711, -0.03636]])
'''
import numpy as np
# Knowing the number of rows
sample = np.array(sample)
numTensors = int(np.shape(sample)[0])
# Normalizing the tensors with the first invariant.
sampleNorm = np.zeros((0, 6))
for k in range(0, numTensors):
a = sample[[k], :] / sum(sample[[k], 0:3][0])
b = np.array(list(a[0])) * 3
sampleNorm = np.vstack([sampleNorm, b])
return sampleNorm
[docs]def meantensor(sample, normalized=False):
'''Estimates the mean tensor :math:`\\boldsymbol{k}` of a randomly chosen
sample of :math:`n` specimens by using the equation (11) of
`Jelínek (1978) <https://doi.org/10.1007/BF01613632>`_ after being
normalized the specimens through the ``normalizeTensors`` function.
Parameters:
sample (`numpy.ndarray`): :math:`(n \\times 6)` array that
contains the values of the tensors after being imported with the
``extractdata`` function.
normalize (`bool`): Logical variable to indicate if the tensors in
the ``sample`` variable are already normalized by using the
equation (11) of
(`Jelínek (1978) <https://doi.org/10.1007/BF01613632>`_). ``False``
is the default value. In the case they are not normalized,
they will be.
Returns:
Three elements are returned; they are described below.
- **meanTensorVect** (`numpy.ndarray`): Mean tensor in vector\
form.
- **meanTensorMtx** (`numpy.ndarray`): Mean tensor in matrix\
form.
- **numTensors** (`int`): Number of tensors.
Examples:
>>> from jelinekstat.tools import dataFromFile
>>> from jelinekstat.jelinekstat import meantensor
>>> sample, numTensors = dataFromFile('inputDataExample.txt')
>>> meanTensorVect, meanTensorMtx, numTensors = meantensor(
sample, normalized=False)
>>> meanTensorVect
array([ 1.02533293, 1.01891542, 0.95575165, -0.01714126, -0.03435001,
>>> -0.03794001])
>>> meanTensorMtx
array([[ 1.02533293, -0.01714126, -0.03794001],
[-0.01714126, 1.01891542, -0.03435001],
[-0.03794001, -0.03435001, 0.95575165]])
>>> numTensors
8
'''
import numpy as np
from jelinekstat.tools import tensorvect2matrixform
sample = sample.T
numTensors = np.shape(sample)[1]
tensorsList = list()
normalizedTensorsList = list()
for i in range(numTensors):
T = sample[:, i]
tensorsList.append(T)
if normalized:
normalizedTensorsList.append(T)
else: # normalization with respect first invariant
inv1 = sum(T[0:3])/3
normalizedTensorsList.append(T/inv1)
meanTensorVect = np.mean(normalizedTensorsList, axis=0)
# mean tensor matricial form
meanTensorMtx = tensorvect2matrixform(meanTensorVect)
return meanTensorVect, meanTensorMtx, numTensors
[docs]def covMtx2PPlane(covMtx, meanTensor, numTensors, tensorVectForm=True):
'''Obtains the covariance matrix
:math:`\\boldsymbol{\\mathrm{V^\\mathrm{P}}}` of the
:math:`\\boldsymbol{k^\\mathrm{P}}`'s elements (*i.e.* going over to a
Cartesian system determined by
:math:`\\boldsymbol{p}_1, \\boldsymbol{p}_2\ \&\ \\boldsymbol{p}_3` as is
shown equation (19) of
`Jelínek (1978) <https://doi.org/10.1007/BF01613632>`_), by using the
equation (20) and (21) of the same reference.
Parameters:
covMtx (`numpy.ndarray`): :math:`(6 \\times 6)` array that estimates\
the unbiased covariance matrix :math:`\\boldsymbol{\\mathrm{V}}`\
of the tensors in the sample. It can be obtained by using the\
equation (13) of
`Jelínek (1978) <https://doi.org/10.1007/BF01613632>`_) or the
``cov`` numpy function (as is shown in the example).
meanTensor (`numpy.ndarray`): mean tensor :math:`\\boldsymbol{k}` of\
he sample either in vector or matrix form.
numTensors (`int`): Number of tensors in the sample.
tensorVectForm (`bool`): Logical variable to indicate if the input \
mean tensor is in vector form. ``True`` is the default value.
Returns:
(`numpy.ndarray`): :math:`(6 \\times 6)` covariance matrix \
:math:`\\boldsymbol{\\mathrm{V^\\mathrm{P}}}` of \
:math:`\\boldsymbol{k^\\mathrm{P}}`'s elements.
Examples:
>>> from numpy import cov
>>> from jelinekstat.tools import dataFromFile
>>> from jelinekstat.jelinekstat import meantensor, covMtx2PPlane
>>> sample, numTensors = dataFromFile('inputDataExample.txt')
>>> normTensors = normalizeTensors(sample)
>>> meanTensorVect, meanTensorMtx, numTensors = meantensor(
>>> normTensors, normalized=True)
>>> covMtx = cov(normTensors.T, bias=False)
>>> covMtx2PPlane(
>>> covMtx, meanTensorVect, numTensors)
array([[ 1.02065332e-04, 4.74951816e-05, -1.49560514e-04,
-1.65998575e-06, 3.81988803e-05, -2.24249042e-05],
[ 4.74951816e-05, 5.55684188e-05, -1.03063600e-04,
-2.03320518e-05, -5.50208910e-06, -1.06377394e-05],
[ -1.49560514e-04, -1.03063600e-04, 2.52624114e-04,
2.19920375e-05, -3.26967912e-05, 3.30626436e-05],
[ -1.65998575e-06, -2.03320518e-05, 2.19920375e-05,
2.95042484e-05, 2.05384181e-05, -7.98602106e-06],
[ 3.81988803e-05, -5.50208910e-06, -3.26967912e-05,
2.05384181e-05, 9.00227092e-05, -7.40633382e-05],
[ -2.24249042e-05, -1.06377394e-05, 3.30626436e-05,
-7.98602106e-06, -7.40633382e-05, 9.59108639e-05]])
'''
import numpy as np
from jelinekstat.tools import tensorvect2matrixform, getEigSorted
if tensorVectForm:
meanTensor = tensorvect2matrixform(meanTensor)
eigVal_mT, eigVec_mT = getEigSorted(meanTensor)
t = eigVec_mT
# --- transformation matrix D (equation (21) of Jelinek, 1978) --- #
D = np.array([[t[0, 0]**2, t[1, 0]**2, t[2, 0]**2,
2*t[0, 0]*t[1, 0], 2*t[1, 0]*t[2, 0], 2*t[2, 0]*t[0, 0]],
[t[0, 1]**2, t[1, 1]**2, t[2, 1]**2,
2*t[0, 1]*t[1, 1], 2*t[1, 1]*t[2, 1], 2*t[2, 1]*t[0, 1]],
[t[0, 2]**2, t[1, 2]**2, t[2, 2]**2,
2*t[0, 2]*t[1, 2], 2*t[1, 2]*t[2, 2], 2*t[2, 2]*t[0, 2]],
[t[0, 0]*t[0, 1], t[1, 0]*t[1, 1], t[2, 0]*t[2, 1],
t[0, 0]*t[1, 1]+t[1, 0]*t[0, 1],
t[1, 0]*t[2, 1]+t[2, 0]*t[1, 1],
t[2, 0]*t[0, 1]+t[0, 0]*t[2, 1]],
[t[0, 1]*t[0, 2], t[1, 1]*t[1, 2], t[2, 1]*t[2, 2],
t[0, 1]*t[1, 2]+t[1, 1]*t[0, 2],
t[1, 1]*t[2, 2]+t[2, 1]*t[1, 2],
t[2, 1]*t[0, 2]+t[0, 1]*t[2, 2]],
[t[0, 2]*t[0, 0], t[1, 2]*t[1, 0], t[2, 2]*t[2, 0],
t[0, 2]*t[1, 0]+t[1, 2]*t[0, 0],
t[1, 2]*t[2, 0]+t[2, 2]*t[1, 0],
t[2, 2]*t[0, 0]+t[0, 2]*t[2, 0]]
])
# Here D should be multiplied by (N-1)/N, being N the number of total
# tensor data but the difference is despreciable for big N, ie for N >= 10.
D = (numTensors - 1) / numTensors * D
pCovMtx = np.dot(D, np.dot(covMtx, D.T)) # equation (20) of Jelinek (1978)
return pCovMtx
[docs]def localCovMtxs(meanTensor, pCovMtx, tensorVectForm=True):
'''Determines the covariance matrix :math:`\\boldsymbol{\\mathrm{W}_i}` of
the random variables :math:`\\left(\\mathrm{d}p_{ji},
\\mathrm{d}p_{ki}\\right)` from the local Cartesian System
:math:`\\mathrm{d}\\boldsymbol{p}_i` that define the :math:`\\mathscr{P}-`
plane where each confidence area of the mean tensor's princial direcions
are drawn by using the equation (27) of
`Jelínek (1978) <https://doi.org/10.1007/BF01613632>`_.
Parameters:
meanTensor (`numpy.ndarray`): mean tensor :math:`\\boldsymbol{k}` of
the sample either in vector or matrix form.
pCovMtx (`numpy.ndarray`): :math:`(6 \\times 6)` covariance matrix
:math:`\\boldsymbol{\\mathrm{V^\\mathrm{P}}}` of
:math:`\\boldsymbol{k^\\mathrm{P}}`'s elements. It is obtained
with the ``covMtx2PPlane`` function.
tensorVectForm (`bool`): Logical variable to indicate if the input
:math:`\\boldsymbol{k}` is in vector form. ``True`` is the default
value.
Returns:
Three elements are returned; they are described below.
- **W** (`list`): list of three :math:`(2 \\times 2)` arrayes\
with the covariance matrix :math:`\\boldsymbol{\\mathrm{W}_i}`\
described above.
- **eigVal_W** (`list`): list with the three couples of \
:math:`\\boldsymbol{\\mathrm{W}_i}`'s eigenvalues obtained \
with the ``getEigSorted`` function.
- **eigVec_W** (`list`): list with the three :math:`(2 \\times 2)`\
arrayes of :math:`\\boldsymbol{\\mathrm{W}_i}`'s eigevectors\
obtained with the ``getEigSorted`` function.
Examples:
>>> from numpy import cov
>>> from jelinekstat.tools import dataFromFile
>>> from jelinekstat.jelinekstat import meantensor, covMtx2PPlane
>>> sample, numTensors = dataFromFile('inputDataExample.txt')
>>> normTensors = normalizeTensors(sample)
>>> meanTensorVect, meanTensorMtx, numTensors = meantensor(
>>> normTensors, normalized=True)
>>> covMtx = cov(normTensors.T, bias=False)
>>> pCovMtx = covMtx2PPlane(
>>> covMtx, meanTensorVect, numTensors)
>>> W, eigVal_W, eigVec_W = localCovMtxs(
>>> meanTensorVect, pCovMtx)
>>> W
[array([[ 0.41635005, -0.00798796],
[-0.00798796, 0.00679994]]),
array([[ 0.00739345, -0.02211065],
[-0.02211065, 0.41635005]]),
array([[ 0.00679994, -0.00565157],
[-0.00565157, 0.00739345]])]
>>> eigVal_W
[array([ 0.41650579, 0.0066442 ]),
array([ 0.41754201, 0.00620149]),
array([ 0.01275605, 0.00143733])]
>>> eigVec_W
[array([[ 0.99980999, 0.01949312],
[-0.01949312, 0.99980999]]),
array([[ 0.05383071, -0.99855008],
[-0.99855008, -0.05383071]]),
array([[ 0.68831829, -0.7254088 ],
[-0.7254088 , -0.68831829]])]
'''
import numpy as np
from jelinekstat.tools import tensorvect2matrixform, getEigSorted
if tensorVectForm:
meanTensor = tensorvect2matrixform(meanTensor)
eigVal_mT, eigVec_mT = getEigSorted(meanTensor)
t = eigVal_mT
W = list()
eigVal_W = list()
eigVec_W = list()
for l, r, s in [(0, 1, 2), (1, 2, 0), (2, 0, 1)]:
Wl = np.array(
[[pCovMtx[3 + l, 3 + l] / (t[l] - t[r])**2,
pCovMtx[3 + l, 3 + s] / ((t[l] - t[r]) * (t[l] - t[s]))],
[pCovMtx[3 + l, 3 + s] / ((t[l] - t[r]) * (t[l] - t[s])),
pCovMtx[3 + s, 3 + s] / (t[l] - t[s])**2]])
W.append(Wl)
eigVal, eigVec = getEigSorted(Wl)
eigVal_W.append(eigVal)
eigVec_W.append(eigVec)
return W, eigVal_W, eigVec_W
[docs]def eigValsIntervals(pCovMtx, numTensors, confLvl=0.95, estimate=True):
'''Determines the limits of the variabilities of :math:`\\boldsymbol{k}`'s
principal values for a confidence level given. Ther are obtained by using
the equation (29) of `Jelínek (1978) <https://doi.org/10.1007/BF01613632>`_
or their estimate values by using the equation (35) of
`Jelínek (1978) <https://doi.org/10.1007/BF01613632>`_.
Parameters:
pCovMtx (`numpy.ndarray`): :math:`(6 \\times 6)` covariance matrix
:math:`\\boldsymbol{\\mathrm{V^\\mathrm{P}}}` of
:math:`\\boldsymbol{k^\\mathrm{P}}`'s elements. It is obtained
with the ``covMtx2PPlane`` function.
numTensors (`int`): Number of tensors in the sample.
confLvl (`float`): Confidence level of the limits of the
variabilities of :math:`\\boldsymbol{k}`'s principal values.
``0.95`` is the default value.
estimate (`bool`): Logical variable to indicate if the output is based
whether on the real or estimate covariance matrix
:math:`\\boldsymbol{\\mathrm{V^\\mathrm{P}}}`. ``True`` is the
default value.
Returns:
(`numpy.ndarray`): Array with the three limits of the variabilities of\
:math:`\\boldsymbol{k}`'s principal values.
Examples:
>>> from numpy import cov
>>> from jelinekstat.tools import dataFromFile
>>> from jelinekstat.jelinekstat import *
>>> sample, numTensors = dataFromFile('inputDataExample.txt')
>>> normTensors = normalizeTensors(sample)
>>> meanTensorVect, meanTensorMtx, numTensors = meantensor(
>>> normTensors, normalized=True)
>>> covMtx = cov(normTensors.T, bias=False)
>>> pCovMtx = covMtx2PPlane(
>>> covMtx, meanTensorVect, numTensors)
>>> eigValsIntervals(
>>> pCovMtx, numTensors, confLvl=0.95, estimate=True)
array([ 0.0084461 , 0.00623205, 0.01328784])
'''
import numpy as np
from scipy.stats import norm, t
if estimate: # t-student distribution
stat = abs(t.ppf((1-confLvl)/2, df=numTensors-1)) # Eq. (35)
else: # normal/gaussian distribution
stat = abs(norm.ppf((1-confLvl)/2)) # Eq. (29)
intervals = list()
variance = np.diag(pCovMtx)/numTensors # Eq. (26)
for i in range(3):
intervals.append(stat * variance[i]**0.5)
return np.array(intervals)
[docs]def eigVectsRegions(W, eigVal_W, eigVec_W, numTensors, confLvl=0.95,
estimate=True):
'''Determines the ellipses' geometric parameters of the confidence regions
that define the limits of the variabilities of the
:math:`\\boldsymbol{k}`'s principal vectors.
The axes lenghts ara obtained by using the equation (32) of
`Jelínek (1978) <https://doi.org/10.1007/BF01613632>`_ or their estimated
values by using the equation (35) of
`Jelínek (1978) <https://doi.org/10.1007/BF01613632>`_ and the inclination
angles by using the equation (41) of the same reference.
Parameters:
W (`list`): list of three :math:`(2 \\times 2)` arrayes with the
covariance matrix :math:`\\boldsymbol{\\mathrm{W}_i}` described
above. It is obtained from the ``eigValsIntervals``
function.
eigVal_W (`list`): list with the three couples of
:math:`\\boldsymbol{\\mathrm{W}_i}`'s eigenvalues obtained with
the ``getEigSorted`` function. It is obtained from the
``eigValsIntervals`` function.
eigVec_W (`list`): list with the three :math:`(2 \\times 2)` arrayes
of :math:`\\boldsymbol{\\mathrm{W}_i}`'s eigevectors obtained with
the ``getEigSorted`` function. It is obtained from the
``eigValsIntervals`` function.
numTensors (`int`): Number of tensors in the sample.
confLvl (`float`): Confidence level of the limits of the
variabilities of :math:`\\boldsymbol{k}`'s principal vectors.
``0.95`` is the default value.
estimate (`bool`): Logical variable to indicate if the output is based
whether on the real or estimate covariance matrix
:math:`\\boldsymbol{\\mathrm{V^\\mathrm{P}}}`. ``True`` is the
default value.
Returns:
Three elements are returned; they are described below.
- **majorAxis** (`numpy.ndarray`): Array with the three lengths of\
the ellipses' major axis that define the confidence region. \
The order is acording to the principal values returned from\
the ``getEigSorted`` function.
- **minorAxis** (`list`): list with the three couples of \
:math:`\\boldsymbol{\\mathrm{W}_i}`'s eigenvalues obtained\
with the ``getEigSorted`` function. The order is acording to \
the principal values returned from the ``getEigSorted`` \
function.
- **theta** (`list`): list with the three ellipse inclinations in\
radians measured from the horizontal axis of the local\
Cartesian System of each ellipse to the respective major axis\
counter clockwise.
Examples: ::
>>> from numpy import cov
>>> from jelinekstat.tools import dataFromFile
>>> from jelinekstat.jelinekstat import *
>>> sample, numTensors = dataFromFile('inputDataExample.txt')
>>> normTensors = normalizeTensors(sample)
>>> meanTensorVect, meanTensorMtx, numTensors = meantensor(
>>> normTensors, normalized=True)
>>> covMtx = cov(normTensors.T, bias=False)
>>> pCovMtx = covMtx2PPlane(
>>> covMtx, meanTensorVect, numTensors)
>>> W, eigVal_W, eigVec_W = localCovMtxs(
>>> meanTensorVect, pCovMtx)
>>> majorAxis, minorAxis, theta = eigVectsRegions(
>>> W, eigVal_W, eigVec_W, numTensors, confLvl=0.95,
>>> estimate=True)
>>> majorAxis
array([ 0.66888885, 0.66949335, 0.13745895])
>>> minorAxis
array([ 0.09950548, 0.09615434, 0.04640122])
>>> theta
array([-0.01949436, -1.51693959, -0.81162812])
'''
import numpy as np
from scipy.stats import chi2, f
majorAxis = list()
minorAxis = list()
theta = list()
if estimate:
stat = f.ppf(confLvl, 2, numTensors-2) # F distribution
T2 = stat * 2 * (numTensors-1) / (numTensors * (numTensors-2))
for i in range(3): # Eq. (37)
majorAxis.append(np.arctan((T2 * eigVal_W[i][0])**0.5))
minorAxis.append(np.arctan((T2 * eigVal_W[i][1])**0.5))
theta.append(np.arctan(eigVec_W[i][1, 0]/eigVec_W[i][0, 0]))
else:
stat = chi2.ppf(confLvl, df=2) # Chi-square distribution
for i in range(3): # Eq. (32)
majorAxis.append((stat * eigVal_W[i][0] / numTensors)**0.5)
minorAxis.append((stat * eigVal_W[i][1] / numTensors)**0.5)
theta.append(np.arctan(eigVec_W[i][1, 0]/eigVec_W[i][0, 0]))
return np.array(majorAxis), np.array(minorAxis), np.array(theta)
[docs]def tensorStat(sample, confLevel=0.95, want2plot=True, plotName='001',
ext='pdf'):
'''Summarizes the `Jelínek (1978) <https://doi.org/10.1007/BF01613632>`_
statisctic proposal for 2nd-order tensors and plots it if is wanted.
Parameters:
sample (`numpy.ndarray`): :math:`(n \times 6)` array that
contains the values obtained from the ``extractdata`` fuction.
confLvl (`float`): Confidence level of the limits of the
variabilities of :math:`\boldsymbol{k}`'s principal vectors and
values. ``0.95`` is the default value.
want2plot (`bool`): Logical variable to indicate if is wanted to plot
the summary. ``True`` is the default value.
plotName (`str`): Sample name for saving the final plot. '01' is the
default value.
ext (`str`): File extension for saving the final plot. 'pdf' is the
default value.
Returns:
(`dict`): Summary of the\
`Jelínek (1978) <https://doi.org/10.1007/BF01613632>`_ statisctic\
proposal for 2nd-order tensors where is stored the data related to\
the mean tensor and its variability expressed as the variability\
of their principal values and vectors.
Examples:
>>> from jelinekstat.tools import dataFromFile
>>> from jelinekstat.jelinekstat import tensorStat
>>> sample, numTensors = dataFromFile('inputDataExample.txt')
>>> jelinekStatsSummary, stereonetPlot = tensorStat(
>>> sample, confLevel=0.95, want2plot=True, plotName='example',
>>> ext='pdf')
>>> jelinekStatsSummary
{'K': array([[ 1.02533293, -0.01714126, -0.03794001],
[-0.01714126, 1.01891542, -0.03435001],
[-0.03794001, -0.03435001, 0.95575165]]),
'k': array([1.02533, 1.01891, 0.95575, -0.01714, -0.03435, -0.03794]),
'n': 8,
'k1': {'value': 1.042393712466853, 'variability': 0.008446101031382},
'k2': {'value': 1.033975634080194, 'variability': 0.006232053391831},
'k3': {'value': 0.92363065345295237, 'variability': 0.0132878448184},
'p1': {'coords': array([-0.92438729, 0.196842 , 0.32674357]),
'plg': 19.071242330117354,
'trd': 167.97880567817268,
'majAx': 0.66888885495049721,
'minAx': 0.099505476893979108,
'incl': -1.1169443953398852},
'p2': {'coords': array([-0.04467227, -0.90653975, 0.41975002]),
'plg': 24.818806179014864,
'trd': 267.17887306022931,
'majAx': 0.66949335258868026,
'minAx': 0.096154338416730253,
'incl': -86.914236097035555},
'p3': {'coords': array([-0.37883047, -0.37341521, -0.8467872 ]),
'plg': 57.863927991299327,
'trd': 44.587546531270618,
'majAx': 0.13745894751107776,
'minAx': 0.0464012210689351,
'incl': -46.502865813326189}}
>>> stereonetPlot.show()
.. figure:: https://rawgit.com/eamontoyaa/jelinekstat/master/examples/figures/jelinekstat_tensorStat_example.svg
:alt: jelinekstat_tensorStat_example
.. only:: html
:download:`example script<../examples/figuresScripts/jelinekstat_tensorStat_example.py>`.
'''
import numpy as np
import matplotlib.pyplot as plt
from jelinekstat.tools import getEigSorted, confRegions2PPlanes, \
eigVects2PlgTrd, proyAllEllipses2LongLat, splitIterables
# Normalization of the imput sample
sample = normalizeTensors(sample)
# Mean tensor from normalized data and sample size (n)
k, K, n = meantensor(sample, True) # k (vector); K (matrix)
# Eigenvalues (kK) and eigenvectors (pK) of the mean tensor
kK, pK = getEigSorted(K)
# Unbiased covariance matrix (V).
V = np.cov(sample.T, bias=False)
# Covariance matrix (V) in the system of the k's principal vectors (pV).
pV = covMtx2PPlane(V, k, n)
# Confidence intervals of eigenvalues of mean tensor (kIntervals).
kIntervals = eigValsIntervals(pV, n, confLevel)
# Local covariance matrices (W) in each P-plane of each confidence region.
W, eigValW, eigVectW = localCovMtxs(k, pV)
# Length and orientation of ellipses semi-axis.
majorAxis, minorAxis, theta = eigVectsRegions(
W, eigValW, eigVectW, n, confLevel)
# Coordiantes of the three ellipses in each P-plane.
x, y, PPlanePlots = confRegions2PPlanes(
majorAxis, minorAxis, theta, want2plot, confLevel)
# Stereographic notation to plot the mean tensor's principal vectors (pK).
pKPlg, pKTrd = eigVects2PlgTrd(k) # Plg (plunge); Trd (trend)
# (plunge,trend) notation to plot principal axis of all tensors.
samplePlgTrd = list(map(eigVects2PlgTrd, sample))
# (lon, lat) notation of each confidence region.
kRegionsLong, kRegionsLat = proyAllEllipses2LongLat(x, y, k)
# Summary of the Jelinek (1978) statistic proposal for 2nd-order tensors.
jelinekStatSummary = {
'K': K,
'k': k,
'n': n,
'k1': {'value': kK[0], 'variability': kIntervals[0]},
'k2': {'value': kK[1], 'variability': kIntervals[1]},
'k3': {'value': kK[2], 'variability': kIntervals[2]},
'p1': {'coords': pK[:, 0], 'plg': pKPlg[0], 'trd': pKTrd[0],
'majAx': majorAxis[0], 'minAx': minorAxis[0],
'incl': np.degrees(theta[0])},
'p2': {'coords': pK[:, 1], 'plg': pKPlg[1], 'trd': pKTrd[1],
'majAx': majorAxis[1], 'minAx': minorAxis[1],
'incl': np.degrees(theta[1])},
'p3': {'coords': pK[:, 2], 'plg': pKPlg[2], 'trd': pKTrd[2],
'majAx': majorAxis[2], 'minAx': minorAxis[2],
'incl': np.degrees(theta[2])}
}
# Plotting.
fig = plt.figure(num='Jelinek plot summary')
if want2plot:
plt.ioff()
markers = ['s', '^', 'o']
labels = ['$k_1 = ' + str(round(kK[0], 3)) + '\pm' +
str(round(kIntervals[0], 3)) + '$',
'$k_2 = ' + str(round(kK[1], 3)) + '\pm' +
str(round(kIntervals[1], 3)) + '$',
'$k_3 = ' + str(round(kK[2], 3)) + '\\pm' +
str(round(kIntervals[2], 3)) + '$']
ax = fig.add_subplot(111, projection='stereonet')
# Eigenvectors of all tensors
for tensor in samplePlgTrd:
for i in range(3):
ax.line(tensor[0][i], tensor[1][i], markers[i], color='0.3',
ms=5, fillstyle='none')
# Eigenvectors of mean tensor
for i in range(3):
ax.line(pKPlg[i], pKTrd[i], markers[i], color='k', ms=7,
label=labels[i])
# Confidence regions
for i in range(3):
kRegionsLongSplitted, kRegionsLatSplitted = splitIterables(
kRegionsLong[i], kRegionsLat[i])
for i in range(len(kRegionsLongSplitted)):
ax.plot(kRegionsLongSplitted[i], kRegionsLatSplitted[i], ':k',
lw=1)
# Empty plot to add the confidence region legends.
confLvl = str(round(confLevel * 100, 1))
ax.line(0, 0, ':k', lw=1, label='$'+confLvl + '\%$ conf. regions')
ax.legend(loc=tuple(np.radians([45, -7])), fontsize='x-small')
ax.grid(True, ls='--', lw=0.5)
fig.savefig(plotName + '.' + ext, bbox_inches='tight')
return jelinekStatSummary, fig