glls.hh 5.57 KB
Newer Older
1
#include <stdexcept>  // std::domain_error
2
#include "radixglls/linalg.hh"
3
#include "radixglls/vector.hh"
4

5
6
#include "radixdl/visibility.hh"

7
8
#pragma once

9
namespace radix {
10

11
12
13
14
15
16
struct GllsResult {
  VectorDouble param;      // parameters determined by GLLS [J]
  GllsMatrix param_covar;  // absolute covariances [J][J]
  double chi2;             // the summary statistic for param
  int dof;             // the chi2 degrees of freedom (I - J + rank(econ_lhs))
  double chi2_thresh;  // the upper bound on chi2 for 95% conf interval
17
18
};

19
20
21
22
struct IconEvaluationResult {
  VectorDouble margin;  // distance from ineq constraint [Ri]
  VectorDouble chi2;    // min chi2 along each ineq constraint [Ri]
  double chi2_thresh;   // chi2_thresh with dof=I
23
24
};

25
26
27
class RADIX_PUBLIC GllsConstraintError : public std::logic_error {
 public:
  using std::logic_error::logic_error;
28
29
  GllsConstraintError(const std::string &msg)
      : std::logic_error(msg){};
30
31
};

Aaron Bevill's avatar
Aaron Bevill committed
32
33
/** Given the data and its absolute covariance; calculated data, parameter
 * estimates, and absolute sensitivities; and equality and inequality
34
35
36
37
 * constraints: return the GLLS parameters, their absolute covar, the
 * goodness-of-fit metric chi-squared, the margin from violating each
 * inequality constraint, and the minimum chi-squared along each inequality
 * constraint.
38
 **/
39
40
41
42
43
44
45
46
47
48
49
50
51
52
class RADIX_PUBLIC GllsModel {
 protected:
  // data
  VectorDouble m_data_measured;  // actual data [I]
  GllsMatrix m_data_covar;       // absolute covariances [I][I]
  // model from LM solution
  VectorDouble m_data_modeled;   // data given parameters_lm [I]
  VectorDouble m_param_modeled;  // the LM-calculated parameters [J]
  GllsMatrix m_sensitivities;    // [I][J]
  // constraints
  GllsMatrix m_econ_lhs;    // equality constr [Re][J]
  VectorDouble m_econ_rhs;  // equality constr [Re]
  GllsMatrix m_icon_lhs;    // ineq constr [Ri][J]
  VectorDouble m_icon_rhs;  // ineq constr [Ri]
53

54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
 public:
  GllsModel(const VectorDouble &data_measured, const GllsMatrix &data_covar,
            const VectorDouble &data_modeled, const VectorDouble &param_modeled,
            const GllsMatrix &sensitivities, const GllsMatrix &econ_lhs,
            const VectorDouble &econ_rhs, const GllsMatrix &icon_lhs,
            const VectorDouble &icon_rhs);
  //
  GllsResult solve() const;
  IconEvaluationResult evaluate_icon(const GllsResult &glls_result) const;
  //
  // helpers
  void linearize(GllsMatrix &a,           // linear model lhs [I][J]
                 VectorDouble &b) const;  // linear model rhs [I]
  void apply_econ(
      const GllsMatrix &a, const VectorDouble &b,
      GllsMatrix &acheck,      // reduced linear model lhs [I][J-R]
      VectorDouble &bcheck,    // reduced linear model rhs [I]
      VectorDouble &yi,        // influential transformed parameters
      GllsMatrix &vti,         // influential right singular vectors
      GllsMatrix &vtn) const;  // non-influential right singular vectors
  void system(const GllsMatrix &acheck, const VectorDouble &bcheck,
              GllsMatrix &lhs,           // a^t sigma^-1 a [J][J]
              VectorDouble &rhs) const;  // a^t sigma^-1 b [J]
  double chi2(const VectorDouble &param) const;
  /// The upper limit on chi2 for a 95% confidence interval
  static double chi2_thresh(int dof);
80
};
81

Aaron Bevill's avatar
Aaron Bevill committed
82
/* Fortran interface:
83
 *
Aaron Bevill's avatar
Aaron Bevill committed
84
  interface
85
86
87
    subroutine glls_fwrap( &
        I, &
        J, &
Aaron Bevill's avatar
Aaron Bevill committed
88
89
        Re, &
        Ri, &
90
91
92
93
94
        data_measured, &
        data_covar_cmaj, &
        data_modeled, &
        param_modeled, &
        sensitivities_cmaj, &
Aaron Bevill's avatar
Aaron Bevill committed
95
96
97
98
        econ_lhs_cmaj, &
        econ_rhs, &
        icon_lhs_cmaj, &
        icon_rhs, &
99
        param_glls, &
Aaron Bevill's avatar
Aaron Bevill committed
100
        params_covar_cmaj, &
Aaron Bevill's avatar
Aaron Bevill committed
101
102
103
        chi2, &
        dof, &
        chi2_thresh, &
104
        icon_margin, &
Aaron Bevill's avatar
Aaron Bevill committed
105
106
        icon_chi2, &
        icon_chi2_thresh, &
107
        errcode) bind(c)
108
      import :: c_int, c_double
109
      integer(kind=c_int) :: I, J, Re, Ri, errcode
110
111
112
113
      real(kind=c_double), dimension(I) :: data_measured, data_modeled
      real(kind=c_double), dimension(I*I) :: data_covar_cmaj
      real(kind=c_double), dimension(J) :: param_modeled
      real(kind=c_double), dimension(I*J) :: sensitivities_cmaj
Aaron Bevill's avatar
Aaron Bevill committed
114
115
116
117
      real(kind=c_double), dimension(Re*J) :: econ_lhs_cmaj
      real(kind=c_double), dimension(Re) :: econ_rhs
      real(kind=c_double), dimension(Ri*J) :: icon_lhs_cmaj
      real(kind=c_double), dimension(Ri) :: icon_rhs
118
119
      real(kind=c_double), dimension(J) :: param_glls
      real(kind=c_double), dimension(J*J) :: params_covar_cmaj
Aaron Bevill's avatar
Aaron Bevill committed
120
121
122
      real(kind=c_double) :: chi2, chi2_thresh, icon_chi2_thresh
      integer(kind=c_int) :: dof
      real(kind=c_double), dimension(Ri) :: icon_margin, icon_chi2
123
124
125
126
127
    end subroutine glls_fwrap
  end interface
 *
 * */
extern "C" {
128
129
130
131
132
133
134
135
136
137
138
139
140
void RADIX_PUBLIC glls_fwrap(int *I,   // number of datapoints
                             int *J,   // number of unknown parameters
                             int *Re,  // number of equality constraints
                             int *Ri,  // number of inequality constraints
                             double data_measured[], double data_covar_cmaj[],
                             double data_modeled[], double param_modeled[],
                             double sensitivities_cmaj[],
                             double econ_lhs_cmaj[], double econ_rhs[],
                             double icon_lhs_cmaj[], double icon_rhs[],
                             double param_glls[], double param_covar_cmaj[],
                             double *chi2, int *dof, double *chi2_thresh,
                             double icon_margin[], double icon_chi2[],
                             double *icon_chi2_thresh, int *errcode);
141
142
}

143
}  // namespace radix