initial fuzzy c-means implementation (#11)

* initial genfis implementation

* update tests and docs

* update changelog

Author: lucaferranti

docs/make.jl

@@ -167,6 +167,7 @@ makedocs(;
                  ],
                  "Membership functions" => "api/memberships.md",
                  "Reading/Writing" => "api/readwrite.md",
+                 "Learning fuzzy models" => "api/genfis.md",
              ],
              "Contributor's Guide" => "contributing.md",
              "Release notes" => "changelog.md",

docs/src/api/genfis.md

@@ -0,0 +1,5 @@
+# Learning fuzzy models
+
+```@docs
+fuzzy_cmeans
+```

docs/src/changelog.md

@@ -4,7 +4,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 
 ## unreleased
 
-- added support Lukasiewicz, drastic, nilpotent and Hamacher T-norms and corresponding S-norms.
+- ![][badge-feature] Added fuzzy c-means
+- ![][badge-enhancement] added support Lukasiewicz, drastic, nilpotent and Hamacher T-norms and corresponding S-norms.
 
 ## v0.1.0 -- 2023-01-10
 
@@ -22,4 +23,13 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
   - linear and constant output for Sugeno
 - initial plotting functionalities
   - plotting variables and membership functions
-  - plotting rules of fuzzy inference system
+  - plotting rules of fuzzy inference system
+
+[badge-breaking]: https://img.shields.io/badge/BREAKING-red.svg
+[badge-deprecation]: https://img.shields.io/badge/deprecation-orange.svg
+[badge-feature]: https://img.shields.io/badge/new%20feature-green.svg
+[badge-enhancement]: https://img.shields.io/badge/enhancement-blue.svg
+[badge-bugfix]: https://img.shields.io/badge/bugfix-purple.svg
+[badge-security]: https://img.shields.io/badge/security-black.svg
+[badge-experimental]: https://img.shields.io/badge/experimental-lightgrey.svg
+[badge-maintenance]: https://img.shields.io/badge/maintenance-gray.svg

src/FuzzyLogic.jl

@@ -11,6 +11,7 @@ include("InferenceSystem.jl")
 include("parser.jl")
 include("evaluation.jl")
 include("plotting.jl")
+include("genfis.jl")
 
 export DifferenceSigmoidMF, LinearMF, GeneralizedBellMF, GaussianMF, ProductSigmoidMF,
        SigmoidMF, TrapezoidalMF, TriangularMF, SShapeMF, ZShapeMF, PiShapeMF,
@@ -19,5 +20,6 @@ export DifferenceSigmoidMF, LinearMF, GeneralizedBellMF, GaussianMF, ProductSigm
        MinImplication, ProdImplication,
        MaxAggregator, ProbSumAggregator, CentroidDefuzzifier, BisectorDefuzzifier,
        @mamfis, MamdaniFuzzySystem, @sugfis, SugenoFuzzySystem,
-       LinearSugenoOutput, ConstantSugenoOutput
+       LinearSugenoOutput, ConstantSugenoOutput,
+       fuzzy_cmeans
 end

src/genfis.jl

@@ -0,0 +1,56 @@
+"""
+Performs fuzzy clustering on th data `X` using `N` clusters.
+
+### Input
+
+- `X` -- ``d × M`` matrix of data, each column is a data point
+- `N` -- number of clusters used.
+
+### Keyword argumes
+
+- `m` -- exponent of the fuzzy membership function, default `2.0`
+- `maxiter` -- maximum number of iterations, default `100`
+- `tol` -- absolute error for stopping condition. Stop if ``|Eₖ - Eₖ₊₁|≤tol``, where ``Eₖ``
+           is the cost function value at the ``k``:th iteration.
+### Output
+
+- `C` -- ``d × N``matrix of centers, each column is the center of a cluster.
+- `U` -- ``M × N`` matrix of membership degrees, `Uᵢⱼ`` tells has the membership degree of
+         the ``j``th point to the ``i``th cluster.
+"""
+function fuzzy_cmeans(X::Matrix{T}, N::Int; m = 2.0, maxiter = 100,
+                      tol = 1e-5) where {T <: Real}
+    m > 1 || throw(ArgumentError("m must be greater than 1"))
+    e = 1 / (m - 1)
+    M = size(X, 2)
+    U = rand(float(T), M, N)
+    C = X * U .^ m ./ sum(U .^ m; dims = 1)
+
+    J = zero(float(T))
+    @inbounds for (j, cj) in enumerate(eachcol(C))
+        for (i, xi) in enumerate(eachcol(X))
+            J += U[i, j]^m * sum(abs2(k) for k in xi - cj)
+        end
+    end
+
+    @inbounds for i in 1:maxiter
+        for (j, cj) in enumerate(eachcol(C))
+            for (i, xi) in enumerate(eachcol(X))
+                U[i, j] = 1 /
+                          sum(sum(abs2.(xi - cj)) / sum(abs2.(xi - ck))
+                              for ck in eachcol(C))^e
+            end
+        end
+        C .= X * U .^ m ./ sum(U .^ m; dims = 1)
+
+        Jnew = zero(float(T))
+        for (j, cj) in enumerate(eachcol(C))
+            for (i, xi) in enumerate(eachcol(X))
+                Jnew += U[i, j]^m * sum(abs2(k) for k in xi - cj)
+            end
+        end
+        abs(J - Jnew) <= tol && break
+        J = Jnew
+    end
+    return C, U
+end

test/runtests.jl

@@ -6,6 +6,7 @@ testfiles = [
     "test_parser.jl",
     "test_evaluation.jl",
     "test_plotting.jl",
+    "test_genfis.jl",
     "test_aqua.jl",
     "test_doctests.jl",
 ]

test/test_genfis.jl

@@ -0,0 +1,9 @@
+using FuzzyLogic, Test
+
+@testset "fuzzy c-means" begin
+    X = [-3 -3 -3 -2 -2 -2 -1 0 1 2 2 2 3 3 3;
+         -2 0 2 -1 0 1 0 0 0 -1 0 1 -2 0 2]
+    C, U = fuzzy_cmeans(X, 2; m = 3)
+    @test sortslices(C; dims = 2)≈[-2.02767 2.02767; 0 0] atol=1e-3
+    @test_throws ArgumentError fuzzy_cmeans(X, 3; m = 1)
+end