Determinant Equation W(s)=sI−H(s)¶
The function H(s) is an integral defined as:
It is possible the function H as well as its determinant using the DeterminantEq objects. This is the object used when solving for the approximate missed-events likelihood. The determinant equation is initialized in one of two ways, either from a matrix or QMatrix.
| python: | from dcprogs.likelihood import DeterminantEq, QMatrix
# Define parameters.
matrix = [ [-3050, 50, 3000, 0, 0],
[2./3., -1502./3., 0, 500, 0],
[ 15, 0, -2065, 50, 2000],
[ 0, 15000, 4000, -19000, 0],
[ 0, 0, 10, 0, -10] ]
qmatrix = QMatrix(matrix, 2)
tau = 1e-4
det0 = DeterminantEq(qmatrix, 1e-4);
det1 = DeterminantEq(matrix, 2, 1e-4);
print("{0!s}\n\n{1!s}".format(det0, det1))
|
|---|---|
| c++11: | #include <iostream>
#include <likelihood/determinant_equation.h>
#include <likelihood/root_finder.h>
int main() {
// Define parameters.
DCProgs::t_rmatrix matrix(5 ,5);
matrix << -3050, 50, 3000, 0, 0,
2./3., -1502./3., 0, 500, 0,
15, 0, -2065, 50, 2000,
0, 15000, 4000, -19000, 0,
0, 0, 10, 0, -10;
DCProgs::QMatrix qmatrix(matrix, /*nopen=*/2);
// Create determinant using a QMatrix and a matrix+nopen.
DCProgs::DeterminantEq det0(qmatrix, 1e-4);
DCProgs::DeterminantEq det1(matrix, 2, 1e-4);
std::cout << det0 << "\n\n" << det1 << "\n";
<< " * d[sI-H(s)]/ds for s=-1e2:\n" << DCProgs::numpy_io(det0.s_derivative(-1e2)) << "\n\n";
return 0;
}
|
With an object in hand, it is possible to compute detW(s) for any s. In the following we demonstrate that the two initialization methods are equivalent and that the determinant is zero at the roots of W(s), per definition.
| python: | The python bindings accept both scalars and arrays as input. x = [0, -1, -1e2]
assert all(abs(det0(x) - det1(x)) < 1e-8)
roots = array([-3045.285776037674, -162.92946543451328])
assert all(abs(det0(roots)) < 1e-6 * abs(roots))
|
|---|---|
| c++11: | if( std::abs(det0(0) - det1(0)) > 1e-6
or std::abs(det0(-1) - det1(-1)) > 1e-6
or std::abs(det0(-1e2) - det1(-1e2)) > 1e-6)
throw DCProgs::errors::Runtime("instanciations differ.");
if( std::abs( det0(-3045.285776037674) ) > 1e-6 * 3e3
or std::abs( det0(-162.92946543451328) ) > 1e-6 * 2e2 )
throw DCProgs::errors::Runtime("Roots are not roots.");
|
There exists a convenience function to transform a determinant equation into its “transpose”, e.g. one where A states become F states and F states become A states:
| python: | transpose = det0.transpose()
roots = array( [-17090.192769236815, -2058.0812921673496, -0.24356535498785126],
dtype=internal_dtype )
assert all(abs(transpose(roots)) < 1e-6 * abs(roots))
Note Here we choose to create an input which has same internal type as the dcprogs package. This may result in faster code since no conversion are required. |
|---|---|
| c++11: | DCProgs::DeterminantEq transpose = det0.transpose();
if( std::abs( transpose(-17090.192769236815) ) > 1e-6 * 2e5
or std::abs( transpose(-2058.0812921673496) ) > 1e-6 * 2e3
or std::abs( transpose(-0.24356535498785126) ) > 1e-6 )
throw DCProgs::errors::Runtime("Roots are not roots.");
|
Finally, it is possible to compute H(s) directly, as well as ∂W(s)∂s, as demonstrated below.
| python: | print(" * H({0}):\n{1}\n".format(0, det0.H(0)))
print(" * H({0}):\n{1}\n".format(-1e2, det0.H(-1e2)))
print(" * d[sI-H(s)]/ds for s={0}:\n{1}\n".format(0, det0.s_derivative(0)))
print(" * d[sI-H(s)]/ds for s={0}:\n{1}\n".format(-1e2, det0.s_derivative(-1e2)))
|
|---|---|
| c++11: | std::cout << " * H(0):\n" << DCProgs::numpy_io(det0.H(0)) << "\n\n"
<< " * H(-1e2):\n" << DCProgs::numpy_io(det0.H(-1e2)) << "\n\n";
std::cout << " * d[sI-H(s)]/ds for s=0:\n" << DCProgs::numpy_io(det0.s_derivative(0)) << "\n\n"
<< " * d[sI-H(s)]/ds for s=-1e2:\n" << DCProgs::numpy_io(det0.s_derivative(-1e2)) << "\n\n";
|