Benchmarks¶

This page compares bulanik with other open-source libraries and with public data. It is there to show how the library behaves, including where it differs and where it does poorly. The tables are read from the result files in tests/benchmarks/data/ when the docs are built. The scripts that produce them are in benchmarks/; its README gives the commands. Offline tests in tests/benchmarks/ keep the frozen values honest.

In short

On matching formulas, every crisp MCDM method gives the same scores as pymcdm and pyrepo-mcda, to rounding, on 209 matrices. Most real differences are choices of variant; pymcdm's COPRAS uses another formula.
bulanik is stricter than the libraries about rounded weights, and it has no min-max option for TOPSIS.
Mamdani output equals simpful's and differs from an exact integral by a discretisation error.
ANFIS gets nearly all of its fit from the least-squares step. At default settings the premise step changes almost nothing on any dataset here, so gradient training adds little. This is the weakest result on the page.

1. MCDM against pymcdm and pyrepo-mcda¶

Method. bulanik and each library run on the same inputs: nine decision matrices taken from the studies named in the libraries' test suites, 200 random matrices from a seeded generator (4 to 12 alternatives, 3 to 8 criteria, mixed benefit and cost), and for AHP 100 random reciprocal matrices on Saaty's scale plus one from a study. A few of the study weight vectors were adjusted in the last entry so that they sum to 1 and every library accepts them. Ranks are compared with a tie tolerance of \(10^{-9}\). The /vector and /sum rows set a variant so that the two formulas match.

Data and licence. Matrices and expected values come from the pymcdm 1.4.0 and pyrepo-mcda 0.1.15 source distributions (both MIT; see benchmarks/THIRD_PARTY_NOTICES.md).

Produced with bulanik 0.3.0, pymcdm 1.4.0, pyrepo-mcda 0.1.15, numpy 2.5.4, scipy 1.18.1, Python 3.12.8; seed 20261105.

Method	Library	Matrices	Max abs. score diff.	Min. Spearman	Identical rankings	Classification
ahp	pymcdm	101	6.63e-13	1.0000	101 of 101	agree
ahp	pyrepo-mcda	101	6.62e-13	1.0000	101 of 101	agree
cocoso	pymcdm	208	7.11e-15	1.0000	208 of 208	agree
cocoso	pyrepo-mcda	208	7.11e-15	1.0000	208 of 208	agree
codas	pymcdm	209	1.78e-15	1.0000	209 of 209	agree
codas	pyrepo-mcda	209	1.78e-15	1.0000	209 of 209	agree
copras	pymcdm	209	8.66e-01	-1.0000	4 of 209	formula form
copras	pyrepo-mcda	209	4.44e-16	1.0000	209 of 209	agree
copras/sum	pymcdm	209	3.33e-16	1.0000	209 of 209	agree
critic	pymcdm	209	1.67e-16	1.0000	209 of 209	agree
critic	pyrepo-mcda	209	1.67e-16	1.0000	209 of 209	agree
critic/cost-types	pymcdm	209	2.50e-01	-1.0000	8 of 209	parameter default
critic/cost-types	pyrepo-mcda	209	2.50e-01	-1.0000	8 of 209	parameter default
edas	pymcdm	209	3.33e-16	1.0000	209 of 209	agree
edas	pyrepo-mcda	209	3.33e-16	1.0000	209 of 209	agree
entropy	pymcdm	209	6.99e-15	1.0000	209 of 209	agree
entropy	pyrepo-mcda	209	2.34e-14	1.0000	209 of 209	agree
marcos	pymcdm	209	3.33e-16	1.0000	209 of 209	agree
marcos	pyrepo-mcda	209	3.33e-16	1.0000	209 of 209	agree
topsis	pymcdm	209	2.26e-01	-0.2000	47 of 209	normalisation variant
topsis	pyrepo-mcda	209	2.26e-01	-0.2000	47 of 209	normalisation variant
topsis/vector	pymcdm	209	2.22e-16	1.0000	209 of 209	agree
topsis/vector	pyrepo-mcda	209	2.22e-16	1.0000	209 of 209	agree
vikor	pymcdm	209	3.33e-16	1.0000	209 of 209	agree
vikor	pyrepo-mcda	209	3.33e-16	1.0000	209 of 209	agree
waspas	pymcdm	209	1.11e-16	1.0000	209 of 209	agree
waspas	pyrepo-mcda	209	1.11e-16	1.0000	209 of 209	agree

