Adding radixglls subpackage.

.gitignore

+TriBITS
+googletest
+testframework

CMakeLists.txt

@@ -10,6 +10,8 @@ IF (NOT TRIBITS_PROCESSING_PACKAGE)
   CMAKE_MINIMUM_REQUIRED(VERSION 2.8.11 FATAL_ERROR)
   INCLUDE("${CMAKE_CURRENT_SOURCE_DIR}/ProjectName.cmake")
   PROJECT(${PROJECT_NAME} NONE)
+  # enable C++11 by default
+  SET(radix_ENABLE_CXX11 ON CACHE BOOL "Compile using the C++11 standard" FORCE)
   # disable generating HTML dependencies webpage and xml files
   SET(${PROJECT_NAME}_DEPS_XML_OUTPUT_FILE OFF CACHE BOOL "")
   # disable TriBITS export system to save time configuring

PackagesList.cmake

@@ -12,9 +12,16 @@
 ##---------------------------------------------------------------------------##
 
 TRIBITS_REPOSITORY_DEFINE_PACKAGES(
+  googletest googletest/googletest PT
+  testframework testframework PT
   radix  . PT
 )
 
+TRIBITS_ALLOW_MISSING_EXTERNAL_PACKAGES(
+   googletest
+   testframework
+)
+
 ##---------------------------------------------------------------------------##
 ##                     end of radix/PackagesList.cmake
 ##---------------------------------------------------------------------------##

TPLsList.cmake

@@ -54,6 +54,7 @@
 #
 
 SET(radix_TPLS_FINDMODS_CLASSIFICATIONS
+ LAPACK "TriBITS/tribits/common_tpls/FindTPLLAPACK.cmake" PT
 )
 
 ##---------------------------------------------------------------------------##

cmake/CallbackSetupExtraOptions.cmake

@@ -12,6 +12,10 @@ ENDIF()
 ##---------------------------------------------------------------------------##
 
 MACRO(TRIBITS_REPOSITORY_SETUP_EXTRA_OPTIONS)
+  #
+  # Include the testframework setup
+  #
+  ADD_SUBDIRECTORY(${radix_SOURCE_DIR}/testframework/setup)
 
   # Set CXX11 to be enabled by default.
   SET(${PROJECT_NAME}_ENABLE_CXX11_DEFAULT TRUE)

cmake/Dependencies.cmake

