-
Notifications
You must be signed in to change notification settings - Fork 3
/
corr_cubic_matrix.cpp
215 lines (193 loc) · 5.81 KB
/
corr_cubic_matrix.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
//#include <Rcpp.h>
#include <RcppArmadillo.h>
using namespace Rcpp;
//using namespace arma;
// This is a simple example of exporting a C++ function to R. You can
// source this function into an R session using the Rcpp::sourceCpp
// function (or via the Source button on the editor toolbar). Learn
// more about Rcpp at:
//
// http://www.rcpp.org/
// http://adv-r.had.co.nz/Rcpp.html
// http://gallery.rcpp.org/
//
// [[Rcpp::export]]
NumericMatrix corr_cubic_matrixC(NumericMatrix x, NumericMatrix y, NumericVector theta) {
int nrow = x.nrow(), ncol = y.nrow();
int nsum = x.ncol();
NumericMatrix out(nrow, ncol);
for (int i = 0; i < nrow; i++) {
for (int j = 0; j < ncol; j++) {
double total = 1;
for(int k = 0; k < nsum; ++k) {
// total += theta[k] * pow((x(i,k) - y(j,k)), 2.0);
double d = fabs(x(i,k) - y(j,k)) / theta[k];
double r = 0;
if (d <= .5) {
r = 1-6*pow(d, 2.0)+6*pow(d, 3.0);
} else if (d <= 1) {
r = 2*pow(1-d, 3.0);
} else {
r = 0;
}
total *= r;
}
out(i, j) = total;
}
}
return out;
}
//' Correlation Cubic matrix in C (symmetric)
//' @param x Matrix x
//' @param theta Theta vector
//' @return Correlation matrix
//' @export
//' @examples
//' corr_cubic_matrix_symC(matrix(c(1,0,0,1),2,2),c(1,1))
// [[Rcpp::export]]
NumericMatrix corr_cubic_matrix_symC(NumericMatrix x, NumericVector theta) {
int nrow = x.nrow();
int nsum = x.ncol();
NumericMatrix out(nrow, nrow);
for (int i = 0; i < nrow - 1; i++) {
for (int j = i + 1; j < nrow; j++) {
double total = 1;
for(int k = 0; k < nsum; ++k) {
// total += theta[k] * pow((x(i,k) - x(j,k)), 2.0);
double d = fabs(x(i,k) - x(j,k)) / theta[k];
double r = 0;
if (d <= .5) {
r = 1-6*pow(d, 2.0) + 6*pow(d, 3.0);
} else if (d <= 1) {
r = 2*pow(1-d, 3.0);
} else {
r = 0;
}
total *= r;
}
out(i, j) = total;
out(j, i) = total; // since symmetric
}
}
for (int i = 0; i < nrow; i++) {
out(i, i) = 1;
}
return out;
}
// [[Rcpp::export]]
NumericVector corr_cubic_matrixvecC(NumericMatrix x, NumericVector y,
NumericVector theta) {
int nrow = x.nrow(); //, ncol = y.nrow();
int nsum = x.ncol();
NumericVector out(nrow);
for (int i = 0; i < nrow; i++) {
double total = 1;
for(int k = 0; k < nsum; ++k) {
// total += theta[k] * pow((x(i,k) - y(k)), 2.0);
double d = fabs(x(i,k) - y(k)) / theta[k];
double r = 0;
if (d <= .5) {
r = 1-6*pow(d, 2.0)+6*pow(d, 3.0);
} else if (d <= 1) {
r = 2*pow(1-d, 3.0);
} else {
r = 0;
}
total *= r;
}
out(i) = total;
}
return out;
}
//' Derivative of cubic kernel covariance matrix in C
//' @param x Matrix x
//' @param theta Theta vector
//' @param C_nonug cov mat without nugget
//' @param s2_est whether s2 is being estimated
//' @param beta_est Whether theta/beta is being estimated
//' @param lenparams_D Number of parameters the derivative is being calculated for
//' @param s2_nug s2 times the nug
//' @param s2 s2
//' @return Correlation matrix
//' @export
// [[Rcpp::export]]
arma::cube kernel_cubic_dC(arma::mat x, arma::vec theta, arma::mat C_nonug,
bool s2_est, bool beta_est, int lenparams_D,
double s2_nug, double s2) {
int nrow = x.n_rows;
int nsum = x.n_cols;
arma::cube dC_dparams(lenparams_D, nrow, nrow);
double log10 = log(10.0);
if (s2_est) {
// dC_dparams(lenparams_D,,) = C * log10;
for (int i = 0; i < nrow - 1; i++) {
for (int j = i + 1; j < nrow; j++) {
dC_dparams(lenparams_D - 1,i,j) = C_nonug(i,j) * log10;
dC_dparams(lenparams_D - 1,j,i) = dC_dparams(lenparams_D - 1,i,j);
}
dC_dparams(lenparams_D - 1, i, i) = (C_nonug(i,i) + s2_nug) * log10;
}
dC_dparams(lenparams_D - 1, nrow - 1, nrow - 1) = (
C_nonug(nrow - 1, nrow - 1) + s2_nug) * log10;
}
if (beta_est) {
for (int i = 0; i < nrow - 1; i++) {
for (int j = i + 1; j < nrow; j++) {
double total = 1;
NumericVector dvec(nsum), rvec(nsum);
for(int k = 0; k < nsum; ++k) {
// total += theta[k] * pow((x(i,k) - x(j,k)), 2.0);
double d = fabs(x(i,k) - x(j,k)) / theta[k];
dvec[k] = d;
double r = 0;
if (d <= .5) {
r = 1-6*pow(d, 2.0)+6*pow(d, 3.0);
} else if (d <= 1) {
r = 2*pow(1-d, 3.0);
} else {
r = 0;
}
rvec[k] = r;
total *= r;
}
for(int k = 0; k < nsum; ++k) {
double grad = 0;
if (x(i,k) - x(j,k) > 0) {
grad = 1;
} else {
grad = -1;
}
double d = fabs(x(i,k) - x(j,k)) / theta[k];
double dr = 0;
if (d <= .5) {
// tmp2 = 1-6*pow(d, 2)+6*pow(d,3);
dr = -12*d+18*pow(d, 2.0);
} else if (d <= 1) {
dr = -6*pow(1-d, 2.0);
} else {
dr = 0;
}
grad *= log10 * (-(x(i,k) - x(j,k))) / theta[k] * dr;
// if (d > 0) {
// grad *= total / rvec[k];
// } else {
for (int l=0; l < nsum; l++) {
if (k != l) {
grad *= rvec[l];
}
}
grad *= s2;
// }
dC_dparams(k,i,j) = grad;
dC_dparams(k,j,i) = dC_dparams(k,i,j);
}
}
}
for (int k=0; k < nsum; k++) {
for (int i = 0; i < nrow; i++) { //# Get diagonal set to zero
dC_dparams(k,i,i) = 0;
}
}
}
return dC_dparams;
}