Known differences

topsis, pymcdm (normalisation variant). pymcdm's default is min-max; bulanik always uses vector normalisation and has no option for min-max. The topsis/vector row sets pymcdm to vector normalisation and agrees.
copras, pymcdm (formula form). pymcdm 1.4.0 computes Q = S+ + (min(S-) * S-) / (S- * (min(S-) / S-)), which simplifies to S+ + S-: a larger cost raises the score. bulanik and pyrepo-mcda use Q = S+ + sum(S-) / (S- * sum(1/S-)). The copras/sum row rebuilds pymcdm's form from bulanik's extras (S_plus + S_minus, divided by its maximum) and agrees.
topsis, pyrepo-mcda (normalisation variant). pyrepo-mcda's default is min-max; bulanik always uses vector normalisation and has no option for min-max. The topsis/vector row sets pyrepo-mcda to vector normalisation and agrees.
critic/cost-types, pymcdm (parameter default). bulanik can flip cost columns in the min-max step (criteria_types=...), which changes the sign of their correlations; both libraries ignore criterion types here. Left at its default (all benefit) bulanik agrees, see the critic rows.
cocoso, pymcdm. bulanik refused one random matrix, which is why the row has 208 matrices: random-150: CoCoSo: alternative 'A2' is worst on every criterion, so k_b divides by zero
AHP consistency ratio, pymcdm (parameter default). Largest gap 0.0209; 1.0e-12 with the same random index. pymcdm uses another random-index table (0.5799, 0.8921, ... ); bulanik and pyrepo-mcda use Saaty's 0.58, 0.90, ... With the same table the ratios agree.
Ranks of tied alternatives (tie handling). For scores [0.437495, 0.437495, 0.953267, 0.440563]: bulanik's positions [3, 4, 1, 2], pymcdm [3.5, 3.5, 1.0, 2.0], pyrepo-mcda [3, 3, 1, 2]. bulanik returns a strict order in which ties keep input order; pymcdm gives tied alternatives their average rank; pyrepo-mcda gives them the same integer rank and does not skip the next one.

The formulas of the pymcdm COPRAS row come from reading its source; the study its test cites was not read here, so this page does not say which form the study prints or which is right.

Expected values from the libraries' test suites. bulanik against the numbers printed in those tests, at the rounding each test uses. "Normalised" means bulanik rejected the rounded weights until normalize_weights=True was set. The studies are the ones the test suites name; they were not read here.