@@ -2,6 +2,7 @@ TRIBITS_PACKAGE_DEFINE_DEPENDENCIES(
    SUBPACKAGES_DIRS_CLASSIFICATIONS_OPTREQS
      bug  radixbug      SS OPTIONAL
      math radixmath     SS OPTIONAL
+     glls radixglls     SS OPTIONAL
      io   radixio       SS OPTIONAL
      geometry radixgeometry SS OPTIONAL
      plot radixplot     SS OPTIONAL

cmake/ExtraRepositoriesList.cmake

+TRIBITS_PROJECT_DEFINE_EXTRA_REPOSITORIES(
+  GTEST_REPO  "googletest" GIT https://github.com/lefebvre/googletest "NOPACKAGES" Continuous
+  TESTF_REPO  "testframework" GIT https://github.com/lefebvre/testframework "NOPACKAGES" Continuous
+)

radixglls/CMakeLists.txt

+TRIBITS_SUBPACKAGE(glls)
+
+SET(SOURCE
+    vector.cc
+    exception.cc
+    linalg.cc
+    glls.cc
+)
+
+TRIBITS_ADD_LIBRARY(radixglls
+  SOURCES  ${SOURCE}
+)
+TRIBITS_ADD_TEST_DIRECTORIES(tests)
+
+TRIBITS_SUBPACKAGE_POSTPROCESS()
+

radixglls/cmake/Dependencies.cmake

+TRIBITS_PACKAGE_DEFINE_DEPENDENCIES(
+LIB_REQUIRED_PACKAGES radixbug
+LIB_OPTIONAL_PACKAGES
+TEST_REQUIRED_PACKAGES testframework
+TEST_OPTIONAL_PACKAGES
+LIB_REQUIRED_TPLS LAPACK
+LIB_OPTIONAL_TPLS
+TEST_REQUIRED_TPLS
+TEST_OPTIONAL_TPLS
+)

radixglls/exception.cc

+#include "radixglls/exception.hh"
+#include <cstring>  // std::strcpy, std::strcat, std::strlen
+
+namespace radix {
+
+void runtime_assert(bool assertion, const char * format, ...)
+{
+    if (not assertion)
+    {
+        // add a newline and a runtime_assert tag
+        char * new_format = new char[16 + std::strlen(format) + 1];
+        std::strcpy(new_format, "runtime_assert: ");
+        std::strcpy(new_format, format);
+        std::strcat(new_format, "\n");
+        // create a formatted string and an exception object
+        va_list arglist;
+        va_start(arglist, format); // get the args after format
+        size_t formatted_size = std::strlen(format) + 200;
+        char * formatted = new char[formatted_size];
+        vsnprintf(formatted, formatted_size, new_format, arglist);
+        std::runtime_error exception(formatted);
+        va_end(arglist);
+        // destruct and throw
+        delete [] new_format;
+        delete [] formatted;
+        throw exception;
+    }
+}
+
+}

radixglls/exception.hh

+#include <stdarg.h> // variadic parameters for debug_assert()
+#include <stdexcept> // std::exception, std::runtime_error
+
+#pragma once
+
+namespace radix {
+
+void runtime_assert(bool assertion, const char * format, ...);
+
+}

radixglls/glls.cc

+#include <numeric>  // std::inner_product
+#include <cmath>  // std::sqrt
+#include "radixglls/glls.hh"
+#include "radixglls/exception.hh"
+
+namespace radix {
+
+void glls(
+    const VectorDouble & data_measured,
+    const VectorVectorDouble & data_covar,
+    const VectorDouble & data_modeled,
+    const VectorDouble & param_modeled,
+    const VectorVectorDouble & sensitivities,
+    VectorDouble & param_glls,
+    VectorVectorDouble & param_covar)
+{
+    VectorVectorDouble a;
+    VectorDouble b;
+    glls_linearize(
+        data_measured,
+        data_covar,
+        data_modeled,
+        param_modeled,
+        sensitivities,
+        a,
+        b);
+    //
+    VectorVectorDouble lhs;
+    VectorDouble rhs;
+    glls_system(a, b, data_covar, lhs, rhs);
+    //
+    glls_solve(lhs, rhs, param_glls, param_covar);
+}
+
+void glls(
+    const VectorDouble & data_measured,
+    const VectorVectorDouble & data_covar,
+    const VectorDouble & data_modeled,
+    const VectorDouble & param_modeled,
+    const VectorVectorDouble & sensitivities,
+    VectorDouble & param_glls,
+    VectorVectorDouble & param_covar,
+    const VectorVectorDouble & constraints_lhs,
+    const VectorDouble & constraints_rhs)
+{
+    if (constraints_lhs.size() == 0)
+    {
+        glls(
+            data_measured,
+            data_covar,
+            data_modeled,
+            param_modeled,
+            sensitivities,
+            param_glls,
+            param_covar);
+        return;
+    }
+    VectorVectorDouble a;
+    VectorDouble b;
+    glls_linearize(
+        data_measured,
+        data_covar,
+        data_modeled,
+        param_modeled,
+        sensitivities,
+        a,
+        b);
+    //
+    VectorVectorDouble acheck, vti, vtn;
+    VectorDouble bcheck, yi;
+    apply_constraints(
+        a,
+        b,
+        constraints_lhs,
+        constraints_rhs,
+        acheck,
+        bcheck,
+        yi,
+        vti,
+        vtn);
+    //
+    VectorVectorDouble lhs;
+    VectorDouble rhs;
+    glls_system(acheck, bcheck, data_covar, lhs, rhs);
+    //
+    VectorDouble yn;
+    VectorVectorDouble covar_yn;
+    glls_solve(lhs, rhs, yn, covar_yn);
+    //
+    // param_glls = vti.T.dot(yi) + vtn.T.dot(yn)
+    param_glls.resize(a[0].size());
+    VectorVectorDouble vi = transpose(vti), vn = transpose(vtn);
+    VectorDouble vi_yi = multiply(vi, yi);
+    VectorDouble vn_yn = multiply(vn, yn);
+    for (size_t j = 0; j < param_glls.size(); ++j)
+    {
+        param_glls[j] = vi_yi[j] + vn_yn[j];
+    }
+    // param_covar = vtn.T.dot(covar_yn).dot(vtn)
+    param_covar = multiply(vn, covar_yn);
+    param_covar = multiply(param_covar, vtn);
+}
+
+void glls_linearize(
+    const VectorDouble & data_measured,
+    const VectorVectorDouble & data_covar,
+    const VectorDouble & data_modeled,
+    const VectorDouble & param_modeled,
+    const VectorVectorDouble & sensitivities,
+    VectorVectorDouble & a,
+    VectorDouble & b)
+{
+    const size_t I = sensitivities.size();  // num measurements
+    const size_t J = param_modeled.size();  // num parameters
+    //
+    // check argument consistency
+    for (size_t i = 0; i < I; ++i)
+    {
+        runtime_assert(
+            sensitivities[i].size() == J,
+            "sensitivities[i].size() == J");
+    }
+    runtime_assert(
+        data_modeled.size() == I,
+        "data_modeled.size() == I");
+    runtime_assert(
+        data_measured.size() == I,
+        "data_measured.size() == I");
+    runtime_assert(
+        data_covar.size() == I,
+        "data_covar.size() == I");
+    for (size_t i = 0; i < I; ++i)
+    {
+        runtime_assert(
+            data_covar[i].size() == I,
+            "data_covar[i].size() == I");
+    }
+    //
+    // build the residual model
+    a = sensitivities;
+    b = multiply(sensitivities, param_modeled);
+    for (size_t i = 0; i < I; ++i)
+    {
+        b[i] += data_measured[i] - data_modeled[i];
+    }
+}
+
+void apply_constraints(
+    const VectorVectorDouble & a,
+    const VectorDouble & b,
+    const VectorVectorDouble & constraints_lhs,
+    const VectorDouble & constraints_rhs,
+    VectorVectorDouble & acheck,
+    VectorDouble & bcheck,
+    VectorDouble & yi,
+    VectorVectorDouble & vti,
+    VectorVectorDouble & vtn)
+{
+    runtime_assert(
+        constraints_lhs.size() <= a[0].size(),
+        "apply_constraints must have fewer constraints than parameters");
+    runtime_assert(
+        constraints_lhs[0].size() == a[0].size(),
+        "apply_constraints constraints_lhs--a shape mismatch");
+    for (size_t i = 0; i < constraints_lhs.size(); ++i)
+    {
+        runtime_assert(
+            constraints_lhs[i].size() == constraints_lhs[0].size(),
+            "apply_constraints constraints_lhs must be square");
+    }
+    runtime_assert(
+        constraints_lhs.size() == constraints_rhs.size(),
+        "apply_constraints constraints_lhs--constraints_rhs shape mismatch");
+    //
+    VectorVectorDouble u;
+    VectorDouble svalues;
+    svd(constraints_lhs, u, svalues, vti);
+    //
+    // determine number of influential vectors
+    size_t rank = 0;
+    for (size_t i = 0; i < svalues.size(); ++i)
+    {
+        if (svalues[i] < 1e-10)
+        {
+            break;
+        }
+        ++rank;
+    }
+    //
+    // split into influential and non-influential vectors (rows)
+    vtn.assign(vti.begin() + rank, vti.end());
+    vti.resize(rank);
+    //
+    // create the transformed problem and solve for the influential parameters
+    // acheck = a.dot(vtn.T)
+    acheck = multiply(a, transpose(vtn));
+    // yi = solve(U diag(s), constraints_rhs) = diag(1 / s) U^t constraints_rhs
+    yi = multiply(transpose(u), constraints_rhs);
+    for (size_t i = 0; i < yi.size(); ++i)
+    {
+            yi[i] /= svalues[i];
+    }
+    // bcheck = b - a.dot(vti.T).dot(yi)
+    VectorVectorDouble a_vi = multiply(a, transpose(vti));
+    bcheck = multiply(a_vi, yi);
+    for (size_t i = 0; i < bcheck.size(); ++i)
+    {
+        bcheck[i] = b[i] - bcheck[i];
+    }
+}
+
+void glls_system(
+    const VectorVectorDouble & a,
+    const VectorDouble & b,
+    const VectorVectorDouble & data_covar,
+    VectorVectorDouble & lhs,
+    VectorDouble & rhs)
+{
+    size_t I = a.size(), J = a[0].size();
+    //
+    // check whether the datapoints are correlated
+    bool correlated = false;
+    for (size_t i1 = 0; i1 < I; ++i1)
+    {
+        for (size_t i2 = 0; i2 < I; ++i2)
+        {
+            if (data_covar[i1][i2] != 0. and i1 != i2)
+            {
+                correlated = true;
+            }
+        }
+    }
+    //
+    lhs.assign(J, VectorDouble(J, 0.));
+    rhs.assign(J, 0.);
+    if (correlated)
+    {
+        FactorizedMatrix covar_fm(data_covar);
+        // lhs = a^t covar^-1 a
+        VectorVectorDouble covarinva(a.begin(), a.end());
+        covar_fm.solve(covarinva);
+        for (size_t jc = 0; jc < J; ++jc)
+        {
+            for (size_t jr = 0; jr < J; ++jr)
+            {
+                for (size_t i = 0; i < I; ++i)
+                {
+                    lhs[jr][jc] += a[i][jr] * covarinva[i][jc];
+                }
+            }
+        }
+        // rhs = a^t covar^-1 b
+        VectorDouble covarinvb(b.begin(), b.end());
+        covar_fm.solve(covarinvb);
+        for (size_t jr = 0; jr < J; ++jr)
+        {
+            for (size_t i = 0; i < I; ++i)
+            {
+                rhs[jr] += a[i][jr] * covarinvb[i];
+            }
+        }
+    }
+    else
+    {
+        // lhs = a^t covar^-1 a
+        for (size_t jc = 0; jc < J; ++jc)
+        {
+            for (size_t jr = 0; jr < J; ++jr)
+            {
+                for (size_t i = 0; i < I; ++i)
+                {
+                    lhs[jr][jc] += a[i][jr] / data_covar[i][i] * a[i][jc];
+                }
+            }
+        }
+        // rhs = a^t covar^-1 b
+        for (size_t jr = 0; jr < J; ++jr)
+        {
+            for (size_t i = 0; i < I; ++i)
+            {
+                rhs[jr] += a[i][jr] / data_covar[i][i] * b[i];
+            }
+        }
+    }
+}
+
+void glls_solve(
+    const VectorVectorDouble & lhs,
+    const VectorDouble & rhs,
+    VectorDouble & param_glls,
+    VectorVectorDouble & param_covar)
+{
+    FactorizedMatrix factorized_lhs(lhs);
+    param_glls = rhs;
+    factorized_lhs.solve(param_glls);
+    param_covar = factorized_lhs.inverse();
+}
+
+}
+
+
+

radixglls/glls.hh

+#include "radixglls/vector.hh"
+#include "radixglls/linalg.hh"
+
+#pragma once
+
+namespace radix {
+
+/** Given the data, parameter estimates, absolute sensitivities, data
+ * modeled at the parameter estimates, and data absolute covar, return
+ * the GLLS parameters and their absolute covar.
+ **/
+void glls(
+    const VectorDouble & data_measured,  // actual data [I]
+    const VectorVectorDouble & data_covar,  // absolute covariances [I][I]
+    const VectorDouble & data_modeled,  // data given parameters_lm [I]
+    const VectorDouble & param_modeled,  // the LM-calculated parameters [J]
+    const VectorVectorDouble & sensitivities,  // [I][J]
+    VectorDouble & param_glls,  // parameters determined by GLLS [J]
+    VectorVectorDouble & param_covar);  // absolute covariances [J][J]
+void glls(
+    const VectorDouble & data_measured,  // actual data [I]
+    const VectorVectorDouble & data_covar,  // absolute covariances [I][I]
+    const VectorDouble & data_modeled,  // data given parameters_lm [I]
+    const VectorDouble & param_modeled,  // the LM-calculated parameters [J]
+    const VectorVectorDouble & sensitivities,  // [I][J]
+    VectorDouble & param_glls,  // parameters determined by GLLS [J]
+    VectorVectorDouble & param_covar,  // absolute covariances [J][J]
+    const VectorVectorDouble & constraints_lhs,  // equality constr [R][J]
+    const VectorDouble & constraints_rhs);  // equality constr [J]
+
+/** Helpers for glls_solve **/
+
+void glls_linearize(
+    const VectorDouble & data_measured,
+    const VectorVectorDouble & data_covar,
+    const VectorDouble & data_modeled,
+    const VectorDouble & param_modeled,
+    const VectorVectorDouble & sensitivities,
+    VectorVectorDouble & a,  // linear model lhs [I][J]
+    VectorDouble & b);  // linear model rhs [I]
+
+void apply_constraints(
+    const VectorVectorDouble & a,
+    const VectorDouble & b,
+    const VectorVectorDouble & constraints_lhs,
+    const VectorDouble & constraints_rhs,
+    VectorVectorDouble & acheck,  // reduced linear model lhs [I][J-R]
+    VectorDouble & bcheck,  // reduced linear model rhs [I]
+    VectorDouble & yi,  // influential transformed parameters
+    VectorVectorDouble & vti,  // influential right singular vectors
+    VectorVectorDouble & vtn);  // non-influential right singular vectors
+
+void glls_system(
+    const VectorVectorDouble & a,
+    const VectorDouble & b,
+    const VectorVectorDouble & data_covar,
+    VectorVectorDouble & lhs,  // a^t sigma^-1 a [J][J]
+    VectorDouble & rhs);  // a^t sigma^-1 b [J]
+
+void glls_solve(
+    const VectorVectorDouble & lhs,
+    const VectorDouble & rhs,
+    VectorDouble & param_glls,
+    VectorVectorDouble & param_covar);
+
+}

radixglls/linalg.cc

+#include <numeric>  // std::inner_product
+#include <cmath>  // std::sqrt
+#include "radixglls/linalg.hh"
+#include "radixglls/exception.hh"
+
+namespace radix {
+
+VectorVectorDouble diag(const VectorDouble & vec)
+{
+    size_t n = vec.size();
+    VectorVectorDouble diag(n, VectorDouble(n, 0.));
+    for (size_t i = 0; i < n; ++i)
+    {
+        diag[i][i] = vec[i];
+    }
+    return diag;
+}
+
+VectorVectorDouble diag_sq(const VectorDouble & vec)
+{
+    size_t n = vec.size();
+    VectorVectorDouble diag(n, VectorDouble(n, 0.));
+    for (size_t i = 0; i < n; ++i)
+    {
+        diag[i][i] = vec[i] * vec[i];
+    }
+    return diag;
+}
+
+VectorVectorDouble diag_inv(const VectorDouble & vec)
+{
+    size_t n = vec.size();
+    VectorVectorDouble diag(n, VectorDouble(n, 0.));
+    for (size_t i = 0; i < n; ++i)
+    {
+        diag[i][i] = 1. / vec[i];
+    }
+    return diag;
+}
+
+VectorVectorDouble transpose(const VectorVectorDouble & matrix)
+{
+    VectorVectorDouble tr(matrix[0].size(), VectorDouble(matrix.size()));
+    for (size_t j = 0; j < matrix[0].size(); ++j)
+    {
+        for (size_t i = 0; i < matrix.size(); ++i)
+        {
+            tr[j][i] = matrix[i][j];
+        }
+    }
+    return tr;
+}
+
+VectorDouble multiply(const VectorVectorDouble & A, const VectorDouble & b)
+{
+    VectorDouble prod(A.size());
+    for (size_t i = 0; i < prod.size(); ++i)
+    {
+        prod[i] = std::inner_product(
+            A[i].begin(),
+            A[i].end(),
+            b.begin(),
+            0.);
+    }
+    return prod;
+}
+
+VectorVectorDouble multiply(
+    const VectorVectorDouble & A,
+    const VectorVectorDouble & B)
+{
+    VectorVectorDouble Btrans = transpose(B);
+    VectorVectorDouble prod(A.size(), VectorDouble(Btrans.size(), 0.));
+    for (size_t i = 0; i < A.size(); ++i)
+    {
+        for (size_t k = 0; k < B[0].size(); ++k)
+        {
+            prod[i][k] = std::inner_product(
+                A[i].begin(),
+                A[i].end(),
+                Btrans[k].begin(),
+                0.);
+        }
+    }
+    return prod;
+}
+
+// LAPACK dsytrf factorizes a double-precision symmetric matrix
+extern "C" void dsytrf_(
+    char * uplo,  // 'U' or 'L' for shape of triangular decomposition
+    int * n,  // a is n by n
+    double * a,  // lhs matrix, serialized column-major (array inout)
+    int * lda,  // leading dimension of A (for slicing into A)
+    int * ipiv,  // pivot array (array out)
+    double * work,  // work(1) is the optimal lwork (array out)
+    int * lwork,  // size of the work array, suggest >= n
+    int * info);  // 0 indicates successful exit (out)
+
+// LAPACK dsytri inverts a factorized matrix
+extern "C" void dsytri_(
+    char * uplo,  // 'U' or 'L' for shape of triangular decomposition
+    int * n,  // a is n by n
+    double * a,  // lhs matrix, serialized column-major (array inout)
+    int * lda,  // leading dimension of A (for slicing into A)
+    int * ipiv,  // pivot array (array out)
+    double * work,  // work(1) is the optimal lwork (array out)
+    int * info);  // 0 indicates successful exit (out)
+
+// LAPACK dsytrs solves a factorized matrix
+extern "C" void dsytrs_(