View all examples

Curve Fitting Problem

Problem

Solving the Curve Fitting Problem consists in constructing a curve, or mathematical function, in order to fit a series of data points as best as possible, and possibly subject to constraints. Curve Fitting can involve either interpolation, where an exact fit to the data is required, or smoothing, in which a “smooth” function is constructed to approximately fit the data. For more details, see Wikipedia.

In the Curve Fitting Problem we consider here, we wish to find the optimal parameters for a defined function to best map a set of inputs to a set of outputs. For example, if we assume the mapping function has the form:

f ( x ) = a sin ⁡ ( b − x ) + c x 2 + d

We want to find the parameters a, b, c, and d for which the mapping function best fits the observations.

Principles learned

Data

Each data file contains:

  • First line: the number of observations
  • Next lines: values of input and output for each observation

Program

The Hexaly model for the Curve Fitting Problem uses float decision variables representing the four parameters a, b, c, and d. Using the ‘sin’ and ‘pow’ nonlinear operators, we can compute the corresponding function predictions for each observation as intermediate expressions.

The objective function is the sum of square errors. Using the intermediate expressions previously computed and the ‘pow’ operator, we can compute the square error for each observation as the squared difference between the output predicted by the mapping function and the observed output.

Download
Execution
hexaly curve_fitting.hxm instances/observations.in [hxTimeLimit=] [solFileName=]
use io;

/* Read instance data */
function input() {

    local usage = "Usage: hexaly curve_fitting.hxm "
        + "inFileName=inputFile [solFileName=outputFile] [hxTimeLimit=timeLimit]";

    if (inFileName == nil) throw usage;

    local inFile = io.openRead(inFileName);
    nbObservations = inFile.readInt();
    for [i in 0...nbObservations] {
        inputs[i] = inFile.readDouble();
        outputs[i] = inFile.readDouble();
    }
}

/* Declare the optimization model */
function model() {

    // Decision variables (parameters of the mapping function)
    a <- float(-100, 100);
    b <- float(-100, 100);
    c <- float(-100, 100);
    d <- float(-100, 100);

    // Minimize square error between prediction and output
    predictions[i in 0...nbObservations] <- a * sin(b - inputs[i]) 
            + c * pow(inputs[i], 2) + d;
    errors[i in 0...nbObservations] <- predictions[i] - outputs[i];
    squareError <- sum[i in 0...nbObservations] (pow(errors[i], 2));
    minimize squareError;
}

/* Parametrize the solver */
function param() {
    if (hxTimeLimit == nil) hxTimeLimit = 3; 
}

/* Write the solution in a file */
function output() {
    if (solFileName == nil) return;
    local solFile = io.openWrite(solFileName);
    solFile.println("Optimal mapping function");
    solFile.println("a = " + a.value);
    solFile.println("b = " + b.value);
    solFile.println("c = " + c.value);
    solFile.println("d = " + d.value);
}
Execution (Windows)
set PYTHONPATH=%HX_HOME%\bin\python
python curve_fitting.py instances\observations.in
Execution (Linux)
export PYTHONPATH=/opt/hexaly_13_0/bin/python
python curve_fitting.py instances/observations.in
import hexaly.optimizer
import sys


def read_float(filename):
    with open(filename) as f:
        return [float(elem) for elem in f.read().split()]

#
# Read instance data
#


def read_instance(instance_file):
    file_it = iter(read_float(instance_file))

    # Number of observations
    nb_observations = int(next(file_it))

    # Inputs and outputs
    inputs = []
    outputs = []
    for i in range(nb_observations):
        inputs.append(next(file_it))
        outputs.append(next(file_it))

    return nb_observations, inputs, outputs


def main(instance_file, output_file, time_limit):
    nb_observations, inputs, outputs = read_instance(instance_file)

    with hexaly.optimizer.HexalyOptimizer() as optimizer:
        #
        # Declare the optimization model
        #
        model = optimizer.model

        # Decision variables (parameters of the mapping function)
        a = model.float(-100, 100)
        b = model.float(-100, 100)
        c = model.float(-100, 100)
        d = model.float(-100, 100)

        # Minimize square error between prediction and output
        predictions = [a * model.sin(b - inputs[i]) + c * inputs[i] ** 2 + d
                       for i in range(nb_observations)]
        errors = [predictions[i] - outputs[i] for i in range(nb_observations)]
        square_error = model.sum(model.pow(errors[i], 2) for i in range(nb_observations))
        model.minimize(square_error)

        model.close()

        # Parameterize the optimizer
        optimizer.param.time_limit = time_limit

        optimizer.solve()

        #
        # Write the solution in a file
        #
        if output_file is not None:
            with open(output_file, 'w') as f:
                f.write("Optimal mapping function\n")
                f.write("a = " + str(a.value) + "\n")
                f.write("b = " + str(b.value) + "\n")
                f.write("c = " + str(c.value) + "\n")
                f.write("d = " + str(d.value) + "\n")


