Skip to content

Commit 1426f54

Browse files
oddish3oddish3
committed
new branch
1 parent 65ca9d5 commit 1426f54

15 files changed

+167
-112
lines changed

‎.github/workflows/R-CMD-check.yaml

+7
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,13 @@ jobs:
4444
message("load_all completed")
4545
'
4646
shell: bash
47+
- name: Enable core dumps (Windows)
48+
if: runner.os == 'Windows'
49+
run: |
50+
reg add "HKLM\SOFTWARE\Microsoft\Windows\Windows Error Reporting\LocalDumps" /v DumpFolder /t REG_EXPAND_SZ /d "%LOCALAPPDATA%\CrashDumps" /f
51+
reg add "HKLM\SOFTWARE\Microsoft\Windows\Windows Error Reporting\LocalDumps" /v DumpCount /t REG_DWORD /d 10 /f
52+
reg add "HKLM\SOFTWARE\Microsoft\Windows\Windows Error Reporting\LocalDumps" /v DumpType /t REG_DWORD /d 2 /f
53+
shell: cmd
4754
- name: Run tests
4855
run: |
4956
Rscript -e '

‎DESCRIPTION

+6-2
Original file line numberDiff line numberDiff line change
@@ -18,13 +18,17 @@ Imports:
1818
devtools,
1919
fixest (>= 0.8.0),
2020
MASS (>= 7.3-53),
21-
progress (>= 1.2.2)
21+
progress (>= 1.2.2),
22+
2223
LinkingTo:
2324
Rcpp,
2425
RcppArmadillo
2526
RoxygenNote: 7.3.2
2627
Suggests:
2728
testthat (>= 3.0.0),
28-
tidyverse
29+
tidyverse,
30+
dplyr,
31+
purrr,
32+
readr
2933
Config/testthat/edition: 3
3034
SystemRequirements: C++11

‎src/RcppExports.cpp

+6-6
Original file line numberDiff line numberDiff line change
@@ -12,20 +12,20 @@ Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
1212
#endif
1313

