Swarm based (uses PSO or DE) experimental non-linear solver

This is a new, simple non-linear solver that uses a swarm
(population) to do global optimization. It uses two algoritms -
Particle Swarm Optimization (PSO) or Differential Evolution (DE)
to find a (non-optimal) solution.

It is experimental as not all functions are implemented and it
needs a lot more testing so that it performs well.

Change-Id: If55dad7eda17394851a9d178ad892de771eca7c9
Reviewed-on: https://gerrit.libreoffice.org/44382Tested-by: Jenkins <ci@libreoffice.org>
Reviewed-by: Tomaž Vajngerl <quikee@gmail.com>

postprocess/Rdb_services.mk

@@ -118,6 +118,7 @@ $(eval $(call gb_Rdb_add_components,services,\
 	$(if $(ENABLE_LPSOLVE), \
 		sccomp/source/solver/lpsolvesolver \
 	) \
+	sccomp/source/solver/swarmsolver \
 	writerfilter/util/writerfilter \
 	writerperfect/source/draw/wpftdraw \
 	writerperfect/source/impress/wpftimpress \

sccomp/CppunitTest_sccomp_swarmsolvertest.mk

+# -*- Mode: makefile-gmake; tab-width: 4; indent-tabs-mode: t -*-
+#
+# This file is part of the LibreOffice project.
+#
+# This Source Code Form is subject to the terms of the Mozilla Public
+# License, v. 2.0. If a copy of the MPL was not distributed with this
+# file, You can obtain one at http://mozilla.org/MPL/2.0/.
+#
+
+$(eval $(call gb_CppunitTest_CppunitTest,swarm_solver_test))
+
+$(eval $(call gb_CppunitTest_add_exception_objects,swarm_solver_test,\
+	sccomp/qa/unit/SwarmSolverTest \
+))
+
+$(eval $(call gb_CppunitTest_use_externals,swarm_solver_test,\
+	boost_headers \
+))
+
+$(eval $(call gb_CppunitTest_use_libraries,swarm_solver_test,\
+	basegfx \
+	comphelper \
+	cppu \
+	cppuhelper \
+	drawinglayer \
+	editeng \
+	for \
+	forui \
+	i18nlangtag \
+	msfilter \
+	oox \
+	sal \
+	salhelper \
+	sax \
+	sb \
+	sc \
+	scqahelper \
+	sfx \
+	sot \
+	subsequenttest \
+	svl \
+	svt \
+	svx \
+	svxcore \
+	test \
+	tk \
+	tl \
+	ucbhelper \
+	unotest \
+	utl \
+	vbahelper \
+	vcl \
+	xo \
+	$(gb_UWINAPI) \
+))
+
+$(eval $(call gb_CppunitTest_set_include,swarm_solver_test,\
+	-I$(SRCDIR)/sc/inc \
+	$$(INCLUDE) \
+))
+
+$(eval $(call gb_CppunitTest_use_sdk_api,swarm_solver_test))
+
+$(eval $(call gb_CppunitTest_use_ure,swarm_solver_test))
+$(eval $(call gb_CppunitTest_use_vcl,swarm_solver_test))
+
+$(eval $(call gb_CppunitTest_use_rdb,swarm_solver_test,services))
+
+$(eval $(call gb_CppunitTest_use_configuration,swarm_solver_test))
+
+# vim: set noet sw=4 ts=4:

sccomp/Library_solver.mk

@@ -22,6 +22,8 @@ $(eval $(call gb_Library_Library,solver))
 $(if $(ENABLE_COINMP),$(eval $(call gb_Library_set_componentfile,solver,sccomp/source/solver/coinmpsolver)))
 $(if $(ENABLE_LPSOLVE),$(eval $(call gb_Library_set_componentfile,solver,sccomp/source/solver/lpsolvesolver)))
 
+$(eval $(call gb_Library_set_componentfile,solver,sccomp/source/solver/swarmsolver))
+
 $(eval $(call gb_Library_use_sdk_api,solver))
 
 $(eval $(call gb_Library_set_include,solver,\
@@ -45,6 +47,7 @@ $(eval $(call gb_Library_use_externals,solver,\
 ))
 
 $(eval $(call gb_Library_add_exception_objects,solver,\
+	sccomp/source/solver/SwarmSolver \
 	sccomp/source/solver/SolverComponent \
 	$(if $(ENABLE_COINMP), sccomp/source/solver/CoinMPSolver) \
 	$(if $(ENABLE_LPSOLVE), sccomp/source/solver/LpsolveSolver) \

sccomp/Module_sccomp.mk

@@ -29,6 +29,7 @@ $(eval $(call gb_Module_add_l10n_targets,sccomp,\
 
 $(eval $(call gb_Module_add_check_targets,sccomp,\
 	CppunitTest_sccomp_solver \
+	CppunitTest_sccomp_swarmsolvertest \
 ))
 
 # vim: set noet sw=4 ts=4:

sccomp/inc/strings.hrc

@@ -24,11 +24,13 @@
 
 #define RID_SOLVER_COMPONENT        NC_("RID_SOLVER_COMPONENT", "%PRODUCTNAME Linear Solver")
 #define RID_COINMP_SOLVER_COMPONENT NC_("RID_COINMP_SOLVER_COMPONENT", "%PRODUCTNAME CoinMP Linear Solver")
+#define RID_SWARM_SOLVER_COMPONENT  NC_("RID_SWARM_SOLVER_COMPONENT", "%PRODUCTNAME Swarm Non-Linear Solver (experimental)")
 #define RID_PROPERTY_NONNEGATIVE    NC_("RID_PROPERTY_NONNEGATIVE", "Assume variables as non-negative")
 #define RID_PROPERTY_INTEGER        NC_("RID_PROPERTY_INTEGER", "Assume variables as integer")
 #define RID_PROPERTY_TIMEOUT        NC_("RID_PROPERTY_TIMEOUT", "Solving time limit (seconds)")
 #define RID_PROPERTY_EPSILONLEVEL   NC_("RID_PROPERTY_EPSILONLEVEL", "Epsilon level (0-3)")
 #define RID_PROPERTY_LIMITBBDEPTH   NC_("RID_PROPERTY_LIMITBBDEPTH", "Limit branch-and-bound depth")
+#define RID_PROPERTY_ALGORITHM      NC_("RID_PROPERTY_ALGORITHM", "Swarm algorithm (0 - Differential Evolution, 1 - Particle Swarm Optimization)")
 #define RID_ERROR_NONLINEAR         NC_("RID_ERROR_NONLINEAR", "The model is not linear.")
 #define RID_ERROR_EPSILONLEVEL      NC_("RID_ERROR_EPSILONLEVEL", "The epsilon level is invalid.")
 #define RID_ERROR_INFEASIBLE        NC_("RID_ERROR_INFEASIBLE", "The model is infeasible. Check limiting conditions.")

sccomp/source/solver/DifferentialEvolution.hxx

+/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+/*
+ * This file is part of the LibreOffice project.
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ *
+ */
+
+#ifndef INCLUDED_SCCOMP_SOURCE_DIFFERENTIALEVOLUTION_HXX
+#define INCLUDED_SCCOMP_SOURCE_DIFFERENTIALEVOLUTION_HXX
+
+#include <vector>
+#include <random>
+#include <limits>
+
+struct Individual
+{
+    std::vector<double> mVariables;
+};
+
+template <typename DataProvider> class DifferentialEvolutionAlgorithm
+{
+    static constexpr double mnDifferentialWeight = 0.5; // [0, 2]
+    static constexpr double mnCrossoverProbability = 0.9; // [0, 1]
+
+    static constexpr double constAcceptedPrecision = 0.000000001;
+
+    DataProvider& mrDataProvider;
+
+    size_t mnPopulationSize;
+    std::vector<Individual> maPopulation;
+
+    std::random_device maRandomDevice;
+    std::mt19937 maGenerator;
+    size_t mnDimensionality;
+
+    std::uniform_int_distribution<> maRandomPopulation;
+    std::uniform_int_distribution<> maRandomDimensionality;
+    std::uniform_real_distribution<> maRandom01;
+
+    Individual maBestCandidate;
+    double mfBestFitness;
+    int mnGeneration;
+    int mnLastChange;
+
+public:
+    DifferentialEvolutionAlgorithm(DataProvider& rDataProvider, size_t nPopulationSize)
+        : mrDataProvider(rDataProvider)
+        , mnPopulationSize(nPopulationSize)
+        , maGenerator(maRandomDevice())
+        , mnDimensionality(mrDataProvider.getDimensionality())
+        , maRandomPopulation(0, mnPopulationSize - 1)
+        , maRandomDimensionality(0, mnDimensionality - 1)
+        , maRandom01(0.0, 1.0)
+        , mfBestFitness(std::numeric_limits<double>::lowest())
+        , mnGeneration(0)
+        , mnLastChange(0)
+    {
+    }
+
+    std::vector<double> const& getResult() { return maBestCandidate.mVariables; }
+
+    int getGeneration() { return mnGeneration; }
+
+    int getLastChange() { return mnLastChange; }
+
+    void initialize()
+    {
+        mnGeneration = 0;
+        mnLastChange = 0;
+        maPopulation.clear();
+        maBestCandidate.mVariables.clear();
+
+        // Initialize population with individuals that have been initialized with uniform random
+        // noise
+        // uniform noise means random value inside your search space
+        for (size_t i = 0; i < mnPopulationSize; ++i)
+        {
+            maPopulation.emplace_back();
+            Individual& rIndividual = maPopulation.back();
+            mrDataProvider.initializeVariables(rIndividual.mVariables, maGenerator);
+        }
+    }
+
+    // Calculate one generation
+    bool next()
+    {
+        bool bBestChanged = false;
+
+        for (size_t agentIndex = 0; agentIndex < mnPopulationSize; ++agentIndex)
+        {
+            // calculate new candidate solution
+
+            // pick random point from population
+            size_t x = agentIndex; // randomPopulation(generator);
+            size_t a, b, c;
+
+            // create a copy of choosen random agent in population
+            Individual& rOriginal = maPopulation[x];
+            Individual aCandidate(rOriginal);
+
+            // pick three different random points from population
+            do
+            {
+                a = maRandomPopulation(maGenerator);
+            } while (a == x);
+
+            do
+            {
+                b = maRandomPopulation(maGenerator);
+            } while (b == x || b == a);
+
+            do
+            {
+                c = maRandomPopulation(maGenerator);
+
+            } while (c == x || c == a || c == b);
+
+            size_t randomIndex = maRandomDimensionality(maGenerator);
+
+            for (size_t index = 0; index < mnDimensionality; ++index)
+            {
+                double randomCrossoverProbability = maRandom01(maGenerator);
+                if (index == randomIndex || randomCrossoverProbability < mnCrossoverProbability)
+                {
+                    double fVarA = maPopulation[a].mVariables[index];
+                    double fVarB = maPopulation[b].mVariables[index];
+                    double fVarC = maPopulation[c].mVariables[index];
+
+                    double fNewValue = fVarA + mnDifferentialWeight * (fVarB - fVarC);
+                    fNewValue = mrDataProvider.boundVariable(index, fNewValue);
+                    aCandidate.mVariables[index] = fNewValue;
+                }
+            }
+
+            double fCandidateFitness = mrDataProvider.calculateFitness(aCandidate.mVariables);
+
+            // see if is better than original, if so replace
+            if (fCandidateFitness > mrDataProvider.calculateFitness(rOriginal.mVariables))
+            {
+                maPopulation[x] = aCandidate;
+
+                if (fCandidateFitness > mfBestFitness)
+                {
+                    if (std::abs(fCandidateFitness - mfBestFitness) > constAcceptedPrecision)
+                    {
+                        bBestChanged = true;
+                        mnLastChange = mnGeneration;
+                    }
+                    mfBestFitness = fCandidateFitness;
+                    maBestCandidate = maPopulation[x];
+                }
+            }
+        }
+        mnGeneration++;
+        return bBestChanged;
+    }
+};
+
+#endif
+
+/* vim:set shiftwidth=4 softtabstop=4 expandtab: */

sccomp/source/solver/ParticelSwarmOptimization.hxx

+/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+/*
+ * This file is part of the LibreOffice project.
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ *
+ */
+
+#ifndef INCLUDED_SCCOMP_SOURCE_PARTICLESWARM_HXX
+#define INCLUDED_SCCOMP_SOURCE_PARTICLESWARM_HXX
+
+#include <vector>
+#include <random>
+#include <limits>
+
+struct Particle
+{
+    Particle(size_t nDimensionality)
+        : mVelocity(nDimensionality)
+        , mPosition(nDimensionality)
+        , mCurrentFitness(std::numeric_limits<float>::lowest())
+        , mBestPosition(nDimensionality)
+        , mBestFitness(std::numeric_limits<float>::lowest())
+    {
+    }
+
+    std::vector<double> mVelocity;
+
+    std::vector<double> mPosition;
+    double mCurrentFitness;
+
+    std::vector<double> mBestPosition;
+    double mBestFitness;
+};
+
+template <typename DataProvider> class ParticleSwarmOptimizationAlgorithm
+{
+private:
+    // inertia
+    static constexpr double constWeight = 0.729;
+    // cognitive coefficient
+    static constexpr double c1 = 1.49445;
+    // social  coefficient
+    static constexpr double c2 = 1.49445;
+
+    static constexpr double constAcceptedPrecision = 0.000000001;
+
+    DataProvider& mrDataProvider;
+
+    size_t mnNumOfParticles;
+
+    std::vector<Particle> maSwarm;
+
+    std::random_device maRandomDevice;
+    std::mt19937 maGenerator;
+    size_t mnDimensionality;
+
+    std::uniform_real_distribution<> maRandom01;
+
+    std::vector<double> maBestPosition;
+    double mfBestFitness;
+    int mnGeneration;
+    int mnLastChange;
+
+public:
+    ParticleSwarmOptimizationAlgorithm(DataProvider& rDataProvider, size_t nNumOfParticles)
+        : mrDataProvider(rDataProvider)
+        , mnNumOfParticles(nNumOfParticles)
+        , maGenerator(maRandomDevice())
+        , mnDimensionality(mrDataProvider.getDimensionality())
+        , maRandom01(0.0, 1.0)
+        , maBestPosition(mnDimensionality)
+        , mfBestFitness(std::numeric_limits<float>::lowest())
+        , mnGeneration(0)
+        , mnLastChange(0)
+    {
+    }
+
+    std::vector<double> const& getResult() { return maBestPosition; }
+
+    int getGeneration() { return mnGeneration; }
+
+    int getLastChange() { return mnLastChange; }
+
+    void initialize()
+    {
+        mnGeneration = 0;
+        mnLastChange = 0;
+        maSwarm.clear();
+
+        mfBestFitness = std::numeric_limits<float>::lowest();
+
+        for (size_t i = 0; i < mnNumOfParticles; i++)
+        {
+            maSwarm.emplace_back(mnDimensionality);
+            Particle& rParticle = maSwarm.back();
+
+            mrDataProvider.initializeVariables(rParticle.mPosition, maGenerator);
+            mrDataProvider.initializeVariables(rParticle.mVelocity, maGenerator);
+
+            for (size_t k = 0; k < mnDimensionality; k++)
+            {
+                rParticle.mPosition[k] = mrDataProvider.clampVariable(k, rParticle.mPosition[k]);
+            }
+
+            rParticle.mCurrentFitness = mrDataProvider.calculateFitness(rParticle.mPosition);
+
+            for (size_t k = 0; k < mnDimensionality; k++)
+            {
+                rParticle.mPosition[k] = mrDataProvider.clampVariable(k, rParticle.mPosition[k]);
+            }
+
+            std::copy(rParticle.mPosition.begin(), rParticle.mPosition.end(),
+                      rParticle.mBestPosition.begin());
+            rParticle.mBestFitness = rParticle.mCurrentFitness;
+
+            if (rParticle.mCurrentFitness > mfBestFitness)
+            {
+                mfBestFitness = rParticle.mCurrentFitness;
+                std::copy(rParticle.mPosition.begin(), rParticle.mPosition.end(),
+                          maBestPosition.begin());
+            }
+        }
+    }
+
+    bool next()
+    {
+        bool bBestChanged = false;
+
+        for (Particle& rParticle : maSwarm)
+        {
+            double fRandom1 = maRandom01(maGenerator);
+            double fRandom2 = maRandom01(maGenerator);
+
+            for (size_t k = 0; k < mnDimensionality; k++)
+            {
+                rParticle.mVelocity[k]
+                    = (constWeight * rParticle.mVelocity[k])
+                      + (c1 * fRandom1 * (rParticle.mBestPosition[k] - rParticle.mPosition[k]))
+                      + (c2 * fRandom2 * (maBestPosition[k] - rParticle.mPosition[k]));
+
+                mrDataProvider.clampVariable(k, rParticle.mVelocity[k]);
+
+                rParticle.mPosition[k] += rParticle.mVelocity[k];
+                rParticle.mPosition[k] = mrDataProvider.clampVariable(k, rParticle.mPosition[k]);
+            }
+
+            rParticle.mCurrentFitness = mrDataProvider.calculateFitness(rParticle.mPosition);
+
+            if (rParticle.mCurrentFitness > rParticle.mBestFitness)
+            {
+                rParticle.mBestFitness = rParticle.mCurrentFitness;
+                std::copy(rParticle.mPosition.begin(), rParticle.mPosition.end(),
+                          rParticle.mBestPosition.begin());
+            }
+
+            if (rParticle.mCurrentFitness > mfBestFitness)
+            {
+                if (std::abs(rParticle.mCurrentFitness - mfBestFitness) > constAcceptedPrecision)
+                {
+                    bBestChanged = true;
+                    mnLastChange = mnGeneration;
+                }
+                std::copy(rParticle.mPosition.begin(), rParticle.mPosition.end(),
+                          maBestPosition.begin());
+                mfBestFitness = rParticle.mCurrentFitness;
+            }
+        }
+        mnGeneration++;
+        return bBestChanged;
+    }
+};
+
+#endif
+
+/* vim:set shiftwidth=4 softtabstop=4 expandtab: */

sccomp/source/solver/swarmsolver.component

+<?xml version="1.0" encoding="UTF-8"?>
+<!--
+ * This file is part of the LibreOffice project.
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ -->
+
+<component loader="com.sun.star.loader.SharedLibrary" environment="@CPPU_ENV@"
+    xmlns="http://openoffice.org/2010/uno-components">
+  <implementation name="com.sun.star.comp.Calc.SwarmSolver" constructor="com_sun_star_comp_Calc_SwarmSolver_get_implementation">
+    <service name="com.sun.star.sheet.Solver"/>
+  </implementation>
+</component>