Merge pull request #19 from rafaelbailo/main

Anisotropic noise

docs/make.jl

@@ -36,6 +36,7 @@ makedocs(;
       "Example objectives" => "example_objectives.md"
     ],
     "Mid-level usage" => [
+      "Noise types" => "noise_types.md"
       "Extended output" => "extended_output.md"
       "Output visualisation" => "output_visualisation.md"
       "Performance and benchmarking" => "performance_benchmarking.md"

docs/src/noise_types.md

@@ -0,0 +1,15 @@
+# Noise types
+
+By default, ConsensusBasedX.jl uses so-called *isotropic noise* (option `noise = :IsotropicNoise`), given by
+```math
+\mathrm{d}X_i(t) = \cdots + \sqrt{2} \sigma \left| X_i(t) - V(t) \right| \mathrm{d}W_i(t),
+```
+where ``W_i`` are independent Brownian processes. The intensity of the noise depends on the distance of each particle to the consensus point, ``\left| X_i(t) - V(t) \right|``.
+
+## Anisotropic noise
+
+ConsensusBasedX.jl also offers *anisotropic noise*, given by 
+```math
+\mathrm{d}X_i(t) = \cdots + \sqrt{2} \sigma \,\mathrm{diag} \left( X_i(t) - V(t) \right) \mathrm{d}W_i(t).
+```
+The intensity of the noise now varies along each dimension. This can be selected with the option `noise = :AnisotropicNoise`.

docs/src/summary_options.md

@@ -33,6 +33,7 @@ See [Stopping criteria](@ref).
 
 ## Advanced options
 
+- `noise = :IsotropicNoise` controls the type of noise, see [Noise types](@ref).
 - `benchmark::Bool = false` controls the benchmark behaviour. `benchmark = true` runs the `ParticleDynamic` on benchmark mode, see [Performance and benchmarking](@ref).
 - `extended_output::Bool = false` controls the output, and by default returns only the computed minimiser. `extended_output = true` returns additional information, see [Extended output](@ref).
 - `parallelisation = :NoParallelisation` controls the parallelisation of the `minimise` routine, switched off by default. `parallelisation=:EnsembleParallelisation` enables parallelisation, see  [Parallelisation](@ref).

src/CBO/CBO_method.jl

@@ -36,12 +36,12 @@ function compute_CBO_consensus!(
 end
 
 function compute_CBO_update!(
-  method::ConsensusBasedOptimisation,
+  method::ConsensusBasedOptimisation{TF, TCorrection, <:TIsotropicNoise},
   method_cache::ConsensusBasedOptimisationCache,
   particle_dynamic::ParticleDynamic,
   particle_dynamic_cache::ParticleDynamicCache,
   m::Int,
-)
+) where {TF, TCorrection}
   @expand particle_dynamic_cache D N X dX Δt root2Δt
   @expand method correction λ σ
   @expand method_cache consensus consensus_energy distance energy
@@ -64,3 +64,26 @@ function compute_CBO_update!(
   end
   return nothing
 end
+
+function compute_CBO_update!(
+  method::ConsensusBasedOptimisation{TF, TCorrection, <:TAnisotropicNoise},
+  method_cache::ConsensusBasedOptimisationCache,
+  particle_dynamic::ParticleDynamic,
+  particle_dynamic_cache::ParticleDynamicCache,
+  m::Int,
+) where {TF, TCorrection}
+  @expand particle_dynamic_cache D N X dX Δt root2Δt
+  @expand method correction λ σ
+  @expand method_cache consensus consensus_energy energy
+
+  for n ∈ 1:N
+    for d ∈ 1:D
+      dX[m][n][d] =
+        (consensus[m][d] - X[m][n][d]) * (
+          Δt * λ * correction(energy[m][n] - consensus_energy[m]) +
+          root2Δt * σ * randn()
+        )
+    end
+  end
+  return nothing
+end

src/CBO/ConsensusBasedOptimisation.jl

@@ -14,25 +14,35 @@ Fields:
 mutable struct ConsensusBasedOptimisation{
   TF,
   TCorrection <: ConsensusBasedXCorrection,
+  TNoise <: Noises,
 } <: ConsensusBasedXMethod
   f::TF
   correction::TCorrection
+  noise::TNoise
   α::Float64
   λ::Float64
   σ::Float64
 end
 
 @config function construct_CBO(
   f,
-  correction::ConsensusBasedXCorrection;
+  correction::ConsensusBasedXCorrection,
+  noise::Noises;
   α::Real = 10,
   λ::Real = 1,
   σ::Real = 1,
 )
   @assert α >= 0
   @assert λ >= 0
   @assert σ >= 0
-  return ConsensusBasedOptimisation(f, correction, float(α), float(λ), float(σ))
+  return ConsensusBasedOptimisation(
+    f,
+    correction,
+    noise,
+    float(α),
+    float(λ),
+    float(σ),
+  )
 end
 
 """

src/ConsensusBasedX.jl

@@ -5,6 +5,7 @@ using DefaultKeywordArguments
 
 include("./types/abstract.jl")
 include("./types/modes.jl")
+include("./types/noise.jl")
 include("./types/parallelisation.jl")
 
 include("./utils/arrays.jl")

src/ConsensusBasedXLowLevel.jl

@@ -2,13 +2,16 @@ module ConsensusBasedXLowLevel
 
 using Reexport
 
-@reexport import ..ConsensusBasedXMethod,
+@reexport import ..AnisotropicNoise,
+  ..ConsensusBasedXMethod,
   ..ConsensusBasedXMethodCache,
   ..ConsensusBasedOptimisation,
   ..ConsensusBasedOptimisationCache,
   ..HeavisideCorrection,
+  ..IsotropicNoise,
   ..NoCorrection,
   ..Modes,
+  ..Noises,
   ..Parallelisations,
   ..ParticleDynamic,
   ..ParticleDynamicCache,

src/interface/minimise.jl

@@ -66,7 +66,7 @@ end
   correction = HeavisideCorrection()
   # correction = RegularisedHeavisideCorrection(1e-3)
 
-  method = construct_CBO(config, f, correction)
+  method = construct_CBO(config, f, correction, config.noise)
   particle_dynamic = construct_particle_dynamic(config, method)
   return particle_dynamic
 end

src/interface/parse_config.jl

@@ -2,6 +2,7 @@ const DEFAULT_PARSED_CONFIG = (;
   N = 20,
   M = 1,
   mode = ParticleMode,
+  noise = IsotropicNoise,
   parallelisation = NoParallelisation,
   verbosity = 0,
 )
@@ -10,6 +11,7 @@ function parse_config(config::NamedTuple)
   check_config_has_D(config)
   config = parse_config_default(config)
   config = parse_config_mode(config)
+  config = parse_config_noise(config)
   config = parse_config_parallelisation(config)
   return config
 end
@@ -43,6 +45,25 @@ function parse_config_mode(config::NamedTuple)
   return config
 end
 
+function parse_config_noise(config::NamedTuple)
+  if (haskey(config, :noise))
+    noise = config.noise
+    noise = (noise isa String) ? Symbol(noise) : noise
+    noise = (noise isa Symbol) ? Val(noise) : noise
+    if noise isa Val
+      if !(noise isa Noises)
+        explanation = "The selected `noise` is not recognised as an instance of `ConsensusBasedX.Noises`."
+        throw(ArgumentError(explanation))
+      end
+    else
+      explanation = "The keyword `noise` should be a `Symbol`, a `String`, or an instance of `ConsensusBasedX.Noises`."
+      throw(ArgumentError(explanation))
+    end
+    return merge(config, (; noise))
+  end
+  return config
+end
+
 function parse_config_parallelisation(config::NamedTuple)
   if (haskey(config, :parallelisation))
     parallelisation = config.parallelisation

src/types/noise.jl

@@ -0,0 +1,12 @@
+const ProtoNoises = (:IsotropicNoise, :AnisotropicNoise)
+
+for Noise ∈ ProtoNoises
+  @eval begin
+    const $Noise = Val($(Meta.quot(Noise)))
+    const $(Symbol(:T, Noise)) = Val{$(Meta.quot(Noise))}
+  end
+end
+
+@eval begin
+  const Noises = Union{$(map(s -> Symbol(:T, s), ProtoNoises)...)}
+end