GreyTheory

English | 繁體中文 | 日本語

GreyTheory is a Python toolkit for grey forecasting, grey relational analysis (GRA), and PCA. Use it to forecast sequences, reduce feature tables, or compare spectra and measurements.

Get started

Requires Python 3.11 or later. Download or clone the repository, then run from its directory:

python3 run.py

To use GreyTheory in your own program, install it from the project directory with that program's Python environment:

python -m pip install .

Time-series forecasting

Pass observations to fit, then call predict(2) for the next two values:

from greytheory import GM11

series = [223.3, 227.3, 230.5, 238.1, 242.9, 251.1]
model = GM11(background_mode="corrected").fit(series)
print(model.predict(2))
print(model.metrics)

Use at least three observations. See Convolution below for two-level window forecasting.

Read coefficients from model.coefficients and solver diagnostics from model.diagnostics. set_data stores observations for a later batch run. reset() clears results; reset(clear_data=True) also clears the data.

Convolution forecasting

Set window_size to enable convolution. Each window uses GM(1,1) to predict its next value; a second GM(1,1) then fits those predictions:

Sequence → window predictions → second-level GM(1,1) → future values

from greytheory import GM11

series = [223.3, 227.3, 230.5, 238.1, 242.9, 251.1]
model = GM11(window_size=4, stride=1, background_mode="corrected").fit(series)
print(model.predict(2))

For the full two-level method, provide enough data for at least three windows. The six observations above, with width four and stride one, produce three windows. coefficients and metrics describe the final fitted level. The last window can include extra tail observations, matching the original method.

Run the comparison to see both ordinary and convolution forecasts:

.venv/bin/python examples/convolution_forecast.py

Example output:

Ordinary GM11, next two values: [256.5532, 263.1223]
Convolution GM11, next two values: [265.9278, 274.0119]

GreyWindowFeatures instead extracts local descriptors for PCA/GRA; it does not run the second forecast level.

Multiple-factor models

GM0N and GM1N model several factors against one target. Each row is an observation and each column is a factor:

from greytheory import GM0N

factors = [[.75, .5], [1.22, 1], [.2, .7], [1, .66], [1, 1], [1, .5]]
targets = [1, 1, 1, 1, 1, 1]
model = GM0N().fit(factors, targets, feature_names=["factor_a", "factor_b"])
print(model.coefficients)
print(model.influence_degrees)

For GM(1,N), use GM1N().fit(factors, targets). Feature names are optional. influence_degrees ranks absolute coefficients, which depend on units and collinearity. Use GRA below to compare sequences with a reference.

Grey relational analysis (GRA)

Supply a reference sequence and candidate sequences. Higher scores mean closer matches under the chosen normalization:

from greytheory.analysis import grey_relation

scores = grey_relation(
    reference=[0.2, 0.5, 0.8],
    candidates=[[0.2, 0.5, 0.8], [0.8, 0.1, 0.2]],
    rho=0.5,
    normalize="none",
)
print(scores)

normalize accepts none, minmax, zscore, initial, or mean. You can also import GreyGRA from greytheory.analysis and use .score(...) or .rank(...).

Scores depend on the candidate batch, so adding or removing candidates can change the ranking. They measure grey relation, not causal importance or prediction accuracy.

PCA and feature selection

GreyFeatureSelector selects original columns. GreyPCA reduces dimensions, then uses GRA to select components against a target. This example uses synthetic data:

import numpy as np
from sklearn.model_selection import train_test_split
from greytheory.analysis import GreyFeatureSelector, GreyPCA

random_generator = np.random.default_rng(7)
factors = random_generator.normal(size=(48, 6))
targets = factors[:, 0] + 0.05 * random_generator.normal(size=48)
train_factors, test_factors, train_targets, _ = train_test_split(
    factors, targets, random_state=7
)

selector = GreyFeatureSelector(k=2).fit(train_factors, train_targets)
selected_factors = selector.transform(test_factors)

analysis = GreyPCA(
    n_components=4, method="sparse", method_params={"alpha": 0.1},
    top_k=2, standardize=True,
).fit(train_factors, train_targets)
components = analysis.transform(test_factors)
print(analysis.selected_components)
print(analysis.component_scores)

