GPSymbRegEvalOp.cpp

#include "GPSymbRegEvalOp.h"
#include "ReadData.h"
#include <vector>
#include <boost/math/special_functions/fpclassify.hpp> // isnan, isinf
 
#include "MeanSquareErrorMetric.h"
#include "MeanAbsoluteErrorMetric.h"
#include "MeanAbsolutePercentageErrorMetric.h"
 
#include "IATree.h"
#include "IATerminal.h"
 
 
void GPSymbRegEvalOp::registerParameters(StateP state)
{
    state->getRegistry()->registerEntry("input_file", (voidP)(new std::string), ECF::STRING);
    state->getRegistry()->registerEntry("linear_scaling", (voidP) (new std::string), ECF::STRING);
    state->getRegistry()->registerEntry("error_weights.file", (voidP) (new std::string), ECF::STRING);
    state->getRegistry()->registerEntry("error_metric", (voidP) (new std::string), ECF::STRING);
}
 
 
// called only once, before the evolution – generates training data
bool GPSymbRegEvalOp::initialize(StateP state)
{
    std::stringstream ss;
    voidP sptr = state->getRegistry()->getEntry("linear_scaling");
    ss << *((std::string*) sptr.get());
    std::string ls = ss.str();
    if (ls.compare("true") == 0) {
        linearScaling = true;
    } else {
        linearScaling = false;
    }
    ss.clear();
 
 
    x.clear();
    y.clear();
    f.clear();
    
    // initialize the Evaluator
    eval.initialize();
 
    // check if the parameters are stated (used) in the conf. file
    if (state->getRegistry()->isModified("input_file")) {
        std::string dataFile = *((std::string*) state->getRegistry()->getEntry("input_file").get());
 
        // read from file into Evaluator data
        if (!readDataFromFile(eval.data, dataFile)) {
            return false;
        }
 
        nSamples = eval.data.size();
        nVariables = eval.data[0].size() - 1;
 
        varNames.clear();
 
        for(uint i = 0; i < nVariables; i++) {
            varNames.push_back ("x" + uint2str(i + 1));
        }
 
        GenotypeP genotype = state->getGenotypes()[0];
        Tree::IATreeP iaTree = boost::dynamic_pointer_cast<Tree::IATree>(genotype);
        
        if (iaTree != 0) {
            for(uint i = 0; i < nVariables; i++) {
                IATerminalP terminal = boost::dynamic_pointer_cast<IATerminal>(iaTree->primitiveSet_->getTerminalByName("x" + uint2str(i + 1)));
                assert(terminal != 0);
                double low = eval.data[0][i];
                double up = low;
                for (uint j = 1; j < nSamples; j++) {
                    double d = eval.data[j][i];
                    if (d > up) {
                        up = d;
                    }
                    if (d < low) {
                        low = d;
                    }
                }
 
                terminal->lowerBound = low;
                terminal->upperBound = up;
            }
        }
    }
 
    errorWeightsDefined = false;
 
    if (state->getRegistry()->isModified("error_weights.file")) {
        std::string dataFile = *((std::string*) state->getRegistry()->getEntry("error_weights.file").get());
 
        std::vector< std::vector<double> > error_weights;
 
        // read from file into Evaluator data
        if (!readDataFromFile(error_weights, dataFile)) {
            return false;
        }
 
        errorWeightsDefined = true;
 
        unsigned int size = error_weights.size();
        for (uint i = 0; i < size; i++) {
            errorWeights.push_back(error_weights[i][0]);
        }
 
        if (size < nSamples) {
            std::cout << "Error weights count less than samples count." << std::endl;
            errorWeightsDefined = false;
            return false;
        }
    }
 
    if (state->getRegistry()->isModified("error_metric")) {
        std::string errorMetricStr = *((std::string*) state->getRegistry()->getEntry("error_metric").get());
 
        if (errorMetricStr.compare("mean_square_error") == 0) {
            errorMetric = new MeanSquareErrorMetric;
        } else if (errorMetricStr.compare("mean_absolute_error") == 0) {
            errorMetric = new MeanAbsoluteErrorMetric;
        } else if (errorMetricStr.compare("mean_absolute_percentage_error") == 0) {
            errorMetric = new MeanAbsolutePercentageErrorMetric;
        } else {
            std::cout << "Error metric unknown." << std::endl;
            return false;
        }
    } else {
        // default error metric
        errorMetric = new MeanSquareErrorMetric;
    }
 
    return true;
}
 
 
bool GPSymbRegEvalOp::csvRead(StateP state, std::string entry, std::vector<double>* vec)
{
    std::ifstream stream;
    std::string line;
    voidP sptr = state->getRegistry()->getEntry(entry);
    std::string fname = *((std::string*) sptr.get());
    stream.open(fname.c_str());
    if (!stream.is_open()) return false;
    while (getline(stream, line)){
        vec->push_back(atof(line.c_str()));
    }
    return true;
}
 
 
FitnessP GPSymbRegEvalOp::evaluateUsingLinearScaling(IndividualP individual)
{
    // we try to minimize the function value, so we use FitnessMin fitness (for minimization problems)
    FitnessP fitness (new FitnessMin);
 
    // get the genotype we defined in the configuration file
    Tree::Tree* tree = (Tree::Tree*) individual->getGenotype().get();
    // (you can also use boost smart pointers:)
    //TreeP tree = boost::static_pointer_cast<Tree::Tree> (individual->getGenotype());
 
    //average value of y and t
    double mean_y = 0, mean_t = 0, sum_yt = 0, sum_sqr_y = 0;
 
    for (uint i = 0 ; i < nSamples ; i++)
    {
        // set only defined variables (x1, x2, ...)
        for(uint term = 0; term < nVariables; term++) {
            tree->setTerminalValue(varNames[term], &eval.data[i][term]);
        }
        // get the y value of the current tree
        double y;
        tree->execute(&y);
 
        mean_y += y;
        double t = eval.data[i][nVariables];
        mean_t += t;
        sum_yt += t*y;
        sum_sqr_y += y*y;
    }
 
    mean_y /= nSamples;
    mean_t /= nSamples;
 
    //coefficients for linear regression
    double b = (sum_yt - nSamples*mean_y*mean_t)/(sum_sqr_y - nSamples*mean_y*mean_y);
    double a = mean_t - b*mean_y;
 
    bool isNaN = boost::math::isnan(a) || boost::math::isnan(b);
    bool isInf = boost::math::isinf(a) || boost::math::isinf(b);
    if (isNaN || isInf) {
        fitness->setValue(1e14);
        return fitness;
    }
 
    double value = 0;
    double var_of_y = 0; //variance of variable result
 
    for(uint i = 0; i < nSamples; i++) {
        // set only defined variables (x1, x2, ...)
        for(uint term = 0; term < nVariables; term++) {
            tree->setTerminalValue(varNames[term], &eval.data[i][term]);
            dataWriter->writeData(eval.data[i][term]);
            dataWriter->writeData('\t');
        }
 
        double y;
        tree->execute(&y);
 
        var_of_y += (y-mean_y)*(y-mean_y);
        // add the difference
        double t = eval.data[i][nVariables];
        double result = a + b*y;
        dataWriter->writeData(result);
        dataWriter->writeData('\n');
 
        double error = errorMetric->calculateError(result, t);
 
        if (errorWeightsDefined) {
            error *= errorWeights[i];
        }
 
        value += error*error;
    }
 
    value /= nSamples;
    var_of_y /= nSamples-1;
 
    //delete expresion (but not explicit, just set fitness to big number)
    // if ( fabs(var_of_y) > 1e7 || fabs(var_of_y) < 1e-7 ) {
    //  value = 1e14;
    // }
 
    if (printParams) {
        std::cout << "\nLinear scaling parameters: scale=" << b << " offset=" << a << std::endl;
    }
 
    fitness->setValue(value);
 
    return fitness;
}
 
 
FitnessP GPSymbRegEvalOp::evaluateWithoutLinearScaling(IndividualP individual)
{
    // we try to minimize the function value, so we use FitnessMin fitness (for minimization problems)
    FitnessP fitness(new FitnessMin);
 
    // get the genotype we defined in the configuration file
    Tree::Tree* tree = (Tree::Tree*) individual->getGenotype().get();
 
    // if using the Evaluator
    //readIndividual(individual);
 
    double value = 0;
    double result;
 
    // evaluating data from input file
    for (uint i = 0; i < nSamples; i++) {
        // set only defined variables (x1, x2, ...)
        for(uint term = 0; term < nVariables; term++) {
            tree->setTerminalValue(varNames[term], &eval.data[i][term]);
            dataWriter->writeData(eval.data[i][term]);
            dataWriter->writeData('\t');
        }
 
        // get the f value of the current tree
        tree->execute(&result);
        dataWriter->writeData(result);
        dataWriter->writeData('\n');
        // Evaluator:
        //result = eval.executeParsedExpression(i);
 
        // add the absolute difference
        //value += fabs(eval.data[i][nVariables] - result);
        // or squared error
        double error = errorMetric->calculateError(result, eval.data[i][nVariables]);
        
        if (errorWeightsDefined) {
            error *= errorWeights[i];
        }
 
        value += error;
        // += fabs(f.at(i*nSamples + j) - result);
    }
 
    value /= pow(1.*nSamples, 2);
    fitness->setValue(value);
    return fitness;
}
 
 
FitnessP GPSymbRegEvalOp::evaluate(IndividualP individual)
{
    if (linearScaling) {
        return evaluateUsingLinearScaling(individual);
    } else {
        return evaluateWithoutLinearScaling(individual);
    }
}
 
 
// read an individual into Evaluator parsedExpression - Tree version
void GPSymbRegEvalOp::readIndividual(IndividualP individual)
{
    TreeP tree = boost::dynamic_pointer_cast<Tree::Tree> (individual->getGenotype());
    static std::string strPrimitive;
    uint nTreeSize;
    uint nTrees = (uint) individual->size();
    eval.parsedExpression.resize(nTrees);
 
    // parse multiple trees (if any)
    for(uint iTree = 0; iTree < nTrees; iTree++) {
        TreeP pTree = boost::dynamic_pointer_cast<Tree::Tree> (individual->getGenotype(iTree)); // pokazivac na tree
        nTreeSize = (uint) pTree->size();
        eval.parsedExpression[iTree].resize(nTreeSize);
 
        // read whole tree expression
        for(uint i = 0; i < nTreeSize; i++) {   
            strPrimitive = (*pTree)[i]->primitive_->getName();
 
            // check function names
            for(uint iPrim = 0; iPrim < eval.funcNames.size(); iPrim++)
            {   // TODO: skip unused functions
                //if(!Nodes[iPrim].active)
                //  continue;
                if(strPrimitive == eval.funcNames[iPrim]) {
                    eval.parsedExpression[iTree][i] = iPrim;
                    break;
                }
            }
 
            // check terminals
            for(uint iPrim = 0; iPrim < eval.termNames.size(); iPrim++)
            {   //if(!Nodes[iPrim].active)
                //  continue;
                if(strPrimitive == eval.termNames[iPrim]) {
                    eval.parsedExpression[iTree][i] = iPrim + Evaluator::TERMINALS;
                    break;
                }
            }
        }
    }
}