I hope you are well!
This is a code that I prepared and which crashes the R session (https://github.com/pachadotdev/fixest2/blob/cpp11_wip/dev/gdb-debug-4.R). The issue is that the C++ function cpp_get_fe_gnl_() doesn't like the matrix obsCluster, which is a matrix of integers. If we cast it as integer in R, the R session crashes with a 'memory not mapped' message. If we pass it as a double with no cast, the function does nothing and returns an error of 'double vs integer'.
# library(fixest2)
devtools::load_all()
# REPRODUCIBLE ERROR ---
# gravity_pois = fepois(Euros ~ log(dist_km) | Origin + Destination + Product + Year, trade)
# fixedEffects = fixef(gravity_pois)
# *** caught segfault ***
# address 0x55cf901cb924, cause 'memory not mapped'
# Traceback:
# 1: .Call(`_fixest2_cpp_get_fe_gnl_`, as.integer(Q), as.integer(N), sumFE, dumMat, as.integer(cluster_sizes), obsCluster)
# 2: cpp_get_fe_gnl(Q, N, S, dumMat, nbCluster, orderCluster)
# 3: fixef.fixest(gravity_pois)
# 4: fixef(gravity_pois)
# Possible actions:
# 1: abort (with core dump, if enabled)
# 2: normal R exit
# 3: exit R without saving workspace
# 4: exit R saving workspace
This seems to be related to an integer overflow that returns negative numbers. The question is why https://github.com/pachadotdev/fixest2/blob/cpp11_wip/src/05_01_misc_helpers.cpp#L590-L591 and the corresponding loop https://github.com/pachadotdev/fixest2/blob/cpp11_wip/src/05_01_misc_helpers.cpp#L593-L626 would introduce negative values. I created an issue for this https://github.com/pachadotdev/fixest2/issues/53.
If I see what the function does line by line, we have
# the problem is in fixef.fixest function L882 Methods.R
# this is the same as to run fixef.fixest line by line
# lines 924-1079 that do not apply for this case
gravity_pois = fepois(Euros ~ log(dist_km) | Origin + Destination + Product + Year, trade)
object <- gravity_pois
notes <- getFixest_notes()
sorted <- TRUE
nthreads <- getFixest_nthreads()
fixef.tol <- 1e-5
fixef.iter <- 10000
S <- object$sumFE
family <- object$family
fixef_names <- object$fixef_vars
fixef_id <- object$fixef_id
Q <- length(fixef_id)
N <- length(S)
id_dummies_vect <- list()
for (i in 1:Q) id_dummies_vect[[i]] <- as.vector(fixef_id[[i]])
is_ref_approx <- FALSE
isSlope <- FALSE
dumMat <- matrix(unlist(id_dummies_vect), N, Q) - 1
orderCluster <- matrix(unlist(lapply(id_dummies_vect, order)), N, Q) - 1
nbCluster <- sapply(fixef_id, max)
# check the data types
print(paste("Q", class(Q)))
print(paste("N", class(N)))
print(paste("S", class(S)))
print(paste("dumMat", class(dumMat)))
print(paste("nbCluster", class(nbCluster)))
print(paste("orderCluster", class(orderCluster)))
# ERROR 1 ----
# THIS IS THE LINE THAT BREAKS FIXEF()
fixef_values <- cpp_get_fe_gnl(Q, N, S, dumMat, nbCluster, orderCluster)
# ERROR MESSAGE:
#
# OK L570OK L5954
# 0
# OK L612OK L6321
#
# *** caught segfault ***
# address 0x5623b3cd0df0, cause 'memory not mapped'
#
# Traceback:
# 1: .Call(`_fixest2_cpp_get_fe_gnl_`, as.integer(Q), as.integer(N), sumFE, as.integer(dumMat), as.integer(cluster_sizes), as.integer(obsCluster))
# 2: cpp_get_fe_gnl(Q, N, S, dumMat, nbCluster, orderCluster)
#
# Possible actions:
# 1: abort (with core dump, if enabled)
# 2: normal R exit
# 3: exit R without saving workspace
# 4: exit R saving workspace
# NOW WE TRY TO STORE THE MATRICES AS INTEGER-TYPE
storage.mode(dumMat) <- "integer"
storage.mode(orderCluster) <- "integer"
# THE FUNCTION BREAKS AGAIN
# fixef_values <- cpp_get_fe_gnl(Q, N, S, dumMat, nbCluster, orderCluster)
# ERROR MESSAGE:
#
# OK L570OK L5954
# 0
# OK L612OK L6321
#
# *** caught segfault ***
# address 0x56144c29706c, cause 'memory not mapped'
#
# Traceback:
# 1: .Call(`_fixest2_cpp_get_fe_gnl_`, as.integer(Q), as.integer(N), sumFE, as.integer(dumMat), as.integer(cluster_sizes), as.integer(obsCluster))
# 2: cpp_get_fe_gnl(Q, N, S, dumMat, nbCluster, orderCluster)
#
# Possible actions:
# 1: abort (with core dump, if enabled)
# 2: normal R exit
# 3: exit R without saving workspace
# 4: exit R saving workspace
To run it, and make the R session crash, I would need this "recipe" to follow the exact same steps I used:
- Clone the project
git clone --depth 1 -b cpp11_wip --single-branch git@github.com:pachadotdev/fixest2.git
Open the fixest2 project in RStudio and press CTRL+SHIFT+B to install
Run the previous example
I also opened an issue for this problem: https://github.com/pachadotdev/fixest2/issues/52.
The C++ code is this (https://github.com/pachadotdev/fixest2/blob/cpp11_wip/src/05_01_misc_helpers.cpp#L540)
[[cpp11::register]] list cpp_get_fe_gnl_(int Q,
int N,
writable::doubles sumFE,
writable::integers_matrix<> dumMat,
writable::integers cluster_sizes,
writable::integers_matrix<> obsCluster)
{
// This function returns a list of the cluster coefficients for each cluster
// dumMat: the matrix of cluster ID for each observation, with cpp index style
// Q, N: nber of clusters / obs
// cluster_sizes: vector of the number of cases per cluster
// obsCluster: the integer vector that is equal to order(dum[[g]])
// RETURN:
// a list for each cluster of the cluster coefficient value
// + the last element is the number of clusters that have been set as references (nb_ref)
// => much optimized version May 2019
int iter = 0, iterMax = 10000;
int iter_loop = 0, iterMax_loop = 10000;
// Creation of the indices to put all the cluster values into a single vector
int nb_coef = 0;
std::vector<int> nb_ref(Q); // nb_ref takes the nb of elements set as ref
for (int q = 0; q < Q; q++)
{
// the total number of clusters (eg if c1: man/woman and c2: 10 countries: total of 12 cases)
nb_coef += cluster_sizes[q];
}
cout << "OK L570";
std::vector<double> cluster_values(nb_coef);
// index of the cluster
std::vector<int*> pindex_cluster(Q);
std::vector<int> index_cluster(nb_coef);
for (int i = 0; i < nb_coef; i++)
{
index_cluster[i] += i;
}
pindex_cluster[0] = index_cluster.data();
for (int q = 1; q < Q; q++)
{
pindex_cluster[q] = pindex_cluster[q - 1] + cluster_sizes[q - 1];
}
// Now we create the vector of observations for each cluster
// we need a starting and an end vector as well
std::vector<int> start_cluster(nb_coef);
std::vector<int> end_cluster(nb_coef);
cout << "OK L595";
int index;
int k;
cout << Q;
cout << "\n";
for (int q = 0; q < Q; q++)
{
cout << q;
cout << "\n";
// table cluster: nber of elements for each cluster class
writable::integers tableCluster;
tableCluster.push_back(cluster_sizes[q]);
for (int i = 0; i < N; i++)
{
// cout << cluster_sizes[q];
k = dumMat(i, q);
tableCluster[k] += 1; // the number of obs per case
}
cout << "OK L612";
// now creation of the start/end vectors
for (int k = 0; k < cluster_sizes[q]; k++)
{
int *pindex = pindex_cluster[q];
index = pindex[k];
// index = start(q) + k;
if (k == 0)
{
start_cluster[index] = 0;
end_cluster[index] = tableCluster[k];
}
else
{
start_cluster[index] = end_cluster[index - 1];
end_cluster[index] = end_cluster[index - 1] + tableCluster[k];
}
}
cout << "OK L632";
}
cout << "OK L635";
// matrix of the clusters that have been computed
writable::integers_matrix<> mat_done(N, Q);
writable::integers rowsums(N);
cout << "OK L639";
// vector of elements to loop over
writable::integers id2do(N);
writable::integers id2do_next(N);
int nb2do = N, nb2do_next = N;
for (int i = 0; i < nb2do; i++)
{
id2do[i] = i;
id2do_next[i] = i;
}
cout << "OK L658";
// Other indices and variables
int qui_max, obs;
int rs, rs_max; // rs: row sum
int id_cluster;
double other_value;
bool first;
cout << "OK L667";
//
// THE MAIN LOOP
//
while (iter < iterMax)
{
iter++;
// Finding the row where to put the 0s
if (iter == 1)
{
// 1st iter, we select the first element
qui_max = 0;
}
else
{
// we find the row that has the maximum of items done
qui_max = 0;
rs_max = 0;
for (int i = 0; i < nb2do; i++)
{
obs = id2do[i];
rs = rowsums[obs];
if (rs == Q - 2)
{
// if rs == Q-2 => its the maximum possible, no need to go further
qui_max = obs;
break;
}
else if (rs < Q && rs > rs_max)
{
// this is to handle complicated cases with more than 3+ clusters
qui_max = obs;
rs_max = rs;
}
else if (qui_max == 0 && rs == 0)
{
// used to initialize qui_max
qui_max = obs;
}
}
}
cout << "OK L716";
// Putting the 0s, ie setting the references
// the first element is spared
first = true;
for (int q = 0; q < Q; q++)
{
if (mat_done(qui_max, q) == 0)
{
if (first)
{
// we spare the first element
first = false;
}
else
{
// we set the cluster to 0
// 1) we find the cluster
id_cluster = dumMat(qui_max, q);
// 2) we get the index of the cluster vector
int *pindex = pindex_cluster[q];
index = pindex[id_cluster];
// 3) we set the cluster value to 0
cluster_values[index] = 0;
// 4) we update the mat_done matrix for the elements of this cluster
for (int i = start_cluster[index]; i < end_cluster[index]; i++)
{
obs = obsCluster(i, q);
mat_done(obs, q) = 1;
rowsums[obs]++;
}
// 5) => we save the information on which cluster was set as a reference
nb_ref[q]++;
}
}
}
// LOOP OF ALL OTHER UPDATES (CRITICAL)
iter_loop = 0;
while (iter_loop < iterMax_loop)
{
iter_loop++;
R_CheckUserInterrupt();
// Selection of indexes (new way) to be updated
// initialisation of the observations to cover (before first loop the two are identical)
if (iter_loop != 1)
{
nb2do = nb2do_next;
for (int i = 0; i < nb2do; i++)
{
id2do[i] += id2do_next[i];
}
}
nb2do_next = 0;
for (int i = 0; i < nb2do; i++)
{
// we compute the rowsum of the obs that still needs to be done
obs = id2do[i];
rs = rowsums[obs];
if (rs < Q - 1)
{
// you need to check it next time!
id2do_next[nb2do_next] = obs;
nb2do_next++;
}
else if (rs == Q - 1)
{
// means: needs to be updated
int q;
for (q = 0; mat_done(obs, q) != 0; q++)
{
// selection of the q that is equal to 0
}
int *pindex = pindex_cluster[q];
int index_select = pindex[dumMat(obs, q)];
// Finding the value of the cluster coefficient
other_value = 0.0;
// Computing the sum of the other cluster values
// and finding the cluster to be updated (q_select)
for (int l = 0; l < Q; l++)
{
// we can loop over all q because cluster_values is initialized to 0
index = pindex_cluster[l][dumMat(obs, l)];
other_value += cluster_values[index];
}
// the index to update
cluster_values[index_select] += sumFE[obs] - other_value;
// Update of the mat_done
for (int j = start_cluster[index_select]; j < end_cluster[index_select]; j++)
{
obs = obsCluster(j, q);
mat_done(obs, q) = 1;
rowsums[obs]++;
}
}
}
if (nb2do_next == nb2do)
break;
}
// out of this loop:
// - nb2do == nb2do_next
// - id2do == id2do_next
// Check that everything is all right
if (iter_loop == iterMax_loop)
{
Rprintf("Problem getting FE, maximum iterations reached (2nd order loop).");
}
// if no more obs to be covered
if (nb2do_next == 0)
break;
}
cout << "OK L843";
if (iter == iterMax)
{
Rprintf("Problem getting FE, maximum iterations reached (1st order loop).");
}
// final formatting and save
writable::list res(Q + 1);
int *pindex;
for (int q = 0; q < Q; q++)
{
writable::doubles quoi(cluster_sizes[q]);
pindex = pindex_cluster[q];
for (k = 0; k < cluster_sizes[q]; k++)
{
index = pindex[k];
// index = start(q) + k;
quoi[k] = cluster_values[index];
}
res[q] = quoi;
}
res[Q] = as_sexp(nb_ref);
return res;
}
Which is called from R here https://github.com/pachadotdev/fixest2/blob/cpp11_wip/R/Methods.R#L882.
My apologies about the lack of a MWE, this seems to be about R and C+ communication so the example needs context.
Best wishes,