n_components sets the decomposition size; top_k sets how many components GRA keeps. Omit the target in fit(train_factors) to run dimensionality reduction alone.

Switch PCA methods

Change method and pass its settings in method_params:

For example, switch to AdvPCA:

analysis = GreyPCA(
    3, method="adv", method_params={"delta": 0.01, "max_iter": 1000},
    top_k=2, standardize=True,
).fit(train_factors, train_targets)
components = analysis.transform(test_factors)
print(analysis.reducer.converged)
print(analysis.reducer.inner_converged)

Check both outer and final inner convergence before evaluating AdvPCA results. See grey_pca_methods.py for a comparison of the methods.

You can also supply a sklearn-compatible model with its own component count:

from sklearn.decomposition import IncrementalPCA

analysis = GreyPCA(
    estimator=IncrementalPCA(n_components=4, batch_size=12), top_k=2,
).fit(train_factors, train_targets)

Custom models need fit, transform, and sklearn clone support, and must return a dense 2D matrix with one row per observation. Use either estimator or method/method_params.

Inspect selected components and rankings

component_scores contains GRA scores. selected_components contains zero-based indices; output columns retain their original component order. To list scores from highest to lowest:

valid = np.flatnonzero(np.isfinite(analysis.component_scores))
ranking = valid[np.argsort(-analysis.component_scores[valid], kind="stable")]
for index in ranking:
    print(f"PC{index + 1}", float(analysis.component_scores[index]))

Original-column scores are in selector.scores. Without a target, component_scores is None. This selection works with continuous targets. Writer IDs and other categorical labels need a separate classification or comparison model.

Choose the component count

For ordinary PCA, you can specify the cumulative variance to retain, such as 85%:

ordinary = GreyPCA(n_components=0.85, method="pca", standardize=True).fit(train_factors)
print(ordinary.n_components_fitted)
print(ordinary.explained_variance_ratio.sum())

report = analysis.reconstruction_report(train_factors, threshold=0.85)
print(report.prefix_coverage)
print(report.first_prefix_at_threshold)
print(report.selected_coverage)

Fractional n_components is supported only for ordinary PCA. For AdvPCA or SparsePCA, inspect reconstruction coverage instead. An 80–85% threshold is a starting point; use held-out data to check prediction performance, since low-variance components can still carry quality signals.

Quality scoring

For example, lower warpage and higher shear force can both be desirable. Convert them to 0–1 benefit values, then compute a weighted GRA grade:

from greytheory.analysis import rank_quality, grade_quality

quality = rank_quality(
    [[80, 10], [40, 20], [20, 25]],  # Rows: samples. Columns: warpage, shear force.
    directions=["smaller", "larger"],
    bounds=[(0, 100), (0, 30)],      # Example bounds; replace with process-specific values.
    weights=[0.5, 0.5],
)
print(quality.normalized_values)
print(quality.scores)
print(quality.ranking)  # [2, 1, 0]: row indices, best first

new_benefits = quality.normalizer.transform([[30, 22]])
new_quality = grade_quality(new_benefits, weights=quality.weights, rho=quality.rho)
print(new_quality.scores)

directions accepts larger, smaller, or target. For thickness near 100, use BenefitNormalizer(["target"], bounds=[(80, 120)], targets=[100]) with fit_transform or transform. Import it from greytheory.analysis.

bounds sets the normalization range. If omitted, it uses training minima and maxima.

Out-of-range values raise by default; out_of_range="clip" clips them to the bounds. Use normalizer.out_of_bounds(data) to locate the original violations.

Making original columns larger-is-better does not make every PCA component larger-is-better; see the component-sign notes above.

See semiconductor_adv_quality.py for process prediction. It keeps process features separate from quality targets, splits by lot, and compares AdvPCA/GRA + Ridge against a mean baseline. Do not feed warpage/shear values used to define the target back into the predictor's inputs.

coefficient = rho / (1 - benefit + rho)
score       = weighted_mean(coefficients)

Trajectory, image, and window features

Build a pipeline from numeric inputs:

Data → trajectory/image features → Grey window features → PCA/GRA → regression or distance comparison

Here is a 2D path with timestamps and stroke IDs. The same class handles 3D and 5D paths:

from greytheory.features import Trajectory, TrajectoryFeatures, GreyWindowFeatures

path = Trajectory(
    points=[[0, 0], [1, 1], [2, 1], [10, 10], [11, 10], [12, 11]],
    timestamps=[0, 0.5, 1, 3, 3.5, 4],
    stroke_ids=[0, 0, 0, 1, 1, 1],
)
sequences = TrajectoryFeatures(n_points=16, include_time=True).fit_transform([path])
windows = GreyWindowFeatures(window_size=8, stride=4, boundary_channel=-1)
features = windows.fit_transform(sequences)
print(features.shape)

paths_3d = [[[0, 0, 0], [1, 2, 1], [2, 2, 3]]]
sequences_3d = TrajectoryFeatures(n_points=16).fit_transform(paths_3d)
paths_5d = [[[0, 0, 0, 0, 0], [1, 2, 1, 3, 2], [2, 2, 3, 4, 5]]]
sequences_5d = TrajectoryFeatures(n_points=16).fit_transform(paths_5d)

To read image files directly, see “Start with image files” below.

To add regression, use sklearn's Pipeline. Here, train_paths and test_paths are lists of Trajectory objects; train_targets contains one continuous value per training path:

from sklearn.linear_model import Ridge
from sklearn.pipeline import Pipeline
from greytheory.analysis import GreyPCA

model = Pipeline([
    ("trajectory", TrajectoryFeatures(n_points=32, include_time=True)),
    ("windows", GreyWindowFeatures(8, stride=4, boundary_channel=-1)),
    ("reducer", GreyPCA(4, top_k=2, standardize=True)),
    ("model", Ridge(alpha=1)),
])
model.fit(train_paths, train_targets)
predictions = model.predict(test_paths)

See grey_feature_pipeline.py for a runnable example. These features can feed a model; handwriting authorship and overpainting detection still need dedicated models and validation data.

Start with image files

Install image support once:

.venv/bin/python -m pip install -e ".[images]"

Supply at least three reference images and a separate query image. Choose any directory; this example uses:

test_data/images/
├── reference/
│   ├── sample_a.png
│   ├── sample_b.jpg
│   └── sample_c.png
└── query.png

load_images accepts one path or a list of paths and returns an RGB batch in [0, 1]. Read PNG/JPG files and rank references by their distance from the query:

from pathlib import Path
import numpy as np
from sklearn.metrics import pairwise_distances
from sklearn.pipeline import Pipeline
from greytheory.io import load_images
from greytheory.features import ImagePatchFeatures, GreyWindowFeatures
from greytheory.analysis import GreyPCA

reference_paths = sorted(
    path for path in Path("test_data/images/reference").iterdir()
    if path.is_file() and path.suffix.lower() in {".png", ".jpg", ".jpeg"}
)
reference_images = load_images(reference_paths, size=(128, 128))
query_image = load_images("test_data/images/query.png", size=(128, 128))

pipeline = Pipeline([
    ("patches", ImagePatchFeatures(grid_shape=(4, 8))),
    ("windows", GreyWindowFeatures(8, stride=8, boundary_channel=-1)),
    ("pca", GreyPCA(2, standardize=True)),
])
reference_components = pipeline.fit_transform(reference_images)
query_components = pipeline.transform(query_image)
distances = pairwise_distances(query_components, reference_components)[0]
for index in np.argsort(distances, kind="stable"):
    print(reference_paths[index].name, float(distances[index]))

Distances are measured in PCA feature space; smaller means closer. Fit on reference images only, then transform the query. These distances do not establish authorship, authenticity, or overpainting.

size=(128, 128) is (width, height) and preserves aspect ratio with white padding. Without it, all image dimensions must match. Loading applies EXIF orientation and composites transparency over white; use mode="L" for grayscale. Align crops and positions beforehand. Keep high-bit-depth scientific intensities as numeric arrays rather than converting them through this display-image loader.

Or run the script directly. The first command generates demo PNGs; the second uses your files:

.venv/bin/python examples/image_files.py

.venv/bin/python examples/image_files.py \
  --references test_data/images/reference \
  --query test_data/images/query.png \
  --size 128 128

Results are written to artifacts/image_files/ by default:

Full example: image_files.py.

Bruker file readers

from greytheory.io.bruker import read_bcf, read_spx

spectrum = read_spx("reference.spx")
image = read_bcf("scan.bcf", lazy=True)
counts = image.get_pixel(y=0, x=0, detector=0)
energy_keV = image.energy_axes[0]
roi = image.get_spectra((0, 0, 2, 2), detector=0)

SPX exposes counts, metadata, and energy calibration when present. BCF keeps each detector's calibration. ROI bounds are (y0, x0, y1, x1), excluding the upper bounds; counts use uint64.

BCF support covers indexed, instructively packed pixels in uncompressed SFS streams. Other encodings, outer compression, and BCF v2 are unsupported. sum_detectors rejects mismatched energy axes. Check readings against the instrument software for your files.

Process and spectroscopy tools

Examples: synthetic spectra, held-out manufacturing lots, AdvPCA quality prediction. Converting XRF peak amplitudes to concentrations requires instrument, standard-sample, and matrix calibration.

Batch execution

Reuse a model with different data. ModelQueue saves a deep copy when each entry is added:

from greytheory import GM11, ModelQueue

if __name__ == "__main__":
    queue = ModelQueue()
    model = GM11(background_mode="corrected")
    for series in ([2, 3, 4, 5, 6], [3, 4, 6, 8, 11]):
        queue.append(model.set_data(series))
    results = queue.run(workers=2)
    print([item.predict(2) for item in results])

queue.run() is sequential by default. Set workers=2 for two processes; results keep input order. You can mix GM0N, GM1N, and GM11. For an existing model list, use run_models(models, workers=2). Keep the if __name__ == "__main__": guard when using processes.

Synthetic XRF and AI handwriting demos

The synthetic inputs and results are included. Download the project and open examples/results/synthetic_demos/index.html in a browser; Python is not needed to view them. GitHub does not render the HTML reports directly; use the images and CSV links below for an online preview.

To rerun both demos, run these commands in the project directory (macOS/Linux):

python3 -m venv .venv
.venv/bin/python -m pip install -e ".[demos]"
.venv/bin/python examples/synthetic_demos.py

On Windows, replace .venv/bin/python with .venv\Scripts\python.exe.

New results go to artifacts/synthetic_demos/; open its index.html. For one demo, append --demo xrf or --demo handwriting to the script command.

To refresh the bundled results, run .venv/bin/python examples/synthetic_demos.py --output examples/results/synthetic_demos.

Or call the modules directly:

from greytheory.demos import run_xrf_demo, run_handwriting_demo

xrf = run_xrf_demo("artifacts/synthetic_demos/xrf", seed=7)
handwriting = run_handwriting_demo("artifacts/synthetic_demos/handwriting")
print(xrf["report"])
print(handwriting["report"])

