Variator.h

Go to the documentation of this file.

//
// Class Variator
//
// Copyright (c) 2010 - 2013, Yves Ineichen, ETH Zürich
// All rights reserved
//
// Implemented as part of the PhD thesis
// "Toward massively parallel multi-objective optimization with application to
// particle accelerators" (https://doi.org/10.3929/ethz-a-009792359)
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#ifndef __VARIATOR_H__
#define __VARIATOR_H__

#include <string>
#include <vector>
#include <map>
#include <utility>

#include "Util/Types.h"
#include "Util/CmdArguments.h"
#include "Optimizer/EA/Population.h"
#include "Optimizer/Optimizer.h"

#include "boost/property_tree/ptree.hpp"
#include "boost/property_tree/json_parser.hpp"

template<
      class ind_t
    , template <class> class CrossoverOperator
    , template <class> class MutationOperator
>
class Variator : public CrossoverOperator<ind_t>,
                 public MutationOperator<ind_t>
{

public:

    Variator(std::vector<std::string> dNames,
             Optimizer::bounds_t dVarBounds, Expressions::Named_t constraints,
             CmdArguments_t args)
        : dNames_m(dNames)
        , dVarBounds_m(dVarBounds)
    {
        // add constraints, if only design variables are needed for evaluation
        for(auto constraint : constraints) {
            bool allDesignVariables = true;
            std::set<std::string> req_vars = constraint.second->getReqVars();
            if (req_vars.empty()) allDesignVariables = false;
            for (std::string req_var : req_vars) {
                // check if it is a design variable
                if (std::find(dNames_m.begin(),dNames_m.end(),req_var) == dNames_m.end()) {
                    allDesignVariables = false;
                    break;
                }
            }
            if (allDesignVariables == true)
                constraints_m.insert(constraint);
        }

        //FIXME: pass population as arg to variator
        //std::shared_ptr< Population<ind_t> >
        population_m.reset(new Population<ind_t>());

        mutationProbability_m =
            args->getArg<double>("mutation-probability", 0.5);

        recombinationProbability_m =
            args->getArg<double>("recombination-probability", 0.5);

        args_ = args;
    }

    ~Variator() {
    }

    //FIXME access population from outside
    std::shared_ptr< Population<ind_t> > population() {
        return population_m;
    }

    void initial_population(size_t sizeInitial, std::string fname) {
        if ( fname.empty() ) {
            for(size_t i = 0; i < sizeInitial; i++) {
                new_individual();
            }
            return;
        }
        // read population from file
        typedef boost::property_tree::ptree ptree_t;
        ptree_t tree;

        boost::property_tree::read_json(fname, tree);

        if ( tree.get<std::string>("name").compare("opt-pilot") ) {
            throw OptPilotException("Variator::initial_population",
                                    "Json file not generated by opt-pilot.");
        }

        std::size_t nDVars = 0;
        for(auto& elem : tree.get_child("dvar-bounds")) {
            // check if it is a design variable
            if (std::find(dNames_m.begin(), dNames_m.end(), elem.first) == dNames_m.end()) {
                throw OptPilotException("Variator::initial_population",
                                        "The design variable '" + elem.first + "' is not in the list.");
            }
            ++nDVars;
        }

        if ( nDVars != dVarBounds_m.size() ) {
            throw OptPilotException("Variator::initial_population",
                                    "The number of design variables do not agree.");
        }

        boost::property_tree::ptree& population = tree.get_child("population");

        std::size_t size = 0;
        Individual::genes_t dvars(nDVars);
        for (auto& ind : population ) {

            if ( size > sizeInitial - 1 )
                break;

            /* just to be sure: It might be that the order of the design variables in the
             * Json file is not the same as it reads in, therefore, we check and take the order
             * of the container dNames_m.
             */
            for (auto& dvar : population.get_child(ind.first).get_child("dvar")) {
                auto it = std::find(dNames_m.begin(), dNames_m.end(), dvar.first);
                std::size_t idx = std::distance(dNames_m.begin(), it);
                //FIXME requires random access of Individual::genes_t
                dvars[idx] = dvar.second.get_value<double>();
            }

            new_individual(dvars);
            ++size;
        }

        // fill with random individuals
        for (std::size_t i = size; i < sizeInitial; ++i) {
            new_individual();
        }
    }

    void infeasible(std::shared_ptr<ind_t> ind) {
        population_m->remove_individual(ind);
        new_individual();
    }

    bool hasMoreIndividualsToEvaluate() {
        return !individualsToEvaluate_m.empty();
    }

    std::shared_ptr<ind_t> popIndividualToEvaluate() {
        unsigned int ind = individualsToEvaluate_m.front();
        individualsToEvaluate_m.pop();
        return population_m->get_staging(ind);
    }

    void variate(std::vector<unsigned int> parents) {

        // copying all individuals from parents
        for(unsigned int parent : parents) {
            new_individual( population_m->get_individual(parent) );
        }

        // only variate new offspring, individuals in staging area have been
        // variated already
        std::queue<unsigned int> tmp(individualsToEvaluate_m);
        while(!tmp.empty()) {

            // pop first individual
            unsigned int idx = tmp.front(); tmp.pop();
            std::shared_ptr<ind_t> a = population_m->get_staging(idx);

            // handle special case where we have an odd number of offspring
            if(tmp.empty()) {
                if (drand(1) <= mutationProbability_m) {
                    // temporary copy in case not successful
                    std::shared_ptr<ind_t> copyA(new ind_t(a));
                    int iter = 0;
                    while (true) {
                        // assign with shared pointer constructor
                        *a = copyA;
                        this->mutate(a, args_);
                        // check if viable offspring
                        if (a->viable()) break;

                        iter++;
                        // if maximum number of tries then create new individual
                        if (iter > 100) {
                            infeasible(a);
                            break;
                        }
                    }
                }
                break;
            }

            // and second if any
            idx = tmp.front(); tmp.pop();
            std::shared_ptr<ind_t> b = population_m->get_staging(idx);

            // create new individuals

            // temporary copy in case not successful
            std::shared_ptr<ind_t> copyA(new ind_t(a));
            std::shared_ptr<ind_t> copyB(new ind_t(b));

            int iter = 0;
            while (true) {
                // assign with shared pointer constructor
                *a = copyA;
                *b = copyB;

                // do recombination
                if(drand(1) <= recombinationProbability_m) {
                    this->crossover(a, b, args_);
                }

                // do mutation
                if (drand(1) <= mutationProbability_m) {
                    this->mutate(a, args_);
                    this->mutate(b, args_);
                }

                // check if viable offspring
                bool viableA = a->viable();
                bool viableB = b->viable();
                if (viableA == true && viableB == true) {
                    break;
                }
                std::cout << "Individuals not viable, I try again: iter= " << iter << std::endl;
                iter++;
                // if maximum number of tries then create new individual(s)
                if (iter > 100) {
                    if (viableA == false) {
                        infeasible(a);
                    }
                    if (viableB == false) {
                        infeasible(b);
                    }
                    break;
                }
            }
        }
    }


protected:
    
    void new_individual(Individual::genes_t& dvars) {
        std::shared_ptr<ind_t> ind(new ind_t(dVarBounds_m, dNames_m, constraints_m));
        std::swap(ind->genes_m, dvars);
        individualsToEvaluate_m.push( population_m->add_individual(ind) );
    }

    void new_individual() {
        std::shared_ptr<ind_t> ind(new ind_t(dVarBounds_m, dNames_m, constraints_m));
        individualsToEvaluate_m.push( population_m->add_individual(ind) );
    }

    void new_individual(std::shared_ptr<ind_t> ind) {
        std::shared_ptr<ind_t> return_ind(new ind_t(ind));
        individualsToEvaluate_m.push(
            population_m->add_individual(return_ind) ) ;
    }

private:

    std::shared_ptr< Population<ind_t> > population_m;

    CmdArguments_t args_;

    std::queue<unsigned int> individualsToEvaluate_m;

    std::vector<std::string> dNames_m;
    Optimizer::bounds_t dVarBounds_m;
    Expressions::Named_t constraints_m;

    double mutationProbability_m;
    double recombinationProbability_m;

    double drand(double range) {
        return (range * (double) rand() / (RAND_MAX + 1.0));
    }
};

#endif