Fixture	Study named by the test	Decimals	Max abs. diff.	Within rounding	Note
`pymcdm-cocoso`	Yazdani, Zarate, Zavadskas & Turskis (2019), Management Decision 57(9), 2501-2519	3	0.00035	yes
`pymcdm-codas`	Badi, Shetwan & Abdulshahed (2017), Proc. 1st ICMNEE, 395-407	4	0.00005	yes
`pymcdm-copras-1`	Zavadskas, Kaklauskas, Peldschus & Turskis (2007), Baltic J. Road Bridge Eng. 2(4), 195-203	4	0.07218	no	pymcdm's test expects its S+ + S- form; see the copras row.
`pymcdm-copras-2`	Zavadskas, Kaklauskas, Peldschus & Turskis (2007), Baltic J. Road Bridge Eng. 2(4), 195-203	4	0.12504	no	pymcdm's test expects its S+ + S- form; see the copras row.
`pymcdm-edas`	Yazdani, Torkayesh, Santibanez-Gonzalez & Otaghsara (2020), Sustainable Operations and Computers 1, 35-42	3	0.00038	yes	Normalised.
`pymcdm-marcos`	Ulutas, Karabasevic, Popovic, Stanujkic, Nguyen & Karakoy (2020), Mathematics 8(10), 1672	4	0.00005	yes	Normalised.
`pymcdm-topsis`	Opricovic & Tzeng (2004), EJOR 156(2), 445-455	3	0.26214	no	pymcdm's test uses its min-max default; bulanik uses vector normalisation.
`pymcdm-vikor-yang`	Yang & Wu (2020), IEEE Access 8, 21261-21271	4	0.00003	yes
`pymcdm-vikor-opricovic`	Opricovic & Tzeng (2004), EJOR 156(2), 445-455	9	0.00000	yes
`pymcdm-waspas`	Chakraborty, Zavadskas & Antucheviciene (2015), ECECSR 49(1), 5-22	4	0.00591	no	The fixture's weights sum to 1.0155 and pymcdm's test turns validation off. bulanik rejects them; with normalize_weights=True the scores shrink by up to 0.0059 and the ranking stays the same.
`pymcdm-entropy`	Elsayed, Dawood & Karthikeyan (2017), Int. J. Eng. Trends Technol. 46(2), 60-66	6	0.00000	yes
`pymcdm-critic`	Tus & Aytac Adali (2019), Opsearch 56(2), 528-538	3	0.00049	yes
`pymcdm-ahp`	Wang & Chin (2011), Information Sciences 181(23), 5240-5248	4	0.00003	yes
`pyrepo-vikor`	Papathanasiou & Ploskas (2018), Multiple Criteria Decision Aid, 31-55	3	0.00041	yes
`pyrepo-topsis`	Ersoy (2021), LogForum 17(3)	2	0.00075	yes
`pyrepo-codas-1`	Badi, Abdulshahed & Shetwan (2018), Decision Making: Appl. Manag. Eng. 1(1), 1-12	4	0.00005	yes
`pyrepo-codas-2`	Badi, Ballem & Shetwan (2018), Int. J. Quality Research 12(3), 609	4	0.00005	yes
`pyrepo-edas`	Ersoy (2021), LogForum 17(3)	3	0.00034	yes
`pyrepo-marcos`	Stevic, Pamucar, Puska & Chatterjee (2020), Computers & Industrial Engineering 140, 106231	3	0.00050	yes	Normalised.
`pyrepo-copras`	Goswami & Mitra (2020), Int. J. Data and Network Science 4(1), 27-42	5	0.00000	yes
`pyrepo-cocoso`	Yazdani, Zarate, Zavadskas & Turskis (2019), Management Decision 57(9), 2501-2519	1	0.00131	yes
`pyrepo-entropy-xu`	Xu (2004), EJOR 156(2), 530-532	4	0.00005	yes
`pyrepo-entropy-podvezko`	Zavadskas & Podvezko (2016), IJITDM 15(2), 267-283	4	0.00004	yes
`pyrepo-entropy-ersoy`	Ersoy (2021), LogForum 17(3)	3	0.00043	yes
`pyrepo-entropy-lee`	Lee & Chang (2018), Renewable and Sustainable Energy Reviews 92, 883-896	3	0.00046	yes
`pyrepo-critic`	Tus & Aytac Adali (2019), Opsearch 56(2), 528-538	3	0.00049	yes

Edge cases. The libraries and bulanik differ most at the edges. These rows show only the cases where the three do not give the same outcome. Numbers are scores; a COPRAS row that differs only as above is left out.