if __name__ == '__main__':
    if len(sys.argv) < 2:
        print("Usage: python curve_fitting.py instance_file [output_file] [time_limit]")
        sys.exit(1)

    instance_file = sys.argv[1]
    output_file = sys.argv[2] if len(sys.argv) >= 3 else None
    time_limit = int(sys.argv[3]) if len(sys.argv) >= 4 else 3
    main(instance_file, output_file, time_limit)
Compilation / Execution (Windows)
cl /EHsc curve_fitting.cpp -I%HX_HOME%\include /link %HX_HOME%\bin\hexaly130.lib
curve_fitting instances\observations.in
Compilation / Execution (Linux)
g++ curve_fitting.cpp -I/opt/hexaly_13_0/include -lhexaly130 -lpthread -o curve_fitting
./curve_fitting instances/observations.in
#include "optimizer/hexalyoptimizer.h"
#include <fstream>
#include <iostream>
#include <vector>

using namespace hexaly;
using namespace std;

class CurveFitting {
public:
    // Number of observations
    int nbObservations;

    // Inputs and outputs
    vector<hxdouble> inputs;
    vector<hxdouble> outputs;

    // Hexaly Optimizer
    HexalyOptimizer optimizer;

    // Decision variables
    HxExpression a, b, c, d;

    // Objective
    HxExpression squareError;

    // Read instance data
    void readInstance(const string& fileName) {
        ifstream infile;
        infile.exceptions(ifstream::failbit | ifstream::badbit);
        infile.open(fileName.c_str());

        infile >> nbObservations;

        inputs.resize(nbObservations);
        outputs.resize(nbObservations);
        for (int i = 0; i < nbObservations; ++i) {
            infile >> inputs[i];
            infile >> outputs[i];
        }
    }

    void solve(int limit) {
        // Declare the optimization model
        HxModel model = optimizer.getModel();

        // Decision variables (parameters of the mapping function)
        a = model.floatVar(-100, 100);
        b = model.floatVar(-100, 100);
        c = model.floatVar(-100, 100);
        d = model.floatVar(-100, 100);

        // Minimize square error between prediction and output
        squareError = model.sum();
        for (int i = 0; i < nbObservations; ++i) {
            HxExpression prediction = a * model.sin(b - inputs[i]) + c * pow(inputs[i], 2) + d;
            HxExpression error = model.pow(prediction - outputs[i], 2);
            squareError.addOperand(error);
        }
        model.minimize(squareError);
        model.close();

        // Parametrize the optimizer
        optimizer.getParam().setTimeLimit(limit);

        optimizer.solve();
    }

    /* Write the solution in a file */
    void writeSolution(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());
        outfile << "Optimal mapping function" << endl;
        outfile << "a = " << a.getDoubleValue() << endl;
        outfile << "b = " << b.getDoubleValue() << endl;
        outfile << "c = " << c.getDoubleValue() << endl;
        outfile << "d = " << d.getDoubleValue() << endl;
    }
};

int main(int argc, char** argv) {
    if (argc < 2) {
        cerr << "Usage: curve_fitting inputFile [outputFile] [timeLimit]" << endl;
        return 1;
    }

    const char* instanceFile = argv[1];
    const char* solFile = argc > 2 ? argv[2] : NULL;
    const char* strTimeLimit = argc > 3 ? argv[3] : "3";

    try {
        CurveFitting model;
        model.readInstance(instanceFile);
        model.solve(atoi(strTimeLimit));
        if (solFile != NULL)
            model.writeSolution(solFile);
        return 0;
    } catch (const exception& e) {
        cerr << "An error occurred: " << e.what() << endl;
        return 1;
    }
}
Compilation / Execution (Windows)
copy %HX_HOME%\bin\Hexaly.NET.dll .
csc CurveFitting.cs /reference:Hexaly.NET.dll
CurveFitting instances\observations.in
using System;
using System.IO;
using System.Collections.Generic;
using System.Globalization;
using Hexaly.Optimizer;

public class CurveFitting : IDisposable
{
    // Number of observations
    int nbObservations;

    // Inputs and outputs
    double[] inputs;
    double[] outputs;

    // Hexaly Optimizer
    HexalyOptimizer optimizer;

    // Decision variables
    HxExpression a, b, c, d;

    // Objective
    HxExpression squareError;

    public CurveFitting()
    {
        optimizer = new HexalyOptimizer();
    }

