HSDS: Harmony Search in Dependent Space

Harmony Search Algorithm for Dependent Design Spaces.

HSDS is a Python library for solving single- and multi-objective optimization problems using the Harmony Search metaheuristic. Its key design principle is search-space first: instead of just minimizing a function, you describe the domain of each variable precisely — including dependencies between variables, discrete grids, catalog lookups, and domain-specific feasibility rules — and let the algorithm handle the rest.

from hsds import DesignSpace, Continuous, Discrete, Minimization

space = DesignSpace()
space.add("h",  Continuous(0.30, 1.20))
space.add("bf", Continuous(lo=lambda ctx: ctx["h"] * 0.5,
                           hi=lambda ctx: ctx["h"] * 2.0))
space.add("n",  Discrete(4, 2, 20))

def objective(harmony):
    h, bf, n = harmony["h"], harmony["bf"], harmony["n"]
    cost    = 1.1 * h * bf + 0.04 * n
    penalty = max(0.0, h - 2 * bf)
    return cost, penalty

result = Minimization(space, objective).optimize(
    memory_size=20, hmcr=0.85, par=0.35, max_iter=5000
)
print(result)

---

Installation

pip install pyhsds

Requires Python 3.8+. No mandatory dependencies beyond the standard library.

For development:

pip install -r requirements-dev.txt
pip install -e .

---

Core concepts

Design variables

Every variable implements three methods that the algorithm calls internally:

Method	Purpose
sample(ctx)	Draw a random feasible value
filter(candidates, ctx)	Keep only feasible values from harmony memory
neighbor(value, ctx)	Return an adjacent feasible value (pitch adjustment)

The ctx argument is a dict of all variable values assigned earlier in the same harmony. This enables dependent bounds — the domain of a variable can depend on previously assigned variables.

Built-in variable types

Type	Domain
Continuous(lo, hi)	ℝ ∩ [lo, hi]
Discrete(lo, step, hi)	{lo, lo+step, …, hi}
Integer(lo, hi)	{lo, lo+1, …, hi}
Categorical(choices)	finite label set

All bounds accept callables for dependent domains:

space.add("d",  Continuous(0.40, 1.20))
space.add("tw", Continuous(lo=lambda ctx: ctx["d"] / 50,
                           hi=lambda ctx: ctx["d"] / 10))

Domain-specific variable spaces

HSDS ships a catalog of ready-made variable types for common engineering and mathematical domains:

from hsds import ACIRebar, SteelSection, ConcreteGrade, PrimeVariable

# ACI 318 ductile bar arrangement — bounds depend on d and fc
space.add("rebar", ACIRebar(d_expr=lambda ctx: ctx["d"],
                            cc_expr=60.0,
                            fc=lambda ctx: ctx["grade"].fck_MPa,
                            fy=420.0))

# Standard steel I-section from built-in catalog (IPE, HEA, HEB, W)
space.add("section", SteelSection(series=["IPE", "HEA"]))

# EN 206 concrete grade (C12/15 to C90/105)
space.add("concrete", ConcreteGrade(min_grade="C25/30", max_grade="C50/60"))

# Prime numbers
space.add("p", PrimeVariable(lo=2, hi=500))

Full catalog:

Category	Types
Mathematical	NaturalNumber, WholeNumber, NegativeInt, NegativeReal, PositiveReal, PrimeVariable, PowerOfTwo, Fibonacci
Structural	ACIRebar, ACIDoubleRebar, SteelSection, ConcreteGrade
Geotechnical	SoilSPT
Seismic	SeismicZoneTBDY

All types are also accessible via the plugin registry:

from hsds import create_variable, list_variable_types
print(list_variable_types())
var = create_variable("aci_rebar", d_expr=0.55, cc_expr=40.0)

Custom variables

Subclass Variable for full control:

from hsds import Variable, register_variable

@register_variable("my_type")
class MyVariable(Variable):
    def sample(self, ctx):              ...
    def filter(self, candidates, ctx):  ...
    def neighbor(self, value, ctx):     ...

Factory function for quick prototyping:

from hsds import make_variable
import random

