help with implementation of Linear programming solver.
CodeGuru Home VC++ / MFC / C++ .NET / C# Visual Basic VB Forums Developer.com
Page 1 of 3 123 LastLast
Results 1 to 15 of 44

Thread: help with implementation of Linear programming solver.

  1. #1
    Join Date
    Sep 2013
    Posts
    36

    help with implementation of Linear programming solver.

    I am trying to implement the linear programming solver.

    This is the header file of the linear programming solver :

    /*!
    \internal

    Representation of a LP constraint like:

    (c1 * X1) + (c2 * X2) + ... = K
    or <= K
    or >= K

    Where (ci, Xi) are the pairs in "variables" and K the real "constant".
    */


    struct QSimplexConstraint
    {
    QSimplexConstraint() : constant(0), ratio(Equal), artificial(0) {}`

    enum Ratio {
    LessOrEqual = 0,
    Equal,
    MoreOrEqual
    };

    QHash<QSimplexVariable *, qreal> variables;
    qreal constant;
    Ratio ratio;

    QPair<QSimplexVariable *, qreal> helper;
    QSimplexVariable * artificial;

    void invert();

    bool isSatisfied() {
    qreal leftHandSide(0);

    QHash<QSimplexVariable *, qreal>::const_iterator iter;
    for (iter = variables.constBegin(); iter != variables.constEnd(); ++iter) {
    leftHandSide += iter.value() * iter.key()->result;
    }

    Q_ASSERT(constant > 0 || qFuzzyCompare(1, 1 + constant));

    if ((leftHandSide == constant) || qAbs(leftHandSide - constant) < 0.0000001)
    return true;

    switch (ratio) {
    case LessOrEqual:
    return leftHandSide < constant;
    case MoreOrEqual:
    return leftHandSide > constant;
    default:
    return false;
    }
    }

    class QSimplex
    {
    public:
    QSimplex();
    virtual ~QSimplex();
    qreal solveMin();
    qreal solveMax();



    bool setConstraints(const QList<QSimplexConstraint *> constraints);
    void setObjective(QSimplexConstraint *objective);

    void dumpMatrix();

    private:
    // Matrix handling
    qreal valueAt(int row, int column);
    void setValueAt(int row, int column, qreal value);
    void clearRow(int rowIndex);
    void clearColumns(int first, int last);
    void combineRows(int toIndex, int fromIndex, qreal factor);

    // Simplex
    bool simplifyConstraints(QList<QSimplexConstraint *> *constraints);
    int findPivotColumn();
    int pivotRowForColumn(int column);
    void reducedRowEchelon();
    bool iterate();

    // Helpers
    void clearDataStructures();
    void solveMaxHelper();
    enum solverFactor { Minimum = -1, Maximum = 1 };
    qreal solver(solverFactor factor);
    void collectResults();

    QList<QSimplexConstraint *> constraints;
    QList<QSimplexVariable *> variables;
    QSimplexConstraint *objective;

    int rows;
    int columns;
    int firstArtificial;

    qreal *matrix;
    };

    inline qreal QSimplex::valueAt(int rowIndex, int columnIndex)
    {
    return matrix[rowIndex * columns + columnIndex];
    }

    inline void QSimplex::setValueAt(int rowIndex, int columnIndex, qreal value)
    {
    matrix[rowIndex * columns + columnIndex] = value;
    }




    I have a text file with some data like:

    Maximize:

    obj: 3e-06 A - 3e-06 B + 2.7e-01 F

    constraints:

    RXN1: -1 A + 1 B -1 C + 1 D -1 E -1 F = 0

    RXN2: -1 A + 1 B -1 C + 1 D -1 E -1 F = 0

    RXN3: -1 A + 1 B -1 C + 1 D -1 E -1 F = 0

    RXN4: -1 A + 1 B -1 C + 1 D -1 E -1 F = 0

    ... many constraints like this

    Bounds:

    A >= 0
    B <= 100
    C >= 0

    .....
    ...........many bounds like this.

    I wrote a small program to get the data from this file:

    #include <iostream>
    #include <fstream>
    #include <string>


    using namespace std;

    int main()
    {
    ifstream ifs("Metabolism.txt");
    string line;
    while(getline(ifs,line)) {
    cout << "[ " << line << " ]" << endl;
    }
    system("PAUSE");
    }

    I want some help to parse all the constraints and bounds as inputs and get the maximum value of the objective function as the output using the above lpsolver header file.I have also attached the sample file below.

    Thanks.
    Attached Files Attached Files

  2. #2
    Join Date
    Apr 1999
    Posts
    27,426

    Re: help with implementation of Linear programming solver.

    Quote Originally Posted by navy1991 View Post
    I want some help to parse all the constraints and bounds as inputs
    First, use code tags when posting code. The code you posted is practically unreadable without code tags.

    What is the exact problem do you have in parsing the input? Where are the data structures in your code to store the parsed input?
    Code:
    string line;
    while(getline(ifs,line)) {
    Parsing is much more than just reading a string -- you have to store that data in some structure and in some format, let alone devise a way to interpret the string in the way that satisfies your conditions. Did you have this designed in any way (not coded, but designed)? If not, you need to start with the design first.
    and get the maximum value of the objective function as the output using the above lpsolver header file.I have also attached the sample file below.
    First things first. Design and write the code to parse the input and place in proper structures and format so that the rest of the code you have knows what to do with the data.

    Regards,

    Paul McKenzie

  3. #3
    Join Date
    Sep 2013
    Posts
    36

    Re: help with implementation of Linear programming solver.

    Hi Paul,
    Point noted. I will use code tags.I am trying to design the data structure so that I can manipulate it. Please give me some time to update the design.

    Thanks,
    navy1991.

  4. #4
    Join Date
    Sep 2013
    Posts
    36

    Re: help with implementation of Linear programming solver.

    Designing the data structure:

    As the representation of a linear programming problem, I need to collect the pairs in Variables (ci, Xi) (RHS) and then K as real constant (LHS) and then the ratio between. So 3 things : RHS, LHS , Ratio between.
    Code:
    /*!
    \internal
    
    Representation of a LP constraint like:
    
    (c1 * X1) + (c2 * X2) + ... = K
    or <= K
    or >= K
    
    Where (ci, Xi) are the pairs in "variables" and K the real "constant".
    */
    Referring to the code below:

    Code:
    struct QSimplexConstraint
    {
    QSimplexConstraint() : constant(0), ratio(Equal), artificial(0) {}`
    
    enum Ratio {
    LessOrEqual = 0,
    Equal,
    MoreOrEqual
    };
    I understand that they have designed a struct named Qsimplexconstraint . Inside this data structure, they have QsimplexConstraint class (bit confusing -- I guess constant refer to 'K' (RHS) ; ratio refers to the connection between RHS and LHS of the equation. (= , >= , <=) ; artificial refers to the RHS (pairs in variables) ) and enum Ratio (enumerator to choose = , >= , <=).

    QList<QSimplexConstraint *> constraints;

    This helps me to understand that they have created QList Constraints (I guess it is a special type in c++ Qt ). Could you please explain more about QList ?



    From
    Code:
    QHash<QSimplexVariable *, qreal> variables;
    qreal constant;
    Ratio ratio;
    I understand that they have created a special QHash variables (Qhash means a way to search by splitting the data into equal halves and to continue)

    From
    Code:
    QPair<QSimplexVariable *, qreal> helper;
    QSimplexVariable * artificial;
    I understand they have used Qpair to pair the QsimplexVariable pointer named artificial (RHS -- (c1, x1)) to qreal (LHS).

    In summary, I am clear about three parts :

    1.) RHS (right handside)
    2.) LHS
    3.) Ratio between.

    This is what I have understood till now. Could you please tell me what should I do next.

    Thanks.

  5. #5
    Join Date
    Apr 1999
    Posts
    27,426

    Re: help with implementation of Linear programming solver.

    Did you write that template class yourself? If not, then shouldn't there have been examples of how to use it? What good is a class if no one knows how it's used?
    This is what I have understood till now. Could you please tell me what should I do next.
    Contact the author(s) of the class. Again, no one writes something like this without showing how to use it.

    Regards,

    Paul McKenzie

  6. #6
    Join Date
    Sep 2013
    Posts
    36

    Re: help with implementation of Linear programming solver.

    Hi paul,
    This is an assignment problem in one of the research lab. I did not write the LP solver.It comes under Operations Research. I just want to use it to solve my problem. My problem is called Flux Balance Analysis. I managed to get the source file of the solver too (.cpp file). I have pasted it below :

    Code:
    QT_BEGIN_NAMESPACE
     
    QSimplex::QSimplex() : objective(0), rows(0), columns(0), firstArtificial(0), matrix(0)
    {
    }
     
    QSimplex::~QSimplex()
    {
        clearDataStructures();
    }
     
    void QSimplex::clearDataStructures()
    {
        if (matrix == 0)
            return;
     
        // Matrix
        rows = 0;
        columns = 0;
        firstArtificial = 0;
        free(matrix);
        matrix = 0;
     
        // Constraints
        for (int i = 0; i < constraints.size(); ++i) {
            delete constraints[i]->helper.first;
            constraints[i]->helper.first = 0;
            constraints[i]->helper.second = 0.0;
            delete constraints[i]->artificial;
            constraints[i]->artificial = 0;
        }
        constraints.clear();
     
        // Other
        variables.clear();
        objective = 0;
    }
     
    void QSimplex::setConstraints(const QList<QSimplexConstraint *> newConstraints)
    {
        clearDataStructures();
     
        if (newConstraints.isEmpty())
            return;
        constraints = newConstraints;
     
        // Set Variables direct mapping
        QSet<QSimplexVariable *> variablesSet;
        for (int i = 0; i < constraints.size(); ++i)
            variablesSet += \
                QSet<QSimplexVariable *>::fromList(constraints[i]->variables.keys());
        variables = variablesSet.toList();
     
        // Set Variables reverse mapping
        for (int i = 0; i < variables.size(); ++i) {
            // The variable "0" goes at the column "1", etc...
            variables[i]->index = i + 1;
        }
     
        // Normalize Constraints
        int variableIndex = variables.size();
        QList <QSimplexVariable *> artificialList;
     
        for (int i = 0; i < constraints.size(); ++i) {
            QSimplexVariable *slack;
            QSimplexVariable *surplus;
            QSimplexVariable *artificial;
     
            Q_ASSERT(constraints[i]->helper.first == 0);
            Q_ASSERT(constraints[i]->artificial == 0);
     
            switch(constraints[i]->ratio) {
            case QSimplexConstraint::LessOrEqual:
                slack = new QSimplexVariable;
                slack->index = ++variableIndex;
                constraints[i]->helper.first = slack;
                constraints[i]->helper.second = 1.0;
                break;
            case QSimplexConstraint::MoreOrEqual:
                surplus = new QSimplexVariable;
                surplus->index = ++variableIndex;
                constraints[i]->helper.first = surplus;
                constraints[i]->helper.second = -1.0;
                // fall through
            case QSimplexConstraint::Equal:
                artificial = new QSimplexVariable;
                constraints[i]->artificial = artificial;
                artificialList += constraints[i]->artificial;
                break;
            }
        }
     
        firstArtificial = variableIndex + 1;
        for (int i = 0; i < artificialList.size(); ++i)
            artificialList[i]->index = ++variableIndex;
        artificialList.clear();
     
        // Matrix
     
        // One for each variable plus the Basic and BFS columns (first and last)
        columns = variableIndex + 2;
        // One for each constraint plus the objective function
        rows = constraints.size() + 1;
     
        matrix = (qreal *)malloc(sizeof(qreal) * columns * rows);
        if (!matrix) {
            qWarning() << "QSimplex: Unable to allocate memory!";
            return;
        }
        for (int i = columns * rows - 1; i >= 0; --i)
            matrix[i] = 0.0;
     
        // Fill Matrix
        for (int i = 1; i <= constraints.size(); ++i) {
            QSimplexConstraint *c = constraints[i - 1];
     
            if (c->artificial) {
                // Will use artificial basic variable
                setValueAt(i, 0, c->artificial->index);
                setValueAt(i, c->artificial->index, 1.0);
     
                if (c->helper.second != 0.0) {
                    // Surplus variable
                    setValueAt(i, c->helper.first->index, c->helper.second);
                }
            } else {
                // Slack is used as the basic variable
                Q_ASSERT(c->helper.second == 1.0);
                setValueAt(i, 0, c->helper.first->index);
                setValueAt(i, c->helper.first->index, 1.0);
            }
     
            QHash<QSimplexVariable *, qreal>::const_iterator iter;
            for (iter = c->variables.constBegin();
                 iter != c->variables.constEnd();
                 ++iter) {
                setValueAt(i, iter.key()->index, iter.value());
            }
     
            setValueAt(i, columns - 1, c->constant);
        }
     
        // Set temporary objective: -1 * sum_of_artificial_vars
        for (int j = firstArtificial; j < columns - 1; ++j)
            setValueAt(0, j, 1.0);
     
        // Maximize our objective (artificial vars go to zero)
        solveMaxHelper();
     
        if (valueAt(0, columns - 1) != 0.0) {
            qWarning() << "QSimplex: No feasible solution!";
            clearDataStructures();
            return;
        }
     
        // Remove artificial variables
        clearColumns(firstArtificial, columns - 2);
    }
     
    void QSimplex::solveMaxHelper()
    {
        reducedRowEchelon();
        while (iterate()) ;
    }
     
    void QSimplex::setObjective(QSimplexConstraint *newObjective)
    {
        objective = newObjective;
    }
     
    void QSimplex::clearRow(int rowIndex)
    {
        qreal *item = matrix + rowIndex * columns;
        for (int i = 0; i < columns; ++i)
            item[i] = 0.0;
    }
     
    void QSimplex::clearColumns(int first, int last)
    {
        for (int i = 0; i < rows; ++i) {
            qreal *row = matrix + i * columns;
            for (int j = first; j <= last; ++j)
                row[j] = 0.0;
        }
    }
     
    void QSimplex::dumpMatrix()
    {
        printf("---- Simplex Matrix ----\n");
     
        printf("       ");
        for (int j = 0; j < columns; ++j)
            printf("  <% 2d >", j);
        printf("\n");
     
        for (int i = 0; i < rows; ++i) {
            printf("Row %2d:", i);
     
            qreal *row = matrix + i * columns;
            for (int j = 0; j < columns; ++j) {
                printf("  % 2.2f", row[j]);
            }
            printf("\n");
        }
        printf("------------------------\n\n");
    }
     
    void QSimplex::combineRows(int toIndex, int fromIndex, qreal factor)
    {
        if (!factor)
            return;
     
        qreal *from = matrix + fromIndex * columns;
        qreal *to = matrix + toIndex * columns;
     
        for (int j = 1; j < columns; ++j) {
            qreal value = from[j];
     
            // skip to[j] = to[j] + factor*0.0
            if (value == 0.0)
                continue;
     
            to[j] += factor * value;
     
            // ### Avoid Numerical errors
            if (qAbs(to[j]) < 0.0000000001)
                to[j] = 0.0;
        }
    }
     
    int QSimplex::findPivotColumn()
    {
        qreal min = 0;
        int minIndex = -1;
     
        for (int j = 0; j < columns-1; ++j) {
            if (valueAt(0, j) < min) {
                min = valueAt(0, j);
                minIndex = j;
            }
        }
     
        return minIndex;
    }
     
    int QSimplex::pivotRowForColumn(int column)
    {
        qreal min = 999999999999.0; // ###
        int minIndex = -1;
     
        for (int i = 1; i < rows; ++i) {
            qreal divisor = valueAt(i, column);
            if (divisor <= 0)
                continue;
     
            qreal quotient = valueAt(i, columns - 1) / divisor;
            if (quotient < min) {
                min = quotient;
                minIndex = i;
            }
        }
     
        return minIndex;
    }
     
    void QSimplex::reducedRowEchelon()
    {
        for (int i = 1; i < rows; ++i) {
            int factorInObjectiveRow = valueAt(i, 0);
            combineRows(0, i, -1 * valueAt(0, factorInObjectiveRow));
        }
    }
     
    bool QSimplex::iterate()
    {
        // Find Pivot column
        int pivotColumn = findPivotColumn();
        if (pivotColumn == -1)
            return false;
     
        // Find Pivot row for column
        int pivotRow = pivotRowForColumn(pivotColumn);
        if (pivotRow == -1) {
            qWarning() << "QSimplex: Unbounded problem!";
            return false;
        }
     
        // Normalize Pivot Row
        qreal pivot = valueAt(pivotRow, pivotColumn);
        if (pivot != 1.0)
            combineRows(pivotRow, pivotRow, (1.0 - pivot) / pivot);
     
        // Update other rows
        for (int row=0; row < rows; ++row) {
            if (row == pivotRow)
                continue;
     
            combineRows(row, pivotRow, -1 * valueAt(row, pivotColumn));
        }
     
        // Update first column
        setValueAt(pivotRow, 0, pivotColumn);
     
        //    dumpMatrix();
        //    printf("------------ end of iteration --------------\n");
        return true;
    }
     
    /*!
      \internal
     
      Both solveMin and solveMax are interfaces to this method.
     
      The enum solverFactor admits 2 values: Minimum (-1) and Maximum (+1).
     */
    qreal QSimplex::solver(solverFactor factor)
    {
        // Remove old objective
        clearRow(0);
     
        // Set new objective
        QHash<QSimplexVariable *, qreal>::const_iterator iter;
        for (iter = objective->variables.constBegin();
             iter != objective->variables.constEnd();
             ++iter) {
            setValueAt(0, iter.key()->index, -1 * factor * iter.value());
        }
     
        solveMaxHelper();
        collectResults();
     
        return factor * valueAt(0, columns - 1);
    }
     
    qreal QSimplex::solveMin()
    {
        return solver(Minimum);
    }
     
    qreal QSimplex::solveMax()
    {
        return solver(Maximum);
    }
     
    void QSimplex::collectResults()
    {
        // All variables are zero unless overridden below.
     
        // ### Is this really needed? Is there any chance that an
        // important variable remains as non-basic at the end of simplex?
        for (int i = 0; i < variables.size(); ++i)
            variables[i]->result = 0;
     
        // Basic variables
        // Update the variable indicated in the first column with the value
        // in the last column.
        for (int i = 1; i < rows; ++i) {
            int index = valueAt(i, 0) - 1;
            if (index < variables.size())
                variables[index]->result = valueAt(i, columns - 1);
        }
    }
     
    QT_END_NAMESPACE
    and here goes the header file :

    Code:
    QT_BEGIN_NAMESPACE
    
    struct QSimplexVariable
    {
        QSimplexVariable() : result(0), index(0) {}
    
        qreal result;
        int index;
    };
    
    
    /*!
      \internal
    
      Representation of a LP constraint like:
    
        (c1 * X1) + (c2 * X2) + ...  =  K
                                 or <=  K
                                 or >=  K
    
        Where (ci, Xi) are the pairs in "variables" and K the real "constant".
    */
    struct QSimplexConstraint
    {
        QSimplexConstraint() : constant(0), ratio(Equal), artificial(0) {}
    
        enum Ratio {
            LessOrEqual = 0,
            Equal,
            MoreOrEqual
        };
    
        QHash<QSimplexVariable *, qreal> variables;
        qreal constant;
        Ratio ratio;
    
        QPair<QSimplexVariable *, qreal> helper;
        QSimplexVariable * artificial;
    
        void invert();
    
        bool isSatisfied() {
            qreal leftHandSide(0);
    
            QHash<QSimplexVariable *, qreal>::const_iterator iter;
            for (iter = variables.constBegin(); iter != variables.constEnd(); ++iter) {
                leftHandSide += iter.value() * iter.key()->result;
            }
    
            Q_ASSERT(constant > 0 || qFuzzyCompare(1, 1 + constant));
    
            if ((leftHandSide == constant) || qAbs(leftHandSide - constant) < 0.0000001)
                return true;
    
            switch (ratio) {
            case LessOrEqual:
                return leftHandSide < constant;
            case MoreOrEqual:
                return leftHandSide > constant;
            default:
                return false;
            }
        }
    
    #ifdef QT_DEBUG
        QString toString() {
            QString result;
            result += QString::fromAscii("-- QSimplexConstraint %1 --").arg(quintptr(this), 0, 16);
    
            QHash<QSimplexVariable *, qreal>::const_iterator iter;
            for (iter = variables.constBegin(); iter != variables.constEnd(); ++iter) {
                result += QString::fromAscii("  %1 x %2").arg(iter.value()).arg(quintptr(iter.key()), 0, 16);
            }
    
            switch (ratio) {
            case LessOrEqual:
                result += QString::fromAscii("  (less) <= %1").arg(constant);
                break;
            case MoreOrEqual:
                result += QString::fromAscii("  (more) >= %1").arg(constant);
                break;
            default:
                result += QString::fromAscii("  (eqal) == %1").arg(constant);
            }
    
            return result;
        }
    #endif
    };
    
    class QSimplex
    {
    public:
        QSimplex();
        virtual ~QSimplex();
    
        qreal solveMin();
        qreal solveMax();
    
        bool setConstraints(const QList<QSimplexConstraint *> constraints);
        void setObjective(QSimplexConstraint *objective);
    
        void dumpMatrix();
    
    private:
        // Matrix handling
        qreal valueAt(int row, int column);
        void setValueAt(int row, int column, qreal value);
        void clearRow(int rowIndex);
        void clearColumns(int first, int last);
        void combineRows(int toIndex, int fromIndex, qreal factor);
    
        // Simplex
        bool simplifyConstraints(QList<QSimplexConstraint *> *constraints);
        int findPivotColumn();
        int pivotRowForColumn(int column);
        void reducedRowEchelon();
        bool iterate();
    
        // Helpers
        void clearDataStructures();
        void solveMaxHelper();
        enum solverFactor { Minimum = -1, Maximum = 1 };
        qreal solver(solverFactor factor);
        void collectResults();
    
        QList<QSimplexConstraint *> constraints;
        QList<QSimplexVariable *> variables;
        QSimplexConstraint *objective;
    
        int rows;
        int columns;
        int firstArtificial;
    
        qreal *matrix;
    };
    
    inline qreal QSimplex::valueAt(int rowIndex, int columnIndex)
    {
        return matrix[rowIndex * columns + columnIndex];
    }
    
    inline void QSimplex::setValueAt(int rowIndex, int columnIndex, qreal value)
    {
        matrix[rowIndex * columns + columnIndex] = value;
    }
    
    QT_END_NAMESPACE

  7. #7
    Join Date
    Apr 1999
    Posts
    27,426

    Re: help with implementation of Linear programming solver.

    Quote Originally Posted by navy1991 View Post
    Hi paul,
    This is an assignment problem in one of the research lab. I did not write the LP solver.
    I didn't write the C++ standard library, but I have plenty of examples of how to use it. So it doesn't matter if you wrote it or not. It also doesn't matter what this class is supposed to do, as your problem is much more higher level than this.

    Where is the simple main() program showing basic usage of this class? How are you supposed to know what to do if you don't have a simple example of how to use the class? A class just by itself with no guidance on how to use the public interface of the class is useless.

    Regards,

    Paul McKenzie

  8. #8
    Join Date
    Sep 2013
    Posts
    36

    Re: help with implementation of Linear programming solver.

    They have attached more comments in this source file :

    Code:
    /*!
      \internal
      \class QSimplex
    
      The QSimplex class is a Linear Programming problem solver based on the two-phase
      simplex method.
    
      It takes a set of QSimplexConstraints as its restrictive constraints and an
      additional QSimplexConstraint as its objective function. Then methods to maximize
      and minimize the problem solution are provided.
    
      The two-phase simplex method is based on the following steps:
      First phase:
      1.a) Modify the original, complex, and possibly not feasible problem, into a new,
           easy to solve problem.
      1.b) Set as the objective of the new problem, a feasible solution for the original
           complex problem.
      1.c) Run simplex to optimize the modified problem and check whether a solution for
           the original problem exists.
    
      Second phase:
      2.a) Go back to the original problem with the feasibl (but not optimal) solution
           found in the first phase.
      2.b) Set the original objective.
      3.c) Run simplex to optimize the original problem towards its optimal solution.
    */
    
    /*!
      \internal
    */
    QSimplex::QSimplex() : objective(0), rows(0), columns(0), firstArtificial(0), matrix(0)
    {
    }
    
    /*!
      \internal
    */
    QSimplex::~QSimplex()
    {
        clearDataStructures();
    }
    
    /*!
      \internal
    */
    void QSimplex::clearDataStructures()
    {
        if (matrix == 0)
            return;
    
        // Matrix
        rows = 0;
        columns = 0;
        firstArtificial = 0;
        free(matrix);
        matrix = 0;
    
        // Constraints
        for (int i = 0; i < constraints.size(); ++i) {
            delete constraints[i]->helper.first;
            delete constraints[i]->artificial;
            delete constraints[i];
        }
        constraints.clear();
    
        // Other
        variables.clear();
        objective = 0;
    }
    
    /*!
      \internal
      Sets the new constraints in the simplex solver and returns whether the problem
      is feasible.
    
      This method sets the new constraints, normalizes them, creates the simplex matrix
      and runs the first simplex phase.
    */
    bool QSimplex::setConstraints(const QList<QSimplexConstraint *> newConstraints)
    {
        ////////////////////////////
        // Reset to initial state //
        ////////////////////////////
        clearDataStructures();
    
        if (newConstraints.isEmpty())
            return true;    // we are ok with no constraints
    
       // Make deep copy of constraints. We need this copy because we may change
       // them in the simplification method.
        for (int i = 0; i < newConstraints.size(); ++i) {
           QSimplexConstraint *c = new QSimplexConstraint;
           c->constant = newConstraints[i]->constant;
           c->ratio = newConstraints[i]->ratio;
           c->variables = newConstraints[i]->variables;
           constraints << c;
       }
    
        // Remove constraints of type Var == K and replace them for their value.
        if (!simplifyConstraints(&constraints)) {
            qWarning() << "QSimplex: No feasible solution!";
            clearDataStructures();
           return false;
       }
        ///////////////////////////////////////
        // Prepare variables and constraints //
        ///////////////////////////////////////
    
        // Set Variables direct mapping.
        // "variables" is a list that provides a stable, indexed list of all variables
        // used in this problem.
        QSet<QSimplexVariable *> variablesSet;
        for (int i = 0; i < constraints.size(); ++i)
            variablesSet += \
                QSet<QSimplexVariable *>::fromList(constraints[i]->variables.keys());
        variables = variablesSet.toList();
    
        // Set Variables reverse mapping
        // We also need to be able to find the index for a given variable, to do that
        // we store in each variable its index.
        for (int i = 0; i < variables.size(); ++i) {
            // The variable "0" goes at the column "1", etc...
            variables[i]->index = i + 1;
        }
    
        // Normalize Constraints
        // In this step, we prepare the constraints in two ways:
        // Firstly, we modify all constraints of type "LessOrEqual" or "MoreOrEqual"
        // by the adding slack or surplus variables and making them "Equal" constraints.
        // Secondly, we need every single constraint to have a direct, easy feasible
        // solution. Constraints that have slack variables are already easy to solve,
        // to all the others we add artificial variables.
        //
        // At the end we modify the constraints as follows:
        //  - LessOrEqual: SLACK variable is added.
        //  - Equal: ARTIFICIAL variable is added.
        //  - More or Equal: ARTIFICIAL and SURPLUS variables are added.
        int variableIndex = variables.size();
        QList <QSimplexVariable *> artificialList;
    
        for (int i = 0; i < constraints.size(); ++i) {
            QSimplexVariable *slack;
            QSimplexVariable *surplus;
            QSimplexVariable *artificial;
    
            Q_ASSERT(constraints[i]->helper.first == 0);
            Q_ASSERT(constraints[i]->artificial == 0);
    
            switch(constraints[i]->ratio) {
            case QSimplexConstraint::LessOrEqual:
                slack = new QSimplexVariable;
                slack->index = ++variableIndex;
                constraints[i]->helper.first = slack;
                constraints[i]->helper.second = 1.0;
                break;
            case QSimplexConstraint::MoreOrEqual:
                surplus = new QSimplexVariable;
                surplus->index = ++variableIndex;
                constraints[i]->helper.first = surplus;
                constraints[i]->helper.second = -1.0;
                // fall through
            case QSimplexConstraint::Equal:
                artificial = new QSimplexVariable;
                constraints[i]->artificial = artificial;
                artificialList += constraints[i]->artificial;
                break;
            }
        }
    
        // All original, slack and surplus have already had its index set
        // at this point. We now set the index of the artificial variables
        // as to ensure they are at the end of the variable list and therefore
        // can be easily removed at the end of this method.
        firstArtificial = variableIndex + 1;
        for (int i = 0; i < artificialList.size(); ++i)
            artificialList[i]->index = ++variableIndex;
        artificialList.clear();
    
        /////////////////////////////
        // Fill the Simplex matrix //
        /////////////////////////////
    
        // One for each variable plus the Basic and BFS columns (first and last)
        columns = variableIndex + 2;
        // One for each constraint plus the objective function
        rows = constraints.size() + 1;
        matrix = (qreal *)malloc(sizeof(qreal) * columns * rows);
        if (!matrix) {
            qWarning() << "QSimplex: Unable to allocate memory!";
            return false;
       }
        for (int i = columns * rows - 1; i >= 0; --i)
            matrix[i] = 0.0;
    
        // Fill Matrix
        for (int i = 1; i <= constraints.size(); ++i) {
            QSimplexConstraint *c = constraints[i - 1];
    
            if (c->artificial) {
                // Will use artificial basic variable
                setValueAt(i, 0, c->artificial->index);
                setValueAt(i, c->artificial->index, 1.0);
    
                if (c->helper.second != 0.0) {
                    // Surplus variable
                    setValueAt(i, c->helper.first->index, c->helper.second);
                }
           } else {
                // Slack is used as the basic variable
                Q_ASSERT(c->helper.second == 1.0);
                setValueAt(i, 0, c->helper.first->index);
               setValueAt(i, c->helper.first->index, 1.0);
            }
    
           QHash<QSimplexVariable *, qreal>::const_iterator iter;
           for (iter = c->variables.constBegin();
                 iter != c->variables.constEnd();
                 ++iter) {
                setValueAt(i, iter.key()->index, iter.value());
           }
    
           setValueAt(i, columns - 1, c->constant);
        }
    
        // Set objective for the first-phase Simplex.
        // Z = -1 * sum_of_artificial_vars
        for (int j = firstArtificial; j < columns - 1; ++j)
            setValueAt(0, j, 1.0);
    
        // Maximize our objective (artificial vars go to zero)
        solveMaxHelper();
    
        // If there is a solution where the sum of all artificial
        // variables is zero, then all of them can be removed and yet
        // we will have a feasible (but not optimal) solution for the
        // original problem.
        // Otherwise, we clean up our structures and report there is
        // no feasible solution.
        if ((valueAt(0, columns - 1) != 0.0) && (qAbs(valueAt(0, columns - 1)) > 0.00001)) {
           qWarning() << "QSimplex: No feasible solution!";
            return false;
        }
    
        // Remove artificial variables. We already have a feasible
        // solution for the first problem, thus we don't need them
        // anymore.
        clearColumns(firstArtificial, columns - 2);
    
        return true;
    }
    
    /*!
      \internal
    
      Run simplex on the current matrix with the current objective.
    
      This is the iterative method. The matrix lines are combined
      as to modify the variable values towards the best solution possible.
      The method returns when the matrix is in the optimal state.
    */
    void QSimplex::solveMaxHelper()
    {
        reducedRowEchelon();
        while (iterate()) ;
    }
    
    /*!
      \internal
    */
    void QSimplex::setObjective(QSimplexConstraint *newObjective)
    {
        objective = newObjective;
    }
    
    /*!
      \internal
    */
    void QSimplex::clearRow(int rowIndex)
    {
        qreal *item = matrix + rowIndex * columns;
        for (int i = 0; i < columns; ++i)
            item[i] = 0.0;
    }
    
    /*!
      \internal
    */
    void QSimplex::clearColumns(int first, int last)
    {
        for (int i = 0; i < rows; ++i) {
            qreal *row = matrix + i * columns;
            for (int j = first; j <= last; ++j)
                row[j] = 0.0;
        }
    }
    
    /*!
      \internal
    */
    void QSimplex::dumpMatrix()
    {
       qDebug("---- Simplex Matrix ----\n");
    
        QString str(QLatin1String("       "));
        for (int j = 0; j < columns; ++j)
            str += QString::fromAscii("  <%1 >").arg(j, 2);
        qDebug("%s", qPrintable(str));
        for (int i = 0; i < rows; ++i) {
            str = QString::fromAscii("Row %1:").arg(i, 2);
    
            qreal *row = matrix + i * columns;
            for (int j = 0; j < columns; ++j)
                str += QString::fromAscii("%1").arg(row[j], 7, 'f', 2);
            qDebug("%s", qPrintable(str));
        }
        qDebug("------------------------\n");
    }
    
    /*!
      \internal
    */
    void QSimplex::combineRows(int toIndex, int fromIndex, qreal factor)
    {
        if (!factor)
            return;
    
        qreal *from = matrix + fromIndex * columns;
        qreal *to = matrix + toIndex * columns;
        for (int j = 1; j < columns; ++j) {
            qreal value = from[j];
    
            // skip to[j] = to[j] + factor*0.0
            if (value == 0.0)
                continue;
    
            to[j] += factor * value;
    
            // ### Avoid Numerical errors
            if (qAbs(to[j]) < 0.0000000001)
                to[j] = 0.0;
        }
    }
    
    /*!
      \internal
    */
    int QSimplex::findPivotColumn()
    {
        qreal min = 0;
        int minIndex = -1;
    
        for (int j = 0; j < columns-1; ++j) {
            if (valueAt(0, j) < min) {
                min = valueAt(0, j);
                minIndex = j;
            }
        }
    
        return minIndex;
    }
    
    /*!
      \internal
    
      For a given pivot column, find the pivot row. That is, the row with the
      minimum associated "quotient" where:
    
      - quotient is the division of the value in the last column by the value
        in the pivot column.
      - rows with value less or equal to zero are ignored
      - if two rows have the same quotient, lines are chosen based on the
        highest variable index (value in the first column)
    
      The last condition avoids a bug where artificial variables would be
      left behind for the second-phase simplex, and with 'good'
      constraints would be removed before it, what would lead to incorrect
      results.
    */
    int QSimplex::pivotRowForColumn(int column)
    {
        qreal min = qreal(999999999999.0); // ###
        int minIndex = -1;
    
        for (int i = 1; i < rows; ++i) {
            qreal divisor = valueAt(i, column);
            if (divisor <= 0)
                continue;
    
            qreal quotient = valueAt(i, columns - 1) / divisor;
            if (quotient < min) {
                min = quotient;
                minIndex = i;
            } else if ((quotient == min) && (valueAt(i, 0) > valueAt(minIndex, 0))) {
                minIndex = i;
            }
        }
    
        return minIndex;
    }
    
    /*!
      \internal
    */
    void QSimplex::reducedRowEchelon()
    {
        for (int i = 1; i < rows; ++i) {
            int factorInObjectiveRow = valueAt(i, 0);
            combineRows(0, i, -1 * valueAt(0, factorInObjectiveRow));
        }
    }
    
    /*!
      \internal
    
      Does one iteration towards a better solution for the problem.
      See 'solveMaxHelper'.
    */
    bool QSimplex::iterate()
    {
        // Find Pivot column
        int pivotColumn = findPivotColumn();
        if (pivotColumn == -1)
            return false;
    
        // Find Pivot row for column
        int pivotRow = pivotRowForColumn(pivotColumn);
        if (pivotRow == -1) {
            qWarning() << "QSimplex: Unbounded problem!";
            return false;
        }
    
        // Normalize Pivot Row
        qreal pivot = valueAt(pivotRow, pivotColumn);
        if (pivot != 1.0)
            combineRows(pivotRow, pivotRow, (qreal(1.0) - pivot) / pivot);
    
        // Update other rows
        for (int row=0; row < rows; ++row) {
          if (row == pivotRow)
                continue;
    
            combineRows(row, pivotRow, -1 * valueAt(row, pivotColumn));
        }
        // Update first column
        setValueAt(pivotRow, 0, pivotColumn);
    
        //    dumpMatrix();
        //    qDebug("------------ end of iteration --------------\n");
        return true;
    }
    
    /*!
      \internal
    
      Both solveMin and solveMax are interfaces to this method.
    
      The enum solverFactor admits 2 values: Minimum (-1) and Maximum (+1).
    
      This method sets the original objective and runs the second phase
      Simplex to obtain the optimal solution for the problem. As the internal
      simplex solver is only able to _maximize_ objectives, we handle the
      minimization case by inverting the original objective and then
      maximizing it.
    */
    qreal QSimplex::solver(solverFactor factor)
    {
        // Remove old objective
        clearRow(0);
    
        // Set new objective in the first row of the simplex matrix
        qreal resultOffset = 0;
        QHash<QSimplexVariable *, qreal>::const_iterator iter;
        for (iter = objective->variables.constBegin();
             iter != objective->variables.constEnd();
             ++iter) {
    
            // Check if the variable was removed in the simplification process.
            // If so, we save its offset to the objective function and skip adding
            // it to the matrix.
            if (iter.key()->index == -1) {
                resultOffset += iter.value() * iter.key()->result;
                continue;
            }
    
            setValueAt(0, iter.key()->index, -1 * factor * iter.value());
      }
        solveMaxHelper();
        collectResults();
    
    #ifdef QT_DEBUG
        for (int i = 0; i < constraints.size(); ++i) {
            Q_ASSERT(constraints[i]->isSatisfied());
       }
    #endif
    
       // Return the value calculated by the simplex plus the value of the
        // fixed variables.
        return (factor * valueAt(0, columns - 1)) + resultOffset;
    }
    
    /*!
      \internal
      Minimize the original objective.
    */
    qreal QSimplex::solveMin()
    {
        return solver(Minimum);
    
    
    /*!
      \internal
      Maximize the original objective.
    */
    qreal QSimplex::solveMax()
    {
        return solver(Maximum);
    }
    /*!
      \internal
    
      Reads results from the simplified matrix and saves them in the
      "result" member of each QSimplexVariable.
    */
    void QSimplex::collectResults()
    {
        // All variables are zero unless overridden below.
    
        // ### Is this really needed? Is there any chance that an
        // important variable remains as non-basic at the end of simplex?
        for (int i = 0; i < variables.size(); ++i)
            variables[i]->result = 0;
    
        // Basic variables
        // Update the variable indicated in the first column with the value
        // in the last column.
        for (int i = 1; i < rows; ++i) {
            int index = valueAt(i, 0) - 1;
            if (index < variables.size())
                variables[index]->result = valueAt(i, columns - 1);
        }
    }
    
    /*!
      \internal
    
      Looks for single-valued variables and remove them from the constraints list.
    */
    bool QSimplex::simplifyConstraints(QList<QSimplexConstraint *> *constraints)
    {
        QHash<QSimplexVariable *, qreal> results;   // List of single-valued variables
        bool modified = true;                       // Any chance more optimization exists?
    
        while (modified) {
            modified = false;
    
            // For all constraints
            QList<QSimplexConstraint *>::iterator iter = constraints->begin();
            while (iter != constraints->end()) {
                QSimplexConstraint *c = *iter;
                if ((c->ratio == QSimplexConstraint::Equal) && (c->variables.count() == 1)) {
                    // Check whether this is a constraint of type Var == K
                    // If so, save its value to "results".
                    QSimplexVariable *variable = c->variables.constBegin().key();
                    qreal result = c->constant / c->variables.value(variable);
    
                    results.insert(variable, result);
                    variable->result = result;
                   variable->index = -1;
                    modified = true;
    
                }
    
                // Replace known values among their variables
                QHash<QSimplexVariable *, qreal>::const_iterator r;
                for (r = results.constBegin(); r != results.constEnd(); ++r) {
                    if (c->variables.contains(r.key())) {
                        c->constant -= r.value() * c->variables.take(r.key());
                        modified = true;
                   }
                }
    
                // Keep it normalized
                if (c->constant < 0)
                    c->invert();
    
                if (c->variables.isEmpty()) {
                    // If constraint became empty due to substitution, delete it.
                    if (c->isSatisfied() == false)
                        // We must ensure that the constraint soon to be deleted would not
                        // make the problem unfeasible if left behind. If that's the case,
                        // we return false so the simplex solver can properly report that.
                        return false;
    
                    delete c;
                    iter = constraints->erase(iter);
                } else {
                    ++iter;
                }
            }
        }
    
        return true;
    }
    
    void QSimplexConstraint::invert()
    {
        constant = -constant;
        ratio = Ratio(2 - ratio);
    
        QHash<QSimplexVariable *, qreal>::iterator iter;
        for (iter = variables.begin(); iter != variables.end(); ++iter) {
            iter.value() = -iter.value();
        }
    }
    QT_END_NAMESPACE

  9. #9
    Join Date
    Sep 2013
    Posts
    36

    Re: help with implementation of Linear programming solver.

    Hi paul,

    I have not seen the examples of how to use this LP solver. It would have been easier for me to understand if they had give me an example of implementation of this class. Any ways, Please give me some time . I will post some examples if I am able to find any.

    Regards,
    Naveen.

  10. #10
    Join Date
    Sep 2013
    Posts
    36

    Re: help with implementation of Linear programming solver.

    Hi paul,
    I could not find any examples of implementation. My logic is to convert the string (lines that I parsed) into struct Qsimplexconstraint and then call the functions available in the LPsolver. I need to iterate this to parse all the lines of the file. Atlast, I should call collectResults() function to get the results. Is my logic correct?

    Regards,
    Naveen.

  11. #11
    Join Date
    Apr 1999
    Posts
    27,426

    Re: help with implementation of Linear programming solver.

    Quote Originally Posted by navy1991 View Post
    Hi paul,
    I could not find any examples of implementation. My logic is to convert the string (lines that I parsed) into struct Qsimplexconstraint
    All you need to do first is this:
    Code:
    #include <iostream>
    #include <string>
    #include <map>
    
    typedef double qreal;
    
    struct QSimplexVariable
    {
        QSimplexVariable() : result(0), index(0) {}
        qreal result;
        int index;
    };
    
    struct QSimplexConstraint
    {
       QSimplexConstraint() : constant(0), ratio(Equal), artificial(0) {}
       enum Ratio {LessOrEqual = 0,Equal,MoreOrEqual};
       float constant;
       Ratio ratio;
       std::pair<QSimplexVariable *, qreal> helper;
       QSimplexVariable * artificial;
    };
    
    
    QSimplexConstraint parseMe(const std::string& str)
    {
        // you fill this in and return a QSimplexConstraint    
    }
    
    using namespace std;
    
    int main()
    {
       std::string constraintString;
       cout << "Enter a constraint string: ";
       getline(cin, constraintString);     
       QSimplexConstraint qc = parseMe( constraintString );
       //...
    }
    Since I don't know anything about QTPair or qreal (these are not C++ types), I used standard types in place of them. Note that the entire code above compiles without any knowledge of Qt or anything to do with linear programming. The problem is now focused on your original question, regardless of all of the other stuff.

    So your job is to take that string, call parseMe, and then return a QSimplexConstraint that contains the proper information. Once you have the simple case working, then and only then do you now consider file I/O and what to do next.

    Regards,

    Paul McKenzie
    Last edited by Paul McKenzie; September 18th, 2013 at 01:20 PM.

  12. #12
    Join Date
    Sep 2013
    Posts
    36

    Re: help with implementation of Linear programming solver.

    Hi Paul,

    Code:
    QSimplexConstraint parseMe(const std::string& str)
    {
        // you fill this in and return a QSimplexConstraint    
    }
    What is this parseMe() function ? I tried to search about it in google but could not learn much.

    The data in the file is of the format :



    c1 * X1 + c2 * X2 + ... = K
    or <= K
    or >= K

    Where (ci, Xi) are the pairs in "variables" and K the real "constant".

    Now I want to parse

    constant = K (RHS)
    ratio = symbol in between RHS and LHS
    QsimplexVariable = (ci,Xi) (LHS)

    Could you please tell me how I can do this and get a QSimplexConstraint?

    Regards,
    Naveen.

  13. #13
    Join Date
    Sep 2013
    Posts
    36

    Re: help with implementation of Linear programming solver.

    Hi Paul,
    As you have mentioned , I implemented your code as a stand alone program and it worked perfectly fine. In the place of " // you fill this in and return a QSimplexConstraint " I wrote "QSimplexConstraint() " . I got a message : "Enter a constraint string: " . After this , I entered a constraint like "+ a -b = 0" . It is similar to the constraint in the file.After that , the console window disappeared. I am thinking how should I use the parseMe() function. Unfortunately, I have not yet found a way.

    Regards,
    Naveen.

  14. #14
    Join Date
    Apr 1999
    Posts
    27,426

    Re: help with implementation of Linear programming solver.

    Quote Originally Posted by navy1991 View Post
    What is this parseMe() function ?
    I just made it up.

    I don't think you're understanding the general point. The issue is that you have a string, and now you want to take this string and turn it into something else, i.e. transform the string into something that fits into the problem domain.

    For example, if the string represented a mathematical expression, the probable goal is to take that string, parse it, and place the tokens in some sort of syntax tree. If that string represented English words, you would parse the string and maybe place the indivdual words in a table, etc.

    All that other stuff concerning linear programming is moot and not important.
    I tried to search about it in google but could not learn much.
    No. The goal is for you to drop everything and concentrate on taking that string, doing soemthing to it, and returning a QSimplexConstraint based on that string. If you can't do that, then there is no need to write the rest of the program.
    The data in the file is of the format :
    Forget about the file. It doesn't matter where the data comes from, and introducing file I/O is a waste of your time at this point. For one, you now have to spend precious hours making sure the file reading work. What good is reading from a file when you can't do a basic parsing of the string?

    Just enter a string that reprsents one line of data. Then take that entered string and see if you can create a QSimplexConstraint from that string. That is how any programmer would start out -- I think your problem is that you want to complete this whole assignment in one shot, and not unit test or develop each piece.
    Code:
    The data in the file is of the format :
    
    c1 * X1 + c2 * X2 + ... = K
    or <= K
    or >= K
    Please post the actual file, not an interpretation. Are you saying that the word "or" is in the file?
    Could you please tell me how I can do this and get a QSimplexConstraint?
    First post the real file, not a mock up. Second, parsing input may require you to first define in a general sense the rules of the syntax for the data. Then you may have to write a recursive descent parser to parse the data correctly, or if the rules are simple, a simple "linear" parser.

    No one can just say "use this code to parse the file" -- that is impossible. You have to ascertain what constitutes a data line, and then assess what type of parser is needed to parse the data and store the data into the program's data structures.

    Regards,

    Paul McKenzie

  15. #15
    Join Date
    Apr 1999
    Posts
    27,426

    Re: help with implementation of Linear programming solver.

    Quote Originally Posted by navy1991 View Post
    I am thinking how should I use the parseMe() function. Unfortunately, I have not yet found a way.
    That is because there is no way to do it by calling just a function. Parsing expressions is not trivial -- you have to know the general syntax of what to expect, probably then write a grammar that matches that syntax using production rules, and then write the recursive parser. I don't know if you knew that it isn't trivial or "automatic" to do these things.

    Just because it looks easy by eyesight doesn't mean it's easy to write a program for what you see. Can you write an algebraic problem solver? We can solve algrebaic problems easy "by hand", but how hard is it to write a program to do it? Very difficult.

    So just because it "looks easy" from a non-programmer perspective to interpret that data, it isn't easy to write a program to do what we can do easily with pencil and paper.

    Regards,

    Paul McKenzie

Page 1 of 3 123 LastLast

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •  


Azure Activities Information Page

Windows Mobile Development Center


Click Here to Expand Forum to Full Width

This is a CodeGuru survey question.


Featured


HTML5 Development Center