export_sef_report() no longer renders the bundled RMarkdown template in
place. rmarkdown::render() writes its intermediate files (*.knit.md,
figure directories) next to the input document, which meant the installed
package directory was written to during rendering. The template is now
copied into a per-call temporary directory and rendered there, so nothing is
written outside the R session's temporary directory (CRAN Policy). This also
fixes PDF/DOCX rendering failing outright when the R library is read-only,
as on the CRAN Debian check machines.bic + 2 * sum(entropy) (was bic - 2 * sum(entropy)). With the
lower-is-better BIC convention the entropy term must be added so that
overlapping/fuzzy partitions are penalised (Biernacki, Celeux & Govaert
2000); the previous sign rewarded them. compare_k() ICL values change
accordingly (BIC, PDI and the case study are unaffected; ICL is not used in
any manuscript figure).archaeo_sim(): the simulated date bounds date_min and date_max are now
returned as integer-valued years (via floor()), consistent with
archaeological dating conventions (reviewer suggestion, PCI Archaeo #1019).harris_from_contexts() gains exclude_contexts: a character vector of
context labels exempt from the verticality (depth-rank) penalty, for the
fills of cuts and multi-period contexts where depth does not track
chronology. Default NULL leaves the behaviour unchanged. The help page now
documents the verticality assumption and its failure modes (reviewer
suggestion, PCI Archaeo #1019).compare_k() and pdi() help pages now explain the BIC and
PDI model-selection metrics and reference Schwarz (1978); the manuscript adds
the same definitions and stresses combining statistical criteria with
archaeological judgement.read_pyarchinit() gains reperti_table (default "pyarchinit_reperti") and
reperti_geometry: when a material find is plotted as a point in
pyarchinit_reperti (matched by site + numero_inventario = id_rep, the
pyarchinit_reperti_view join), its own point x, y (and z from quota
or the geometry Z) are used instead of the US-polygon centroid. Finds without
a point fall back to the US centroid as before. Read via sf/GDAL so SpatiaLite
and PostGIS geometries are parsed correctly. This enables find-level spatial
resolution where piece-plotting is available.read_pyarchinit(): when taf is NULL (the default) and the
chronology_table (palimpsest_chronology) carries a per-US taf column,
those values now populate taf_score (units without a value default to 0.5).
This lets the pyArchInit QGIS algorithms (fit / intrusions / report) honour a
per-US taphonomic score stored alongside the absolute chronology, without any
change to the calling .rsx.read_pyarchinit() gains chronology_table (default "palimpsest_chronology"):
an optional per-US table of absolute dating with columns sito, area, us,
start, end (calendar years, BCE negative) — e.g. OxCal-calibrated ranges
from chronology_from_oxcal(). When present it overrides the free-text
datazione for the matching units (envelope of multiple dates per US), joined
on (sito, area, us) with a (sito, us) fallback.r:palimpsestrfit, r:palimpsestrintrusions,
r:palimpsestrreport) accept an optional PG_connection parameter (a libpq
DSN) and read from PostgreSQL/PostGIS when it is set, falling back to the
SQLite/Spatialite Database_file otherwise. read_pyarchinit() itself was
already database-agnostic (plain DBI). The pyArchInit tab passes its active
connection string, so PostgreSQL pyArchInit projects are supported.export_sef_report(): assembles a complete, narrated report from a
sef_fit — the report_sef() interpretive text, all applicable gg_*
diagnostic plots, and diagnostic tables (US summary, phase-transition matrix,
per-US phase assignments, model statistics) — rendered to PDF and/or DOCX from
a shipped RMarkdown template (inst/rmarkdown/sef_report.Rmd), in Italian or
English. A markdown narrative sidecar (<file>.md) is always written; when
pandoc/LaTeX are unavailable the function degrades gracefully to that
markdown plus PNG figures. tinytex added to Suggests.read_pyarchinit() now sources elevation per archaeological practice: US
elevation from pyarchinit_quote (mean quota_q per unit, with a
(sito, us) join fallback), material elevation from quota_usm, and pottery
inheriting its US elevation. The deprecated us_table.quota_abs form field is
no longer used.source = "both"|"materials"|"pottery" argument: reads
inventario_materiali_table and/or pottery_table (pottery class taken
from the first non-empty of ware/material/form); a find_source column
is added when source = "both".r:palimpsestrreport algorithm (palimpsestr_report_db.rsx): database →
fit → narrated PDF/DOCX report, with Source / Language / Format parameters.r:palimpsestrfit and r:palimpsestrintrusions gain a Source parameter.villa_romana example database builder now populates quota_usm,
pottery_table, and pyarchinit_quote to exercise the corrected sourcing.read_pyarchinit(): reads the inventario_materiali_table (finds) and
us_table (stratigraphic units) of a pyArchInit database (SQLite/Spatialite
or PostgreSQL/PostGIS) and assembles a fit_sef()-ready data.frame. Each
find inherits its stratigraphic unit's coordinates (centroid of the US
polygon supplied via us_geometry, e.g. the pyunitastratigrafiche layer),
elevation, chronology, and context. Free-text period strings (e.g. "II sec.
a.C.", "I/II sec. d.C.", "età romana") are resolved to numeric date_min/
date_max by a built-in archaeological-date parser (extendable via
date_labels); units without a digitised polygon fall back to synthetic
per-unit coordinates so the analysis can still run at unit resolution.
Validated against real pyArchInit Spatialite databases.chronology_from_oxcal(): OxCal counterpart of chronology_from_rcarbon().
Accepts either an oxcAAR::oxcalCalibrate() result
(oxcAARCalibratedDatesList, reduced by HPD / median+IQR / weighted mean) or
a generic data.frame of start/end calibrated ranges, and returns the
date_min/date_max/date_mid columns expected by fit_sef(). oxcAAR
is in Suggests.detect_intrusions() gains an intrusion_threshold argument and an
intrusion_type factor column that couples the intrusion magnitude (the
noise-component posterior) with the chronological direction: not_flagged,
residual (flagged and older-than-context), latent_feature (flagged and
younger-than-context), or outlier_in_context. Existing columns are
unchanged.gg_outliers() now colours each find by intrusion_type.gg_phase_composition(): per-phase class profile from the multinomial model
(cat_prob), as a heatmap (default) or stacked bars.gg_direction(): diverging lollipop of the directional intrusion diagnostic
(residual / older-than-context vs latent-feature / younger-than-context),
with a graceful placeholder when every find is in context.gg_outliers(): non-spatial ranking of the intrusion probability, using the
noise-component posterior when available, otherwise the heuristic composite.gg_unit_coherence(): within-stratigraphic-unit phase coherence, with an
optional comparison fit (e.g. the legacy Gaussian-class model).All four are exposed in a new "Model profile" tab of the Shiny app and are
compatible with as_plotly().
Bug fixes for noise-component model selection (addresses code review on the 0.17.0 pull request).
cv_sef() and
optimize_weights() scored only the Gaussian phases of a noise = TRUE
fit, so the uniform component was omitted and outlier held-out points were
mis-scored (their mass was missing and the Gaussian weights no longer summed
to one). The fit now stores the noise mixing weight (noise_weight), its
log-density (noise_logdens), and the global class frequencies
(cat_global), and both evaluators add the noise component to the held-out
likelihood.compare_k(..., noise = TRUE) no longer gives noise models a free parameter
relative to non-noise fits.New fit_sef(noise = FALSE, noise_prior = 0.05). When TRUE, a uniform
background component is added to the mixture (Fraley & Raftery, 1998). Finds
that fit no Gaussian phase accumulate posterior probability on it, so the fit
gains a noise_prob vector that is a genuine posterior probability of
being an outlier/intrusion — unlike the previous composite score, it is not
forced onto [0, 1] by rescaling, so a clean assemblage need not contain a
probability-1 intrusion. As a side effect the noise component absorbs extreme
finds, keeping them out of the phase centroid and variance estimates
(robust phase estimation). The reported phase_prob is then the distribution
over phases conditional on a find not being noise.
detect_intrusions() now returns the model-based posterior in its
intrusion_prob column whenever the fit was made with noise = TRUE,
falling back to the heuristic composite (rescaled entropy + energy +
inverse local SEI) otherwise. The directional columns are unchanged.
The option is propagated through bootstrap_sef() and exposed in the Shiny
app.
This release adds two model-based treatments of evidence that the EM engine previously simplified away: per-find dating uncertainty, and the stratigraphic constraint. Both are opt-in and off by default, so existing fits are unchanged.
fit_sef(chrono_uncertainty = FALSE). When TRUE, a find's dating
interval width is propagated into the likelihood instead of collapsing the
interval to its mid-point. The mid-date is modelled as observed with a
measurement variance (date_max - date_min)^2 / 12 (a uniform prior over
the interval), added to the temporal dimension's component variance in the
E-step and removed from the component variance estimate by method-of-moments
deconvolution in the M-step. Finds with wide dating ranges then rely less on
chronology and carry their dating uncertainty into the phase probabilities.
Unlike chrono_precision (which adds 1/tspan as a feature), this is a
model-based treatment of the uncertainty. The internal diag_log_density()
gains an extra_var argument.New fit_sef(strat_dynamic = FALSE, strat_beta = 1). The stratigraphic
context constraint was previously a static penalty frozen from the k-means
initialisation: if the EM moved assignments, the penalty still reflected the
starting clusters. With strat_dynamic = TRUE it becomes a Neighborhood-EM /
hidden-Markov-random-field term recomputed from the current phase posteriors
at every E-step (Ambroise & Govaert, 1997): each find is rewarded for
sharing the phase of its stratigraphic-unit neighbours, so units stay
coherent as the fit evolves. strat_beta sets the field strength; any
harris penalty continues to be applied statically.
Both options are propagated through bootstrap_sef() and exposed in the
Shiny app.
The class label is now modelled as a per-phase categorical distribution
(a Gaussian-times-multinomial mixed-type mixture), controlled by the new
class_model argument to fit_sef(). "multinomial" is the default;
"gaussian" reproduces the previous one-hot behaviour and is kept for
backward compatibility.
Up to v0.14.0 the class label was one-hot encoded into the Gaussian feature block. Because zero/one dummies make finds of different classes almost perfectly separable, this had two undesirable effects: per-component variances on the dummy columns collapsed to the floor (over-confident assignments), and — more importantly — the class label could split a single stratigraphic unit across several phases, even though a unit is one depositional event. Under the multinomial model the numeric evidence (space, depth, chronology) drives the phase assignment and the class label only tilts it softly, so units with finds of mixed classes stay coherent.
New class_smoothing argument (default 1): a Dirichlet pseudo-count,
shrunk towards the global class frequencies, that keeps every per-phase
class probability strictly positive. 0 removes the shrinkage.
New phase_composition() accessor: returns the estimated per-phase
class profile (row j = probability that a find in phase j belongs to
each class). The model-based, directly interpretable counterpart of a
phase-by-class cross-tabulation, stored on the fit as $cat_prob.
cv_sef(), optimize_weights(), and bootstrap_sef() propagate the
class model and include the categorical term in the held-out likelihood;
reorder_phases() reorders the per-phase class profiles; BIC/ICL use the
correct mixed-type parameter count (k(L-1) categorical parameters instead
of treating the dummies as Gaussian).
Note: on datasets recorded at stratigraphic-unit resolution (finds inheriting
US-centroid coordinates and US-tied chronology, e.g. the bundled
villa_romana), phase assignments remain near-certain under either model
because there are only as many distinct numeric points as there are units;
that certainty reflects the data resolution, not the model. The multinomial
model's benefit there is keeping units coherent and yielding interpretable
per-phase class profiles.
This release is a statistical-correctness pass on the core engine. Several
parameters that were silently inert are now active, and three diagnostics
that depended on the measurement scale or were distorted by single outliers
have been made scale-invariant and robust. Fits produced with non-default
weights, and the output of optimize_weights(), will differ from 0.13.0;
default-weight fits are affected only by the more careful BIC/ICL accounting
and the new multi-start default.
Domain weights now enter the model. weights = c(ws, wz, wt, wc)
previously affected only the SEI matrix, leaving fit_sef() — and hence
the cross-validated objective of optimize_weights() — completely
invariant to them. The weights now scale the corresponding standardised
feature dimensions (ws: planar coordinates, wz: depth, wt:
chronology-derived features, wc: the one-hot class block), so they
genuinely influence the mixture likelihood. Weights must be strictly
positive; they are stored on the fit and propagated by bootstrap_sef(),
cv_sef(), and optimize_weights().
New var_structure argument to fit_sef() ("diagonal", default, or
"spherical"). Free diagonal variances absorb any per-dimension rescaling,
so domain weights are not identifiable under a diagonal fit; the
"spherical" structure (one shared variance per component) makes them
identifiable. optimize_weights() now selects weights under
var_structure = "spherical" and corrects each held-out log-likelihood by
the Jacobian of the weighting (n_test * sum(log(w_d))) so configurations
with different weights are directly comparable.
Taphonomic down-weighting is no longer applied twice. fit_sef()
passed weights_obs = 1 - 0.5 * taf and the EM multiplied by
(1 - 0.5 * taf) again, so a find with taf = 1 contributed at weight
0.25 instead of the documented 0.5. It is now applied exactly once.
Removed the inert taf_weight_e term. It subtracted a per-row constant
from every mixture component, which cancels exactly in the softmax and so
had no effect on the posteriors (only on the reported log-likelihood). The
dead taf, taf_weight_m, and taf_weight_e arguments were removed from
the internal EM.
BIC and ICL are now built from the true unpenalized mixture
likelihood. With a stratigraphic penalty the EM objective is a penalized
criterion, not a log-likelihood; using it biased comparisons across K and
across context settings. model_stats$loglik is the unpenalized
likelihood; the penalized EM objective is available as
model_stats$loglik_penalized. The "spherical" structure also gets its
own (smaller) parameter count in the BIC.
SEI and ESE spatial/vertical components are now scale-invariant and
robust to duplicate coordinates. The former 1/d kernel normalised by
its maximum was dominated by the single closest pair, so one near-duplicate
pair (e.g. two finds sharing a grid square) collapsed every other spatial
affinity towards zero. Both components now use a bounded exponential kernel
exp(-d / h) with a data-driven bandwidth h (median pairwise
separation). ese() consequently no longer depends on the measurement unit
of the site grid. The eps and z_floor arguments of sei_matrix() are
deprecated and ignored.
detect_intrusions() directional envelope is now robust. The
leave-one-out unit envelope used a strict min/max over the other finds, so
a single chronological outlier in a unit could mask genuinely
residual/intrusive neighbours. A new envelope argument (default
c(0.05, 0.95)) uses robust quantiles; pass c(0, 1) for the previous
min/max behaviour.
n_init now defaults to 5 (was 1) and k-means initialisation uses
nstart = 5. The EM objective is multimodal, so a single start risked a
poor local optimum.
detect_intrusions() now returns two new columns:
direction: factor classifying each find as older_than_context,
in_context, or younger_than_context based on interval overlap
between the find's chronology and the leave-one-out chronological
envelope of its stratigraphic unit.chrono_gap: signed offset in years (negative = residual,
positive = latent intrusion, 0 = in-context).
Backward compatible: the existing intrusion_prob column is unchanged.New type_longevity() function: per-class temporal envelope of the
phases with mean posterior weight above a configurable threshold
(default 0.1). Returns longevity_min, longevity_max,
longevity_span, dominant_phase, n_finds, and a list-column
weight_matrix with the per-phase posterior weights.
New gg_longevity() Gantt-style plot for the type_longevity() output.
New chronology_from_rcarbon() adapter: converts an
rcarbon::CalDates object into the date_min/date_max/date_mid
columns expected by fit_sef(). Three reduction methods: HPD
(default), median_iqr, weighted_mean. BCE/CE sign convention
by default. rcarbon is in Suggests.
New recommend_setup() helper: inspects a dataset and returns
recommended fit_sef() flags (class_scale, taf_as_feature,
chrono_precision, residuality) based on diagnostic heuristics
(number of classes, singleton classes, US-centroid coordinates,
US-tied chronology, taf_score informativeness). Reports caveats
(e.g. directional intrusions will all be in_context when the
chronology is US-tied) and prints a ready-to-paste fit_sef()
call. Closes #2.
direction and chrono_gap.No breaking changes. All v0.12.0 user code continues to work as before.
gg_cv() (cross-validation diagnostics), gg_bootstrap() (confidence interval forest plot), gg_weights() (weight sensitivity heatmap).villa_romana dataset: realistic 300-find Roman villa with 4 phases, bioturbation, construction cuts, and residual pottery.cv_sef(), optimize_weights()).reorder_phases() applied to each replicate before ARI computation, solving label switching in confidence intervals.neighbourhood is set, energy is divided by actual neighbour count instead of n-1, removing peripheral bias.feature_matrix() now accepts external center/scale parameters and preserves scaling attributes for proper cross-validation.adjusted_rand_index() and confusion_matrix() for comparing estimated vs true phases.us_summary_table() aggregates diagnostics per stratigraphic unit; phase_transition_matrix() reveals vertical phase ordering; export_results() writes all results to CSV.harris_from_contexts() auto-generates penalties from depth ordering; read_harris() imports CSV edge lists; validate_phases_harris() checks phase-stratigraphy consistency.gg_convergence() (EM trace), gg_phase_profile() (depth vs phase), gg_confusion() (confusion matrix heatmap with ARI).bootstrap_sef() for bootstrap confidence intervals on PDI, entropy, energy, loglik, and ARI.cv_sef() for k-fold cross-validation of phase count selection.optimize_weights() for data-driven SEI weight estimation via grid search cross-validation.reorder_phases() relabels phases by mean depth (phase 1 = deepest/oldest) to solve label switching.sei_sparse() wrapper for sparse SEI computation on large datasets (n > 1000).fit_sef() gains n_init parameter for multiple random initialisations (default: 1). Best run by log-likelihood is retained.fit_sef() default em_iter increased from 25 to 100 for better convergence.$converged flag in sef_fit objects; warning issued on non-convergence.fit_sef() and compare_k() model selection metrics.sei_matrix() gains max_dist parameter for sparse computation (zeroes beyond threshold).sei_matrix() and ese() fully vectorized (50-100x faster on large datasets).as_plotly() with enriched archaeological tooltips (ID, context, phase, dating, class, entropy, energy, intrusion probability).gg_* functions now embed tooltip data for seamless plotly conversion.@examples, @family, and @seealso for all exported functions._pkgdown.yml for documentation website.cran-comments.md for CRAN submission.demo_easy, demo_moderate, demo_compressed.summary.sef_fit() and print.sef_fit() (NAMESPACE now uses S3method() instead of export()).as_sf_phase() and as_sf_links() CRS default (NA → NA_integer_) for compatibility with sf >= 1.0.Authors@R in DESCRIPTION.LazyData field from DESCRIPTION.fit_sef().as_sf_links() and enhanced as_sf_phase() with XY/XYZ support.