EvenVar = make_variable(
    sample   = lambda ctx: random.choice(range(2, 101, 2)),
    filter   = lambda cands, ctx: [c for c in cands if c % 2 == 0],
    neighbor = lambda val, ctx: val + random.choice([-2, 2]),
    name     = "even",
)
space.add("n", EvenVar())

---

Optimisers

Minimization

result = Minimization(space, objective).optimize(
    memory_size      = 20,       # Harmony Memory Size (HMS)
    hmcr             = 0.85,     # Harmony Memory Considering Rate
    par              = 0.35,     # Pitch Adjusting Rate
    max_iter         = 5000,
    bw_max           = 0.05,     # Initial bandwidth (5% of domain width)
    bw_min           = 0.001,    # Final bandwidth (exponential decay)
    resume           = "auto",   # "auto" | "new" | "resume"
    checkpoint_path  = "run.json",
    checkpoint_every = 500,
    use_cache        = False,    # Cache identical harmony evaluations
    cache_maxsize    = 4096,
    log_init         = False,    # Write initial memory to CSV
    log_history      = False,    # Write best-per-iteration to CSV
    log_evaluations  = False,    # Write every evaluated harmony to CSV
    history_every    = 1,
    verbose          = True,
    callback         = my_callback,
)
print(result.best_harmony)
print(result.best_fitness)

Maximization

Same interface — negates internally, reports original sign.

result = Maximization(space, objective).optimize(...)

MultiObjective

def objective(harmony):
    f1 = harmony["x"] ** 2
    f2 = (harmony["x"] - 2) ** 2
    return (f1, f2), 0.0    # tuple of objectives, penalty

result = MultiObjective(space, objective).optimize(
    max_iter     = 10_000,
    archive_size = 100,
)

for entry in result.front:
    print(entry.objectives, entry.harmony)

Callback and early stopping

def my_callback(iteration, partial_result):
    print(iteration, partial_result.best_fitness)
    if partial_result.best_fitness < 1e-4:
        raise StopIteration    # stops the loop cleanly

---

Advanced features

Dynamic bandwidth narrowing

The pitch adjustment step size decays exponentially from bw_max to bw_min over the run — wide exploration early, fine convergence late.

result = Minimization(space, objective).optimize(
    bw_max=0.10,   # 10% of domain width at iteration 0
    bw_min=0.001,  # 0.1% at final iteration
    max_iter=5000,
)

Set bw_max == bw_min for constant bandwidth (original HS behavior). Discrete and categorical variables are unaffected by bandwidth.

Resume control

# "auto"   — continue if checkpoint exists, start fresh otherwise (safe default)
# "new"    — always start fresh, overwrite any existing checkpoint
# "resume" — always continue; raises FileNotFoundError if checkpoint missing

result = optimizer.optimize(
    max_iter        = 50_000,
    checkpoint_path = "run.json",
    resume          = "auto",
)

The initial harmony memory is saved immediately at startup — even a run interrupted in the first seconds can be resumed cleanly.

Evaluation cache

Identical harmonies are never re-evaluated when use_cache=True. Particularly valuable for expensive objectives (FEM, CFD, etc.).

result = optimizer.optimize(use_cache=True, cache_maxsize=4096)
print(optimizer._cache.stats())
# EvaluationCache: 412 hits / 1005 total (41.0% hit rate)  size=593/4096

CSV logging

result = optimizer.optimize(
    checkpoint_path  = "run.json",
    log_init         = True,    # → run_init.csv     (initial memory)
    log_history      = True,    # → run_history.csv  (best per iteration)
    log_evaluations  = True,    # → run_evals.csv    (every evaluation)
    history_every    = 10,      # write history every 10 iterations
)

All CSV files are readable directly in Excel or with pandas.read_csv().

---

Decoding engineering variables

Variables like ACIRebar and SteelSection store integer codes in the harmony. Use decode() to get full properties:

rebar_var = ACIRebar(d_expr=0.55, cc_expr=40.0)
code = result.best_harmony["rebar"]
diameter_mm, bar_count = rebar_var.decode(code)
print(rebar_var.describe(code))   # "8 bars of Ø19.00 mm"

section_var = SteelSection(series=["IPE"])
sec = section_var.decode(result.best_harmony["section"])
print(sec.name, sec.Iy_cm4, "cm4")