Case	Method	bulanik	pymcdm	pyrepo-mcda
constant column	vikor	0.511, 0.000, 1.000, 0.785	raises ValueError	0.511, 0.000, 1.000, 0.785
constant column	cocoso	1.691, 3.630, 1.534, 1.957	1.705, 2.662, 1.648, 1.865	1.705, 2.662, 1.648, 1.865
constant column	critic	0.492, 0.000, 0.508	returns non-finite	returns non-finite
zero entry in a benefit column	codas	-2.229, 2.710, -0.549, 0.068	raises ValueError	-2.229, 2.710, -0.549, 0.068
zero entry in a benefit column	waspas	raises ValueError	raises ValueError	0.112, 0.953, 0.441, 0.558
zero entry in a benefit column	marcos	raises ValueError	0.218, 0.926, 0.460, 0.547	0.218, 0.926, 0.460, 0.547
zero entry in a benefit column	entropy	0.729, 0.053, 0.217	raises ValueError	0.729, 0.053, 0.217
alternative worst on every criterion	cocoso	raises ValueError	returns non-finite	returns non-finite
benefit criteria only	copras	0.802, 1.000, 0.739, 0.647	raises ValueError	returns non-finite
all alternatives identical	topsis	0.500, 0.500, 0.500	returns non-finite	returns non-finite
all alternatives identical	vikor	0.000, 0.000, 0.000	raises ValueError	-0.600, -0.600, -0.600
all alternatives identical	edas	0.500, 0.500, 0.500	returns non-finite	returns non-finite
all alternatives identical	cocoso	raises ValueError	1.985, 1.985, 1.985	1.985, 1.985, 1.985
all alternatives identical	entropy	raises ValueError	0.333, 0.333, 0.333	0.333, 0.333, 0.333
all alternatives identical	critic	raises ValueError	returns non-finite	returns non-finite

Where bulanik is stricter: it refuses a zero in MARCOS, which both libraries accept, and in WASPAS, which pyrepo-mcda accepts and pymcdm refuses. For all-identical alternatives it refuses CoCoSo and entropy, which the libraries score. Where an alternative is worst on every CoCoSo criterion it raises an error, and the libraries return non-finite scores. On a constant column CoCoSo sets the normalised value to 0 where both libraries set it to 1, which changes the scores (not the order, in this case).

2. Human Development Index¶

Method. The HDI of a country is the geometric mean of three dimension indices: health, education and income. The published file holds the raw indicators and not the indices, so they are computed from life expectancy, expected and mean years of schooling, and GNI per capita, as \((LE - 20) / 65\), \((\min(EYS, 18) / 18 + \min(MYS, 15) / 15) / 2\) and \((\ln \min(GNI, 75000) - \ln 100) / (\ln 75000 - \ln 100)\). These goalposts were not read from UNDP's technical notes for this benchmark: they are checked by reproducing the published HDI (below).

Countries are then ranked by waspas with equal weights and lam=0, the product form, and compared with UNDP's rank. WASPAS first divides each column by its maximum, so its product form is the HDI divided by a constant, and it should match nearly exactly by construction. This checks the implementation against a published ranking. It does not show that WASPAS learns anything about countries. The sum form (lam=1) and the default blend show how far other forms drift.

Data and licence. HDR 2025 composite indices, United Nations Development Programme, used under CC BY 3.0 IGO as stated in its terms of use. A 2023 extract with computed index columns is committed at tests/benchmarks/data/hdi_2023_extract.csv; benchmarks/THIRD_PARTY_NOTICES.md lists the changes. UNDP does not endorse this work.

