ODFstudio/odfstudio
2
0
Fork 0
Code Issues Packages Activity

SwarmSolver: fix PSO position and velocity handling

Browse Source 
Three problems in the particle swarm code:

The best-position vectors were updated with insert at the front rather
than assign, so each improvement prepended another copy. The vectors
grew on every step and the returned solution ended up longer than the
number of variables.

The velocity was seeded from initializeVariables, which fills it with a
random position. The first step then moved each particle by a whole
coordinate. The velocity now starts at zero.

The result of clamping the velocity was thrown away, so the call did
nothing. Clamping the velocity to the position bounds would in fact be
wrong when the bounds do not straddle zero, so the dead call is removed
and a comment explains why the velocity is left alone.

Also drops a redundant second clamp of the position during setup.

Tests drive the algorithm through a small mock data provider so the
checks do not depend on the random seed: one asserts the result has one
value per variable, the other asserts setup calls initializeVariables
once per particle. The large-objective test now runs Particle Swarm
Optimization too.

Signed-off-by: Tomaž Vajngerl <tomaz.vajngerl@collabora.co.uk>
Change-Id: I8197aa3368d20108c06bd3643ef1138ac6f2254a
Reviewed-on: https://gerrit.collaboraoffice.com/c/online/+/4050
Tested-by: Tomaž Vajngerl <tomaz.vajngerl@collabora.com>
Reviewed-by: Miklos Vajna <vmiklos@collabora.com>
(cherry picked from commit 9353efb6419c16f0bcb104e57c4eb08ddcb4da81)
Reviewed-on: https://gerrit.libreoffice.org/c/core/+/208435
Tested-by: Xisco Fauli <xiscofauli@libreoffice.org>
Reviewed-by: Xisco Fauli <xiscofauli@libreoffice.org>
This commit is contained in:
Tomaž Vajngerl
2026-07-22 21:51:22 +02:00
committed by Xisco Fauli
parent 1b87464d79
commit 67b9784801
3 changed files with 154 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files
sc/CppunitTest_sc_swarmsolvertest.mk
+1
View File
@@ -57,6 +57,7 @@ $(eval $(call gb_CppunitTest_use_libraries,sc_swarmsolvertest,\
$(eval $(call gb_CppunitTest_set_include,sc_swarmsolvertest,\
-I$(SRCDIR)/sc/inc \
-I$(SRCDIR)/sc/source/solver \
$$(INCLUDE) \
))
sc/qa/unit/SwarmSolverTest.cxx
+143 -8
View File
@@ -16,10 +16,49 @@
#include <test/unoapi_test.hxx>
#include <ParticelSwarmOptimization.hxx>
using namespace css;
namespace
{
// Minimal data provider that drives the algorithm template without a document.
class MockProvider
{
size_t mnDimensionality;
public:
int mnInitCalls = 0;
explicit MockProvider(size_t nDimensionality)
: mnDimensionality(nDimensionality)
{
}
size_t getDimensionality() const { return mnDimensionality; }
void initializeVariables(std::vector<double>& rVariables, std::mt19937& rGenerator)
{
std::uniform_real_distribution<double> aRandom(-10.0, 10.0);
rVariables.resize(mnDimensionality);
for (double& rValue : rVariables)
rValue = aRandom(rGenerator);
mnInitCalls++;
}
// Separable, with its maximum where every variable is zero, so the swarm
// always has somewhere to move towards and improvements actually happen.
static double calculateFitness(std::vector<double> const& rVariables)
{
double fFitness = 0.0;
for (double fValue : rVariables)
fFitness -= fValue * fValue;
return fFitness;
}
static double clampVariable(size_t, double fValue) { return fValue; }
};
class SwarmSolverTest : public UnoApiTest
{
void testUnconstrained();
@@ -29,7 +68,10 @@ class SwarmSolverTest : public UnoApiTest
void testMultipleVariables();
void testInfeasibleConstraints_Differential_Evolution();
void testInfeasibleConstraints_Particle_Swarm();
void testLargeObjectiveStillSolvable();
void testLargeObjectiveStillSolvable_Differential_Evolution();
void testLargeObjectiveStillSolvable_Particle_Swarm();
void testParticleSwarmResultLength();
void testParticleSwarmVelocityNotInitializedAsPosition();
public:
SwarmSolverTest()
@@ -45,7 +87,10 @@ public:
CPPUNIT_TEST(testTwoVariables);
CPPUNIT_TEST(testInfeasibleConstraints_Differential_Evolution);
CPPUNIT_TEST(testInfeasibleConstraints_Particle_Swarm);
CPPUNIT_TEST(testLargeObjectiveStillSolvable);
CPPUNIT_TEST(testLargeObjectiveStillSolvable_Differential_Evolution);
CPPUNIT_TEST(testLargeObjectiveStillSolvable_Particle_Swarm);
CPPUNIT_TEST(testParticleSwarmResultLength);
CPPUNIT_TEST(testParticleSwarmVelocityNotInitializedAsPosition);
CPPUNIT_TEST_SUITE_END();
};
@@ -417,11 +462,12 @@ void SwarmSolverTest::testInfeasibleConstraints_Particle_Swarm()
CPPUNIT_ASSERT_MESSAGE("Infeasible model must not report success", !xSolver->getSuccess());
}
void SwarmSolverTest::testLargeObjectiveStillSolvable()
void SwarmSolverTest::testLargeObjectiveStillSolvable_Differential_Evolution()
{
// Regression: a feasible objective below the float range used to lose to the
// infeasible penalty, so the search gave up on a solvable model. Differential
// Evolution only, until Particle Swarm Optimization's result handling is fixed.
// infeasible penalty, so the search gave up on a solvable model. The Particle
// Swarm result also used to grow past the variable count, which the length
// check below guards.
loadFromFile(u"Simple.ods");
@@ -432,14 +478,17 @@ void SwarmSolverTest::testLargeObjectiveStillSolvable()
u"com.sun.star.comp.Calc.SwarmSolver"_ustr, m_xContext),
uno::UNO_QUERY_THROW);
uno::Reference<beans::XPropertySet> xPropSet(xSolver, uno::UNO_QUERY_THROW);
xPropSet->setPropertyValue(u"Algorithm"_ustr, uno::Any(sal_Int32(0)));
table::CellAddress aObjective(0, 1, 1);
uno::Sequence<table::CellAddress> aVariables{ { 0, 1, 0 } };
// The variable is pushed far enough out that the objective
// 10*B1^2 - 60*B1 - 40 reaches the 1e39 range, beyond what a float can
// hold. Minimizing turns that into a fitness near -1e39. The objective cap
// is a non-bounded constraint, so the upper part of the box is infeasible
// while the lower end stays feasible.
// hold. Minimizing turns that into a fitness near -1e39. The objective
// cap is a non-bounded constraint, so the upper part of the box is
// infeasible while the lower end stays feasible.
uno::Sequence<sheet::SolverConstraint> aConstraints{
{ /* [0] Left */ table::CellAddress(0, 1, 0),
/* Operator */ sheet::SolverConstraintOperator_GREATER_EQUAL,
@@ -469,6 +518,92 @@ void SwarmSolverTest::testLargeObjectiveStillSolvable()
CPPUNIT_ASSERT(aSolution[0] <= 2.0e19);
}
void SwarmSolverTest::testLargeObjectiveStillSolvable_Particle_Swarm()
{
// Regression: a feasible objective below the float range used to lose to the
// infeasible penalty, so the search gave up on a solvable model. The Particle
// Swarm result also used to grow past the variable count, which the length
// check below guards.
loadFromFile(u"Simple.ods");
uno::Reference<sheet::XSpreadsheetDocument> xDocument(mxComponent, uno::UNO_QUERY_THROW);
uno::Reference<sheet::XSolver> xSolver;
xSolver.set(m_xContext->getServiceManager()->createInstanceWithContext(
u"com.sun.star.comp.Calc.SwarmSolver"_ustr, m_xContext),
uno::UNO_QUERY_THROW);
uno::Reference<beans::XPropertySet> xPropSet(xSolver, uno::UNO_QUERY_THROW);
xPropSet->setPropertyValue(u"Algorithm"_ustr, uno::Any(sal_Int32(1)));
table::CellAddress aObjective(0, 1, 1);
uno::Sequence<table::CellAddress> aVariables{ { 0, 1, 0 } };
// The variable is pushed far enough out that the objective
// 10*B1^2 - 60*B1 - 40 reaches the 1e39 range, beyond what a float can
// hold. Minimizing turns that into a fitness near -1e39. The objective
// cap is a non-bounded constraint, so the upper part of the box is
// infeasible while the lower end stays feasible.
uno::Sequence<sheet::SolverConstraint> aConstraints{
{ /* [0] Left */ table::CellAddress(0, 1, 0),
/* Operator */ sheet::SolverConstraintOperator_GREATER_EQUAL,
/* Right */ uno::Any(1.0e19) },
{ /* [1] Left */ table::CellAddress(0, 1, 0),
/* Operator */ sheet::SolverConstraintOperator_LESS_EQUAL,
/* Right */ uno::Any(2.0e19) },
{ /* [2] Left */ table::CellAddress(0, 1, 1),
/* Operator */ sheet::SolverConstraintOperator_LESS_EQUAL,
/* Right */ uno::Any(2.0e39) }
};
xSolver->setDocument(xDocument);
xSolver->setObjective(aObjective);
xSolver->setVariables(aVariables);
xSolver->setConstraints(aConstraints);
xSolver->setMaximize(false);
xSolver->solve();
CPPUNIT_ASSERT_MESSAGE("Solvable model must report success", xSolver->getSuccess());
uno::Sequence<double> aSolution = xSolver->getSolution();
CPPUNIT_ASSERT_EQUAL(aVariables.getLength(), aSolution.getLength());
// the returned point must lie in the feasible part of the box
CPPUNIT_ASSERT(aSolution[0] >= 1.0e19);
CPPUNIT_ASSERT(aSolution[0] <= 2.0e19);
}
void SwarmSolverTest::testParticleSwarmResultLength()
{
// Regression: the best-position vectors were updated with insert at the front
// instead of assign, so every improvement prepended another copy and the result
// grew past the variable count.
MockProvider aProvider(3);
ParticleSwarmOptimizationAlgorithm<MockProvider> aAlgorithm(aProvider, 8);
aAlgorithm.initialize();
for (int i = 0; i < 20; ++i)
aAlgorithm.next();
CPPUNIT_ASSERT_EQUAL(size_t(3), aAlgorithm.getResult().size());
}
void SwarmSolverTest::testParticleSwarmVelocityNotInitializedAsPosition()
{
// Regression: the velocity was seeded by initializeVariables, which fills it
// with a random position. Setup must call initializeVariables once per particle
// (for the position only), not twice.
MockProvider aProvider(3);
ParticleSwarmOptimizationAlgorithm<MockProvider> aAlgorithm(aProvider, 8);
aAlgorithm.initialize();
CPPUNIT_ASSERT_EQUAL(8, aProvider.mnInitCalls);
}
CPPUNIT_TEST_SUITE_REGISTRATION(SwarmSolverTest);
}
sc/source/solver/ParticelSwarmOptimization.hxx
+10 -16
View File
@@ -96,7 +96,9 @@ public:
Particle& rParticle = maSwarm.back();
mrDataProvider.initializeVariables(rParticle.mPosition, maGenerator);
mrDataProvider.initializeVariables(rParticle.mVelocity, maGenerator);
// Only the position is seeded. The velocity keeps the zero the
// constructor gave it, so a particle starts at rest and builds up
// speed gradually over the first few steps.
for (size_t k = 0; k < mnDimensionality; k++)
{
@@ -105,20 +107,13 @@ public:
rParticle.mCurrentFitness = mrDataProvider.calculateFitness(rParticle.mPosition);
for (size_t k = 0; k < mnDimensionality; k++)
{
rParticle.mPosition[k] = mrDataProvider.clampVariable(k, rParticle.mPosition[k]);
}
rParticle.mBestPosition.insert(rParticle.mBestPosition.begin(),
rParticle.mPosition.begin(), rParticle.mPosition.end());
rParticle.mBestPosition.assign(rParticle.mPosition.begin(), rParticle.mPosition.end());
rParticle.mBestFitness = rParticle.mCurrentFitness;
if (rParticle.mCurrentFitness > mfBestFitness)
{
mfBestFitness = rParticle.mCurrentFitness;
maBestPosition.insert(maBestPosition.begin(), rParticle.mPosition.begin(),
rParticle.mPosition.end());
maBestPosition.assign(rParticle.mPosition.begin(), rParticle.mPosition.end());
}
}
}
@@ -139,8 +134,9 @@ public:
+ (c1 * fRandom1 * (rParticle.mBestPosition[k] - rParticle.mPosition[k]))
+ (c2 * fRandom2 * (maBestPosition[k] - rParticle.mPosition[k]));
mrDataProvider.clampVariable(k, rParticle.mVelocity[k]);
// The inertia weight below one keeps the velocity bounded on
// its own, so only the position is clamped, which holds the
// particle inside its bounds.
rParticle.mPosition[k] += rParticle.mVelocity[k];
rParticle.mPosition[k] = mrDataProvider.clampVariable(k, rParticle.mPosition[k]);
}
@@ -150,8 +146,7 @@ public:
if (rParticle.mCurrentFitness > rParticle.mBestFitness)
{
rParticle.mBestFitness = rParticle.mCurrentFitness;
rParticle.mBestPosition.insert(rParticle.mBestPosition.begin(),
rParticle.mPosition.begin(),
rParticle.mBestPosition.assign(rParticle.mPosition.begin(),
rParticle.mPosition.end());
}
@@ -162,8 +157,7 @@ public:
bBestChanged = true;
mnLastChange = mnGeneration;
}
maBestPosition.insert(maBestPosition.begin(), rParticle.mPosition.begin(),
rParticle.mPosition.end());
maBestPosition.assign(rParticle.mPosition.begin(), rParticle.mPosition.end());
mfBestFitness = rParticle.mCurrentFitness;
}
}