    // Read instance data
    void ReadInstance(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            nbObservations = int.Parse(input.ReadLine());
            inputs = new double[nbObservations];
            outputs = new double[nbObservations];

            for (int i = 0; i < nbObservations; ++i)
            {
                string[] splittedObs = input.ReadLine().Split(' ');
                inputs[i] = double.Parse(splittedObs[0], CultureInfo.InvariantCulture);
                outputs[i] = double.Parse(splittedObs[1], CultureInfo.InvariantCulture);
            }
        }
    }

    public void Dispose()
    {
        if (optimizer != null)
            optimizer.Dispose();
    }

    void Solve(int limit)
    {
        // Declare the optimization model
        HxModel model = optimizer.GetModel();

        // Decision variables (parameters of the mapping function)
        a = model.Float(-100, 100);
        b = model.Float(-100, 100);
        c = model.Float(-100, 100);
        d = model.Float(-100, 100);

        // Minimize square error between prediction and output
        squareError = model.Sum();
        for (int i = 0; i < nbObservations; ++i)
        {
            HxExpression prediction = a * model.Sin(b - inputs[i]) + c * Math.Pow(inputs[i], 2) + d;
            HxExpression error = model.Pow(prediction - outputs[i], 2);
            squareError.AddOperand(error);
        }
        model.Minimize(squareError);
        model.Close();

        // Parameterize the optimizer
        optimizer.GetParam().SetTimeLimit(limit);

        optimizer.Solve();
    }

    /* Write the solution in a file */
    void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine("Optimal mapping function");
            output.WriteLine("a = " + a.GetDoubleValue());
            output.WriteLine("b = " + b.GetDoubleValue());
            output.WriteLine("c = " + c.GetDoubleValue());
            output.WriteLine("d = " + d.GetDoubleValue());
        }
    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: CurveFitting inputFile [solFile] [timeLimit]");
            Environment.Exit(1);
        }
        string instanceFile = args[0];
        string outputFile = args.Length > 1 ? args[1] : null;
        string strTimeLimit = args.Length > 2 ? args[2] : "3";

        using (CurveFitting model = new CurveFitting())
        {
            model.ReadInstance(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }
}
Compilation / Execution (Windows)
javac CurveFitting.java -cp %HX_HOME%\bin\hexaly.jar
java -cp %HX_HOME%\bin\hexaly.jar;. CurveFitting instances\observations.in
 
Compilation / Execution (Linux)
javac CurveFitting.java -cp /opt/hexaly_13_0/bin/hexaly.jar
java -cp /opt/hexaly_13_0/bin/hexaly.jar:. CurveFitting instances/observations.in
import java.util.*;
import java.io.*;
import com.hexaly.optimizer.*;

public class CurveFitting {
    // Number of observations
    private int nbObservations;

    // Inputs and outputs
    private double[] inputs;
    private double[] outputs;

    // Hexaly Optimizer
    private final HexalyOptimizer optimizer;

    // Decision variables
    private HxExpression a, b, c, d;

    // Objective
    private HxExpression squareError;

    private CurveFitting(HexalyOptimizer optimizer) {
        this.optimizer = optimizer;
    }

    // Read instance data
    private void readInstance(String fileName) throws IOException {
        try (Scanner input = new Scanner(new File(fileName))) {
            input.useLocale(Locale.ROOT);

            nbObservations = input.nextInt();

            inputs = new double[nbObservations];
            outputs = new double[nbObservations];
            for (int i = 0; i < nbObservations; ++i) {
                inputs[i] = input.nextDouble();
                outputs[i] = input.nextDouble();
            }
        }
    }

    private void solve(int limit) {
        // Declare the optimization model
        HxModel model = optimizer.getModel();

        // Decision variables (parameters of the mapping function)
        a = model.floatVar(-100, 100);
        b = model.floatVar(-100, 100);
        c = model.floatVar(-100, 100);
        d = model.floatVar(-100, 100);

        // Minimize square error between prediction and output
        squareError = model.sum();
        for (int i = 0; i < nbObservations; ++i) {
            HxExpression prediction = model.sum(model.prod(a, model.sin(model.sub(b, inputs[i]))),
                    model.prod(c, Math.pow(inputs[i], 2)), d);
            HxExpression error = model.pow(model.sub(prediction, outputs[i]), 2);
            squareError.addOperand(error);
        }

        model.minimize(squareError);
        model.close();

        // Parameterize the optimizer
        optimizer.getParam().setTimeLimit(limit);

        optimizer.solve();
    }

    /* Write the solution in a file */
    private void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(fileName)) {
            output.println("Optimal mapping function");
            output.println("a = " + a.getDoubleValue());
            output.println("b = " + b.getDoubleValue());
            output.println("c = " + c.getDoubleValue());
            output.println("d = " + d.getDoubleValue());
        }
    }

    public static void main(String[] args) {
        if (args.length < 1) {
            System.err.println("Usage: java CurveFitting inputFile [outputFile] [timeLimit]");
            System.exit(1);
        }

        String instanceFile = args[0];
        String outputFile = args.length > 1 ? args[1] : null;
        String strTimeLimit = args.length > 2 ? args[2] : "3";

        try (HexalyOptimizer optimizer = new HexalyOptimizer()) {
            CurveFitting model = new CurveFitting(optimizer);
            model.readInstance(instanceFile);
            model.solve(Integer.parseInt(strTimeLimit));
            if (outputFile != null) {
                model.writeSolution(outputFile);
            }
        } catch (Exception ex) {
            System.err.println(ex);
            ex.printStackTrace();
            System.exit(1);
        }
    }
}

Similar examples

Discover the ease of use and performance of Hexaly