XRF exports emsa/*.msa, spx/*.spx, spectra.csv, and dataset.npz. Bruker documents EMSA exchange, but importing these files still needs testing with your actual M4/ESPRIT version. The SPX files are minimal fixtures for this project's reader, not guaranteed vendor exports. Synthetic mixing weights are not calibrated elemental wt%.

The handwriting sheet was created with imagegen; the image and prompt record are bundled. Labels describe prompt-assigned styles, not human writers. All samples share one generated sheet; style matches demonstrate the workflow, not cross-document or human-authorship accuracy. Grey means window descriptors here; categorical style IDs are never used as continuous GRA targets.

Development and tests

python -m pip install -e . -r requirements-dev.txt
python -m pytest -q
python -m build

Test inputs are generated or included in the repository; no dataset download is needed. scripts/verify_full_install.py checks the installed APIs. Examples use synthetic data to demonstrate the workflow, not establish process, chemical, or forensic accuracy.

To prepare a GitHub source ZIP, run python3 scripts/package_github.py. The archive in artifacts/github_upload/ includes source and synthetic results, excluding local environments, caches, and external instrument datasets.

Benchmark records are in benchmarks/results. Rerun the core benchmark or decomposition comparison. Recorded timings are from macOS arm64 and depend on the environment. The older GreyTheory, GreyGM*, add_pattern/add_patterns, and grey.run.gm* APIs remain available.

License

MIT. Dependencies retain their own licenses; see Third-party notices.

References

Algorithms, upstream tools, and further reading. SCRAMBLE is related research, not an included implementation.

1. Deng, J. L. (1982). Control problems of grey systems. Systems & Control Letters, 1(5), 288–294.

2. Liu, S. & Lin, Y. (2011). Grey Systems: Theory and Applications. Springer. Grey modelling and relational analysis.

3. Kayacan, E., Ulutas, B. & Kaynak, O. (2010). Grey system theory-based models in time series prediction. Expert Systems with Applications, 37(2), 1784–1789.

4. NumPy. Least-squares solver and rank diagnostics; accurate evaluation of exp(x) − 1.

5. Jolliffe, I. T. & Cadima, J. (2016). Principal component analysis: a review and recent developments. Philosophical Transactions A, 374, 20150202.

6. scikit-learn. PCA, SparsePCA and KernelPCA implementations; training/test preprocessing; custom estimator cloning.

7. Candès, E. J., Li, X., Ma, Y. & Wright, J. (2011). Robust Principal Component Analysis?. Journal of the ACM, 58(3), Article 11. Principal Component Pursuit.

8. Lin, Z., Chen, M. & Ma, Y. (2010). The Augmented Lagrange Multiplier Method for Exact Recovery of Corrupted Low-Rank Matrices. Inexact ALM solver for PCP.

9. Erichson, N. B. et al. (2020). Sparse Principal Component Analysis via Variable Projection, author manuscript. SIAM Journal on Applied Mathematics, 80(2), 977–1002. Robust sparse model: equation (17), Algorithm 2.

10. Vävinggren, D. et al. (2026). A Robust Optimization Approach to Sparse Principal Component Analysis. AdvPCA: Algorithms 1–2 and Appendix A.3.

11. Pfeiffer, P., Vana-Gür, L. & Filzmoser, P. (2025). Cellwise robust and sparse principal component analysis. SCRAMBLE; related work, not the algorithm used by RobustSparsePCA.

12. Lee, D. D. & Seung, H. S. (1999). Learning the parts of objects by non-negative matrix factorization. Nature, 401, 788–791. NMF background; this package delegates fitting to scikit-learn.

13. Camp, C. H. Jr. (2019). pyMCR: A Python Library for Multivariate Curve Resolution Analysis with Alternating Regression (MCR-AR). Journal of Research of NIST, 124, 124018.

14. Green, A. A., Berman, M., Switzer, P. & Craig, M. D. (1988). A transformation for ordering multispectral data in terms of image quality with implications for noise removal. IEEE Transactions on Geoscience and Remote Sensing, 26(1), 65–74. Spectral Python MNF API.

15. Anscombe, F. J. (1948). The transformation of Poisson, binomial and negative-binomial data. Biometrika, 35(3–4), 246–254. Variance stabilization used by PoissonPCA.

16. Wentzell, P. D. et al. (1997). Maximum likelihood principal component analysis. Journal of Chemometrics, 11(4), 339–366. fit_mlpca implements the independent-error, fixed-variance case.

17. pybaselines. Asymmetric least-squares baseline correction. The baseline_correct implementation uses Baseline.asls.

18. Wold, S., Sjöström, M. & Eriksson, L. (2001). PLS-regression: a basic tool of chemometrics. Chemometrics and Intelligent Laboratory Systems, 58(2), 109–130.

19. Rasmussen, C. E. & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press. SciPy Latin-hypercube sampling.

20. Jones, D. R., Schonlau, M. & Welch, W. J. (1998). Efficient Global Optimization of Expensive Black-Box Functions. Journal of Global Optimization, 13, 455–492. Expected-improvement acquisition.

21. xraylib. X-ray interaction data and documentation. Used for characteristic line energies and relative rates.

22. RosettaSciIO community. Bruker detector-calibration discussion, issue #516. File-format and per-detector energy calibration notes.

23. NIST/SEMATECH. Multiple responses: the desirability approach. Engineering preference directions and targets; GreyTheory uses separately documented arithmetic GRA grades, not the geometric-mean desirability function.

24. Pillow. Image.open and ImageOps. File loading, EXIF orientation, and aspect-preserving padding; installed with the images extra.