Is cuminc the only way to estimate the cumulative incidence function for an event of interest in the presence of competing risk events? The cuminc() function in R from the package cmprsk (https://cran.r-project.org/web/packages/cmprsk/cmprsk.pdf) allows us to get an estimation of the CIF for an event of interest in the presence of competing risk events.
Is this the only valid (known to this day) way of estimating the CIF when there are competing risk events? If not, is it the one one should go to by default blindly, or are there situations when it is not advised to go for it?

The cuminc() function in the R cmprsk package works with the subdistribution of each risk. That's a somewhat tricky concept. AdamO explains it nicely on this page:

The intepretation of this subdistributional hazard function is the instantaneous risk of death from cause 1 given you are either still alive, or you've already died of something else. In effect, it averages across these two possibilities in such a way that a high risk of dying previously from other causes lowers your hazard for that specific failure.

The cumulative incidence sums up the instantaneous subdistributional hazards. That provides a simple way to evaluate the associations of covariates with each of the competing risks, as each event is modeled separately. A downside is that the sum of predicted probabilities over all events can sometimes exceed 1. See this page and its links for more details.
An alternative is to work directly with cumulative numbers of each event type, if you aren't modeling covariates, or to do a combined competing risks survival model and work with cumulative hazards from that model. The R survival competing risks vignette illustrates that approach and compares it against the Fine-Gray subdistribution approach. The sum of probabilities over all events won't ever exceed 1, but the interpretation of covariate coefficients can be more complex as all events are modeled together.