grade_var = ConcreteGrade()
grade = grade_var.decode(result.best_harmony["concrete"])
print(grade.name, grade.fck_MPa, "MPa", grade.Ecm_GPa, "GPa")

Steel section catalog

The built-in catalog covers IPE 80–600, HEA 100–500, HEB 100–500, and W-sections. Override with your own file:

var = SteelSection(catalog="my_sections.json")  # custom catalog
var = SteelSection(series=["HEA", "HEB"])          # filter series

---

Algorithm background

HSDS implements Harmony Search with several enhancements:

Dynamic bandwidth narrowing — pitch adjustment step size decays exponentially. Early iterations explore broadly; late iterations converge precisely.

Intelligent pitch adjustment — neighbor() is called with the current dependency context so the perturbed value stays feasible. The common incorrect approach of calling sample() on PAR is avoided.

Dependent search spaces — variables are sampled in definition order; each receives a context dict of previously assigned values. Dependent bounds, catalog filters, and feasibility checks can reference earlier variables without any special handling in the optimizer loop.

Deb constraint handling — feasible solutions always rank above infeasible ones; among infeasible solutions ranking is by total penalty.

References

• Geem, Z. W., Kim, J. H., & Loganathan, G. V. (2001). A new heuristic optimization algorithm: Harmony search. Simulation, 76(2), 60–68.
• Lee, K. S., & Geem, Z. W. (2005). A new meta-heuristic algorithm for continuous engineering optimization. Computer Methods in Applied Mechanics and Engineering, 194(36–38), 3902–3933.
• Deb, K. (2000). An efficient constraint handling method for genetic algorithms. Computer Methods in Applied Mechanics and Engineering, 186(2–4), 311–338.
• Ricart, J., Hüttemann, G., Lima, J., & Barán, B. (2011). Multiobjective harmony search algorithm proposals. Electronic Notes in Theoretical Computer Science, 281, 51–67.

---

Testing

pip install -r requirements-dev.txt
pytest tests/ -v

472 tests across 15 test files covering:

• All variable types — sample, filter, neighbor, edge cases, lo > hi validation
• DesignSpace — dependency chains, empty space, 50-variable stress test
• Optimisers — Minimization, Maximization, MultiObjective
• New features — bandwidth decay, resume modes, evaluation cache, CSV logging
• Pareto archive — dominance, crowding distance, serialization
• Engineering physics — EC2 formulas, ACI 318 feasibility, steel section properties
• Determinism — same seed produces identical results
• Numerical correctness — Sphere, Rosenbrock, constrained minimization

---

Project structure

hsds/
├── hsds/
│   ├── variables.py       # Continuous, Discrete, Integer, Categorical
│   ├── space.py           # DesignSpace
│   ├── optimizer.py       # Minimization, Maximization, MultiObjective
│   ├── pareto.py          # Pareto archive, crowding distance
│   ├── registry.py        # register_variable, make_variable
│   ├── logging.py         # EvaluationCache, RunLogger
│   └── spaces/
│       ├── math.py        # Mathematical search spaces
│       └── engineering.py # Engineering domain spaces
├── examples/
│   ├── 01_quickstart.py
│   ├── 02_welded_beam.py
│   ├── 03_rc_beam_design.py
│   ├── 04_custom_variable.py
│   ├── 05_multi_objective.py
│   ├── 06_steel_beam_design.py
│   └── 07_rc_section_full.py
├── tests/                 # 472 tests across 15 files
├── requirements-dev.txt
├── ruff.toml
├── pyproject.toml
└── LICENSE

---

Citation

If you use hsds in your research, please cite it as follows:

APA:

Özcan, A. (2026). hsds: Harmony Search Algorithm for Dependent Design Spaces (Version 2.0.0) [Computer software]. https://doi.org/10.5281/zenodo.19160019

BibTeX:

@software{ozcan_hsds_2026,
  author       = {Özcan, Abdulkadir},
  title        = {hsds: Harmony Search Algorithm for Dependent Design Spaces},
  month        = mar,
  year         = 2026,
  publisher    = {Zenodo},
  version      = {2.0.0},
  doi          = {10.5281/zenodo.19160019},
  url          = {https://doi.org/10.5281/zenodo.19160019}
}

---

License

MIT © Abdulkadir Özcan