193 countries for 2023; left out for missing data: Korea (Democratic People's Rep. of), Monaco. Recomputed HDI equals the published value at 3 decimals for 193 of 193 countries.

Form	lam	Spearman with UNDP rank	Same rank, unrounded	Largest rank gap, unrounded	Same rank after rounding to 3 decimals
product	0.0	0.999982	159 of 193	2.0	193 of 193
blend_0.5	0.5	0.999902	118 of 193	4.0	125 of 193
sum	1.0	0.999765	86 of 193	5.0	88 of 193

Known differences. UNDP ranks on the HDI rounded to 3 decimals, so 34 countries share a rank with another (1, 2, 2, 4, ...). Unrounded scores split those ties, which accounts for the whole gap of the product form: once its values are rounded to 3 decimals, its ranks are UNDP's.

3. ANFIS on Mackey-Glass¶

Method. The series is generated here by RK4 from \(\dot x = 0.2\,x(t-17) / (1 + x(t-17)^{10}) - 0.1\,x(t)\) with step 0.1, \(x(0) = 1.2\) and \(x(t) = 0\) for \(t < 0\). Inputs are \(x(t-18), x(t-12), x(t-6), x(t)\), the target is \(x(t+6)\), with 500 training points from \(t = 118\) and the next 500 for testing. Two terms per input give 16 rules. Training is fit, the hybrid method; it is pure Python and needs no extra (the sezgi extra is only for the optional derivative-free train_anfis). Training is deterministic, so there is no seed. Every step size tried is listed; none was selected. No figure from the literature is quoted, because none was read for this benchmark and the series could differ.

Data and licence. Generated by tests/benchmarks/_bench_support.py; no external data.

Model	lr	Epochs	Train RMSE	Test RMSE	Fit seconds
Last value carried forward			0.18550	0.18476
Linear least squares			0.09516	0.09830
ANFIS gaussian, 16 rules	0.01 (default)	1	0.00217	0.00300	0.17
ANFIS gaussian, 16 rules	0.01 (default)	500	0.00217	0.00300	50.33
ANFIS bell, 16 rules	0.01 (default)	1	0.00281	0.00359	0.18
ANFIS bell, 16 rules	0.01 (default)	500	0.00281	0.00359	52.48
ANFIS gaussian, 16 rules	0.1	1	0.00217	0.00300	0.17
ANFIS gaussian, 16 rules	0.1	500	0.00217	0.00300	50.37
ANFIS bell, 16 rules	0.1	1	0.00281	0.00359	0.17
ANFIS bell, 16 rules	0.1	500	0.00279	0.00359	51.07
ANFIS gaussian, 16 rules	1	1	0.00217	0.00300	0.19
ANFIS gaussian, 16 rules	1	500	0.00217	0.00300	52.11
ANFIS bell, 16 rules	1	1	0.00281	0.00359	0.18
ANFIS bell, 16 rules	1	500	0.00267	0.00352	51.8
ANFIS gaussian, 16 rules	10	1	0.00217	0.00300	0.2
ANFIS gaussian, 16 rules	10	500	0.00217	0.00300	52.89
ANFIS bell, 16 rules	10	1	0.00281	0.00359	0.2
ANFIS bell, 16 rules	10	500	0.00270	0.00354	56.55
ANFIS gaussian, 16 rules	100	1	0.00217	0.00300	0.24
ANFIS gaussian, 16 rules	100	500	0.00217	0.00300	54.22
ANFIS bell, 16 rules	100	1	0.00281	0.00359	0.23
ANFIS bell, 16 rules	100	500	0.00268	0.00353	56.36

What this shows. At the default lr=0.01, the fit after one epoch and after 500 are the same to five decimals: the least-squares consequents do the work and the premise barely moves. The premise gradient is correct (an offline test checks it against finite differences on a 100-point subset, as the gradient of the mean squared error), but its entries there are around \(10^{-7}\) once the error is this small, so a small step does nothing. Larger values of lr help the bell terms a little and the Gaussian terms not at all, and the Gaussian rows have the best test RMSE at every lr. The test set goes slightly outside the training range, so the extrapolation warning fires on every run.

4. Mamdani and Sugeno against simpful and scikit-fuzzy¶

Method. Two rule bases written for this benchmark (tests/benchmarks/_inference_cases.py), with triangular, trapezoidal and Gaussian terms, are evaluated over a grid of input pairs. All libraries use AND as stated, OR max, NOT \(1-\mu\), clip implication and max aggregation.

Scope of scikit-fuzzy. skfuzzy.control imports networkx, which scikit-fuzzy 0.5.0 does not declare and this suite does not install, so its control engine was not run. The scikit-fuzzy column uses the library's membership functions and defuzz, with firing strengths and aggregation written in the benchmark script. No file in the package mentions Sugeno, so it has no Sugeno row.

Data and licence. No dataset. simpful and scikit-fuzzy outputs are frozen for the offline test.

Produced with bulanik 0.3.0, simpful 2.12.0, scikit-fuzzy 0.5.0.

Mamdani rule base	AND	Resolution	Grid points	Max diff. vs simpful	Max diff. vs scikit-fuzzy primitives	Seconds: bulanik, simpful, scikit-fuzzy
climate	min	101	441	4.26e-14	4.85e-01	0.036, 0.245, 0.09
climate	min	1001	441	1.71e-13	4.87e-02	0.301, 2.015, 0.226
climate	min	10001	25	1.56e-13	4.38e-03	0.144, 1.102, 0.089
climate	prod	1001	441	2.77e-13	4.84e-02	0.303, 1.999, 0.226
mixed	min	1001	441	2.42e-13	8.06e-02	0.272, 1.166, 0.2

Defuzzification (climate, resolution 1001)	Max diff. vs scikit-fuzzy	Mean diff.
bisector	3.20e-01	3.43e-02
mom	1.42e-14	2.71e-15
som	0.00e+00	0.00e+00
lom	0.00e+00	0.00e+00

Sugeno rule base	AND	Max diff. vs simpful	Seconds: bulanik, simpful
climate	min	2.84e-14	0.0062, 0.0617
climate	prod	2.84e-14	0.006, 0.0598
mixed	min	2.67e+01	0.0043, 0.0387
mixed, weights removed	min	2.84e-14	0.0042, 0.0388

Known differences

Mamdani against simpful. Same scores within rounding: same operators, same sampled centroid.
Centroid against scikit-fuzzy (formula form (quadrature)). bulanik and simpful take sum(x mu) / sum(mu) over the samples, a rectangle rule that counts the universe endpoints in full; skfuzzy integrates the piecewise-linear set exactly. The gap shrinks about tenfold for each tenfold rise in resolution. Applying the rectangle rule to the aggregate set built with scikit-fuzzy's functions gives bulanik's value to 2.8e-14.
Bisector against scikit-fuzzy. formula form (quadrature): bulanik returns the first sample at which the running sum reaches half; skfuzzy solves for the exact area split. On a low flat stretch the endpoint weight moves the split by more than one grid step.
Weighted Sugeno against simpful (formula form). With rule weights, bulanik computes sum(w a f) / sum(w a), simpful sum(w a f) / sum(a): simpful weights the numerator only. Without weights they agree (the row 'mixed, weights removed'). bulanik has no option for simpful's form; rebuilding it from bulanik's rule activations reproduces simpful. Largest gap of the rebuilt form: 1.4e-14.
When no rule fires. bulanik: raises NoRuleFiredError: no rule fired for output 'fan' and no default was given; the error's .. simpful: returns 0.0 (prints a warning). scikit-fuzzy's defuzz: raises EmptyMembershipError: The membership function area is empty..

For bulanik's centroid the error against the exact integral is largest when a term has non-zero membership at the end of the output universe. At a resolution of 1001 it reached about 0.08 on a universe of width 100 here. Raise resolution if that matters.

5. ANFIS on public regression data¶

Method. Five fixed folds from a seeded shuffle, the same for all models, on raw inputs. ANFIS uses defaults: Gaussian terms, lr=0.01, 20 epochs; ranges come from the training fold. The inputs were chosen before running: for Airfoil, frequency, suction-side displacement thickness and attack angle, plus a run with all five; for Auto MPG, weight, horsepower and model year. scikit-learn's models use all inputs and default settings. Nothing was tuned. The script runs only in benchmarks/.

The plain ANFIS rows are the result. The rows marked "inputs scaled" (inputs mapped to \([0, 1]\) by the training fold's range) and "lr 1e-12" (premise frozen, least squares only) are diagnostics. They were added after seeing the first results, which exposed a bug in fit, since fixed; the figures below are from the fixed code.

Data and licence. Both CC BY 4.0, downloaded at run time and not committed: Brooks, Pope & Marcolini (1989), Airfoil Self-Noise, and Quinlan (1993), Auto MPG, UCI Machine Learning Repository. Auto MPG rows with unknown horsepower are dropped.

airfoil (1503 rows)

Model	Mean RMSE	Std. over folds	Worst fold	Mean fit seconds
ANFIS 3 inputs, 2 terms (8 rules)	3.680	0.188	3.811	1.833
ANFIS 3 inputs, 2 terms (8 rules), inputs scaled	3.680	0.188	3.811	1.75
ANFIS 3 inputs, 2 terms (8 rules), lr 1e-12	3.680	0.188	3.811	1.483
ANFIS 3 inputs, 3 terms (27 rules)	4.277	2.502	8.737	8.592
ANFIS 3 inputs, 3 terms (27 rules), inputs scaled	4.277	2.501	8.735	8.365
ANFIS 3 inputs, 3 terms (27 rules), lr 1e-12	4.277	2.502	8.737	7.615
ANFIS 5 inputs, 2 terms (32 rules)	3.293	0.152	3.423	19.704
ANFIS 5 inputs, 2 terms (32 rules), inputs scaled	3.293	0.152	3.423	19.388
ANFIS 5 inputs, 2 terms (32 rules), lr 1e-12	3.293	0.152	3.423	18.58
scikit-learn LinearRegression (all inputs)	4.824	0.171	5.027	0.001
scikit-learn GradientBoostingRegressor (all inputs)	2.636	0.112	2.767	0.043

auto_mpg (392 rows)

Model	Mean RMSE	Std. over folds	Worst fold	Mean fit seconds	Folds with extrapolation warning
ANFIS 3 inputs, 2 terms (8 rules)	3.332	1.437	5.897	0.403	3
ANFIS 3 inputs, 2 terms (8 rules), inputs scaled	3.332	1.438	5.899	0.456	3
ANFIS 3 inputs, 2 terms (8 rules), lr 1e-12	3.332	1.438	5.899	0.435	3
ANFIS 3 inputs, 3 terms (27 rules)	5.183	2.703	9.706	2.23	3
ANFIS 3 inputs, 3 terms (27 rules), inputs scaled	5.137	2.599	9.457	2.392	3
ANFIS 3 inputs, 3 terms (27 rules), lr 1e-12	5.187	2.705	9.710	2.331	3
scikit-learn LinearRegression (all inputs)	3.372	0.142	3.586	0.0	0
scikit-learn GradientBoostingRegressor (all inputs)	2.853	0.377	3.389	0.026	0

The premise step adds little. The plain rows are within 0.005 of the "lr 1e-12" rows, so the fit comes from least squares. A premise step is kept only if it does not raise the training error, so on data like this, with a large target and a small input range, lr is shrunk rather than harming the model. A small case shows it:

import math
from bulanik.anfis import ANFIS

X = [[0.05 * i / 199] for i in range(200)]
y = [125 + 5 * math.sin(120 * x[0]) for x in X]

history = ANFIS(1, 2).fit(X, y, epochs=10)
print("RMSE per epoch, defaults:", [round(v, 2) for v in history])
print("least squares only (lr=1e-12):", round(ANFIS(1, 2).fit(X, y, epochs=1, lr=1e-12)[0], 3))

RMSE per epoch, defaults: [0.14, 0.14, 0.14, 0.14, 0.14, 0.14, 0.13, 0.13, 0.13, 0.13]
least squares only (lr=1e-12): 0.136

The least-squares step itself looks sound: in a one-off check on three Airfoil folds, it gave the same train and test RMSE as numpy.linalg.solve on the same design matrix, to 9 decimals.

Against scikit-learn. Gradient boosting beats every ANFIS row on both datasets. The best ANFIS rows do beat linear regression on Airfoil; on Auto MPG they are about level with it, and their worst fold is well above its worst.

Known differences. This is not a like-for-like contest: gradient boosting sees every input, and the ANFIS grid cannot afford to, since rules multiply with each input. Runtime is a result too: the five-input grid takes seconds per fold where scikit-learn's models take hundredths of a second or less.