1414
// bspline
15-
List bspline(NumericVector x, int l, int r);
15+
List bspline(NumericVector x, size_t l, size_t r);
1616
RcppExport SEXP _contdid_bspline(SEXP xSEXP, SEXP lSEXP, SEXP rSEXP) {
1717
BEGIN_RCPP
1818
Rcpp::RObject rcpp_result_gen;
1919
Rcpp::RNGScope rcpp_rngScope_gen;
2020
Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP);
21-
Rcpp::traits::input_parameter< int >::type l(lSEXP);
22-
Rcpp::traits::input_parameter< int >::type r(rSEXP);
21+
Rcpp::traits::input_parameter< size_t >::type l(lSEXP);
22+
Rcpp::traits::input_parameter< size_t >::type r(rSEXP);
2323
rcpp_result_gen = Rcpp::wrap(bspline(x, l, r));
2424
return rcpp_result_gen;
2525
END_RCPP
2626
}
2727
// jhat
28-
Rcpp::List jhat(const arma::mat& PP, const arma::mat& BB, const arma::vec& CJ, const arma::vec& CK, const arma::vec& TJ, double M, int n, int nL);
28+
Rcpp::List jhat(const arma::mat& PP, const arma::mat& BB, const arma::vec& CJ, const arma::vec& CK, const arma::vec& TJ, double M, arma::uword n, arma::uword nL);
2929
RcppExport SEXP _contdid_jhat(SEXP PPSEXP, SEXP BBSEXP, SEXP CJSEXP, SEXP CKSEXP, SEXP TJSEXP, SEXP MSEXP, SEXP nSEXP, SEXP nLSEXP) {
3030
BEGIN_RCPP
3131
Rcpp::RObject rcpp_result_gen;
@@ -36,8 +36,8 @@ BEGIN_RCPP
3636
Rcpp::traits::input_parameter< const arma::vec& >::type CK(CKSEXP);
3737
Rcpp::traits::input_parameter< const arma::vec& >::type TJ(TJSEXP);
3838
Rcpp::traits::input_parameter< double >::type M(MSEXP);
39-
Rcpp::traits::input_parameter< int >::type n(nSEXP);
40-
Rcpp::traits::input_parameter< int >::type nL(nLSEXP);
39+
Rcpp::traits::input_parameter< arma::uword >::type n(nSEXP);
40+
Rcpp::traits::input_parameter< arma::uword >::type nL(nLSEXP);
4141
rcpp_result_gen = Rcpp::wrap(jhat(PP, BB, CJ, CK, TJ, M, n, nL));
4242
return rcpp_result_gen;
4343
END_RCPP

‎src/RcppExports.o

11.8 KB
Binary file not shown.

‎src/bspline.cpp

+25-16
Original file line numberDiff line numberDiff line change
@@ -1,77 +1,86 @@
11
#include <RcppArmadillo.h>
22
// [[Rcpp::depends(RcppArmadillo)]]
3-
43
using namespace Rcpp;
54
using namespace arma;
65

76
// [[Rcpp::export]]
8-
List bspline(NumericVector x, int l, int r) {
9-
int N = x.size();
10-
int m = std::pow(2, l) - 1;
7+
List bspline(NumericVector x, size_t l, size_t r) {
8+
Rcpp::Rcout << "Entering bspline function" << std::endl;
9+
Rcpp::Rcout << "Input parameters: l = " << l << ", r = " << r << std::endl;
10+
11+
size_t N = x.size();
12+
size_t m = std::pow(2, l) - 1;
1113
r = r + 1;
1214

15+
Rcpp::Rcout << "N = " << N << ", m = " << m << ", r = " << r << std::endl;
16+
1317
// Define the augmented knot set
1418
std::vector<double> kts;
1519
if (l == 0) {
1620
kts.resize(2 * (r - 1), 0.0);
1721
std::fill(kts.begin() + (r - 1), kts.end(), 1.0);
1822
} else if (l >= 1) {
1923
kts.resize(2 * (r - 2) + m + 2, 0.0);
20-
for (int i = 0; i <= m; ++i) {
24+
for (size_t i = 0; i <= m; ++i) {
2125
kts[r - 2 + i] = static_cast<double>(i) / std::pow(2, l);
2226
}
2327
std::fill(kts.begin() + r - 2 + m + 1, kts.end(), 1.0);
2428
}
2529

30+
Rcpp::Rcout << "Knot set size: " << kts.size() << std::endl;
31+
2632
// Initialize for recursion
2733
arma::cube BB(N, m + 2 * r - 2, r - 1, arma::fill::zeros);
28-
for (int i = 0; i < N; ++i) {
29-
for (int j = r - 1; j <= r + m - 1; ++j) {
34+
for (size_t i = 0; i < N; ++i) {
35+
for (size_t j = r - 1; j <= r + m - 1; ++j) {
3036
if (x[i] >= kts[j] && x[i] < kts[j + 1]) {
3137
BB(i, j - (r - 1), 0) = 1.0;
3238
break;
3339
}
3440
}
3541
}
3642

43+
Rcpp::Rcout << "BB cube dimensions: " << BB.n_rows << "x" << BB.n_cols << "x" << BB.n_slices << std::endl;
44+
3745
// Recursion
38-
for (int j = 1; j < r - 1; ++j) {
39-
for (int i = 0; i < m + 2 * r - 2 - j; ++i) {
46+
for (size_t j = 1; j < r - 1; ++j) {
47+
for (size_t i = 0; i < m + 2 * r - 2 - j; ++i) {
4048
arma::vec a1(N, arma::fill::zeros);
4149
if (kts[i + j] - kts[i] != 0) {
4250
a1 = (x - kts[i]) / (kts[i + j] - kts[i]);
4351
}
44-
4552
arma::vec a2(N, arma::fill::zeros);
4653
if (kts[i + j + 1] - kts[i + 1] != 0) {
47-
a2 = (x - kts[i + 1]) / (kts[i + j + 1] - kts[i + 1]);
54+
a2 = (kts[i + j + 1] - x) / (kts[i + j + 1] - kts[i + 1]);
4855
}
49-
50-
BB.slice(j).col(i) = a1 % BB.slice(j-1).col(i) + (1 - a2) % BB.slice(j-1).col(i+1);
56+
BB.slice(j).col(i) = a1 % BB.slice(j-1).col(i) + a2 % BB.slice(j-1).col(i+1);
5157
}
5258
}
5359

5460
arma::mat XX = BB.slice(r - 2).cols(0, std::pow(2, l) + r - 3);
5561

62+
Rcpp::Rcout << "XX matrix dimensions: " << XX.n_rows << "x" << XX.n_cols << std::endl;
63+
5664
// Calculate derivative
5765
arma::mat DX(N, m + r - 1, arma::fill::zeros);
58-
for (int i = 0; i < m + r - 1; ++i) {
66+
for (size_t i = 0; i < m + r - 1; ++i) {
5967
arma::vec a1(N, arma::fill::zeros);
6068
if (kts[i + r - 2] - kts[i] != 0) {
6169
a1.fill(1.0 / (kts[i + r - 2] - kts[i]));
6270
}
63-
6471
arma::vec a2(N, arma::fill::zeros);
6572
if (kts[i + r - 1] - kts[i + 1] != 0) {
6673
a2.fill(1.0 / (kts[i + r - 1] - kts[i + 1]));
6774
}
68-
6975
if (i < m + r - 2) {
7076
DX.col(i) = (r - 2) * (a1 % BB.slice(r - 3).col(i) - a2 % BB.slice(r - 3).col(i + 1));
7177
} else {
7278
DX.col(i) = (r - 2) * (a1 % BB.slice(r - 3).col(i));
7379
}
7480
}
7581

82+
Rcpp::Rcout << "DX matrix dimensions: " << DX.n_rows << "x" << DX.n_cols << std::endl;
83+
Rcpp::Rcout << "Exiting bspline function" << std::endl;
84+
7685
return List::create(Named("XX") = XX, Named("DX") = DX);
7786
}

‎src/bspline.o

21.4 KB
Binary file not shown.

‎src/contdid.so

71.7 KB
Binary file not shown.

‎src/jhat.cpp

+29-13
Original file line numberDiff line numberDiff line change
@@ -1,16 +1,15 @@
11
#include <RcppArmadillo.h>
22
#include "shat.h"
33
// [[Rcpp::depends(RcppArmadillo)]]
4-
54
using namespace Rcpp;
65
using namespace arma;
76

87
// Helper function to print matrix summary
98
void print_matrix_summary(const arma::mat& matrix, const std::string& name) {
10-
// Rcpp::Rcout << "Matrix " << name << " dimensions: " << matrix.n_rows << "x" << matrix.n_cols << std::endl;
11-
// Rcpp::Rcout << "First few elements of " << name << ":" << std::endl;
12-
for (uword i = 0; i < std::min(10u, matrix.n_rows); ++i) {
13-
for (uword j = 0; j < std::min(5u, matrix.n_cols); ++j) {
9+
Rcpp::Rcout << "Matrix " << name << " dimensions: " << matrix.n_rows << "x" << matrix.n_cols << std::endl;
10+
Rcpp::Rcout << "First few elements of " << name << ":" << std::endl;
11+
for (arma::uword i = 0; i < std::min(static_cast<arma::uword>(10), matrix.n_rows); ++i) {
12+
for (arma::uword j = 0; j < std::min(static_cast<arma::uword>(5), matrix.n_cols); ++j) {
1413
Rcpp::Rcout << std::setprecision(10) << std::setw(15) << matrix(i, j) << " ";
1514
}
1615
Rcpp::Rcout << std::endl;
@@ -21,18 +20,32 @@ void print_matrix_summary(const arma::mat& matrix, const std::string& name) {
2120
// [[Rcpp::export]]
2221
Rcpp::List jhat(const arma::mat& PP, const arma::mat& BB,
2322
const arma::vec& CJ, const arma::vec& CK,
24-
const arma::vec& TJ, double M, int n, int nL) {
23+
const arma::vec& TJ, double M, arma::uword n, arma::uword nL) {
24+
Rcpp::Rcout << "Entering jhat function" << std::endl;
25+
Rcpp::Rcout << "Input dimensions: PP(" << PP.n_rows << "x" << PP.n_cols << "), BB(" << BB.n_rows << "x" << BB.n_cols << ")" << std::endl;
26+
Rcpp::Rcout << "CJ size: " << CJ.n_elem << ", CK size: " << CK.n_elem << ", TJ size: " << TJ.n_elem << std::endl;
27+
Rcpp::Rcout << "M: " << M << ", n: " << n << ", nL: " << nL << std::endl;
28+
2529
arma::vec lb(nL + 1);
2630
arma::vec ub(nL + 1);
2731

28-
for (int ll = 1; ll <= nL + 1; ++ll) {
32+
for (arma::uword ll = 1; ll <= nL + 1; ++ll) {
2933
double s;
3034
try {
31-
arma::mat P_sub = PP.cols(CJ(ll-1), CJ(ll) - 1);
32-
arma::mat B_sub = BB.cols(CK(ll-1), CK(ll) - 1);
35+
arma::uword start_col = static_cast<arma::uword>(CJ(ll-1));
36+
arma::uword end_col = static_cast<arma::uword>(CJ(ll) - 1);
37+
arma::mat P_sub = PP.cols(start_col, end_col);
38+
start_col = static_cast<arma::uword>(CK(ll-1));
39+
end_col = static_cast<arma::uword>(CK(ll) - 1);
40+
arma::mat B_sub = BB.cols(start_col, end_col);
41+
42+
Rcpp::Rcout << "Submatrix dimensions for ll=" << ll << ": P_sub(" << P_sub.n_rows << "x" << P_sub.n_cols
43+
<< "), B_sub(" << B_sub.n_rows << "x" << B_sub.n_cols << ")" << std::endl;
44+
3345
s = shat(P_sub, B_sub);
46+
Rcpp::Rcout << "shat result for ll=" << ll << ": " << s << std::endl;
3447
} catch (const std::exception& e) {
35-
Rcpp::warning("Exception caught in jhat: %s", e.what());
48+
Rcpp::warning("Exception caught in jhat for ll=%d: %s", ll, e.what());
3649
s = 1e-20;
3750
}
3851
double J = TJ(ll-1);
@@ -43,10 +56,10 @@ Rcpp::List jhat(const arma::mat& PP, const arma::mat& BB,
4356
ub(nL) = arma::datum::inf;
4457

4558
double threshold = 2 * M * std::sqrt(static_cast<double>(n));
59+
Rcpp::Rcout << "Threshold: " << threshold << std::endl;
4660

4761
arma::uvec L = arma::find(lb <= threshold && threshold <= ub);
48-
49-
int LL;
62+
arma::uword LL;
5063
int flag;
5164

5265
if (!L.empty()) {
@@ -63,7 +76,10 @@ Rcpp::List jhat(const arma::mat& PP, const arma::mat& BB,
6376
}
6477
}
6578

66-
LL = std::max(LL, 1);
79+
LL = std::max(LL, static_cast<arma::uword>(1));
80+
81+
Rcpp::Rcout << "Final LL: " << LL << ", flag: " << flag << std::endl;
82+
Rcpp::Rcout << "Exiting jhat function" << std::endl;
6783

6884
return Rcpp::List::create(
6985
Rcpp::Named("LL") = LL,

‎src/jhat.o

46.6 KB
Binary file not shown.

‎src/jlep.o

57.5 KB
Binary file not shown.

‎src/npiv.o

12.8 KB
Binary file not shown.

0 commit comments

Comments
 (0)