#include <stdio.h>
#include <stdlib.h>
#include <glpk.h>
#include <float.h>
#include <string.h>
#include <assert.h>

#define MAX( a, b ) ( a > b ) ? ( a ) : ( b )

#define WEIGHT_SLACK_VARIABLE 1000

/** Input Data */


int L   = 250;  // bar size

int m   = 4;                      // nr. of orders
int w[] = { 187, 119, 74, 90 };   // size of bar in order i
int b[] = {   1,   2,  2,  1 };   // number of bars in order i

/* vector of dual variables */
double *pi = NULL;

/* GLPK linear program */
glp_prob *lp = NULL;

/** solves the knapsack problem where repetition of items are allowed
    cap: capacity    n: number o items
    weight[0,...n-1]: vector of weights
    profit[0,...n-1]: vector of profits
    selected[0,...n-1]: solution will be filled here, where
    selected[i] contains an integer indicating how many times item i was included
    work[] is a vector of size (cap+1)*2 where dynamic programming computation will occur **/
int solve_knapsack( int cap, int n, int weight[], int profit[],
                    int selected[], int work[] );


/** creates constraints and adds
    artificial variables ***/
void fill_initial_LP();

int main( int argc, char **argv )
{
    int exit_code = EXIT_SUCCESS;

    lp = glp_create_prob();

    pi = malloc( sizeof(double)*(m+1) );

    int iteration = 0;

    /// dynamic programming related memory structures
    int *dpWork = NULL;
    dpWork = malloc( sizeof(int)*((L+1)*2) );
    int *selected = NULL;
    selected = malloc( sizeof(int)*m );
    int *ipi = NULL;
    ipi = malloc( sizeof(int)*m );
    int nColsLambda = 0;

    int *colIdx = NULL;
    double *colCoef = NULL;
    colIdx = malloc( sizeof(int)*(m+1) );
    colCoef = malloc( sizeof(double)*(m+1) );

    printf("Cutting Stock Solver.\n\n");

    fill_initial_LP();

    int newCols;

SOLVE_RELAXATION:
    ++iteration;
    printf("-> on iteration %d\n\n", iteration );
    glp_simplex( lp, NULL );

    if ( glp_get_status( lp ) != GLP_OPT )
    {
        fprintf( stderr, "ERROR: No optimal solution obtained, check problem data.\n" );
        exit_code = EXIT_FAILURE;
        goto TERMINATE;
    }

    /* getting dual values */
    int i;
    for ( i=1 ; (i<=m) ; ++i )
        pi[i] = glp_get_row_dual( lp, i );

    printf("\npi vector - dual variables for constraints: \n\t");
    for ( i=1 ; (i<=m) ; ++i )
        printf("[%d] %g ", i, pi[i] );
    printf("\n\n");

    /** solving the pricing problem ***/
    for ( i=0 ; (i<m) ; ++i )            // knapsack solver needs
        ipi[i] = ( pi[i+1] * 1000 );     // integer values

    double rc = 1.0 - (((double)( solve_knapsack( L, m, w, ipi, selected, dpWork ) )) / 1000.0);
    newCols = 0;
    if ( rc < -0.000001 )
    {
        newCols++;
        printf("\tattractive column rc:(%g)\n", rc );
        printf("\t\tselected: ");
        int k;
        for ( k=0 ; (k<m) ; ++k )
            if (selected[k])
                printf("%d(%d) ", w[k], selected[k] );
        printf("\n");

        int nSizes = 0;
        for ( k=0 ; (k<m) ; ++k )
            if (selected[k])
            {
                nSizes++;
                colIdx[ nSizes ] = k+1;
                colCoef[ nSizes ] = selected[k];
            }


        /// inserting column
        glp_add_cols( lp, 1 );
        int col = glp_get_num_cols( lp );
        nColsLambda++;
        char colName[ 256 ];
        sprintf( colName, "lbda[%d]", nColsLambda );
        glp_set_col_name( lp, col, colName );
        glp_set_col_bnds( lp, col, GLP_DB, 0.0, 1.0 );
        glp_set_obj_coef( lp, col, 1 );
        glp_set_mat_col( lp, col, nSizes, colIdx, colCoef );
    }

    char pNameIt[256];
    sprintf( pNameIt, "cg_%d.lp", iteration );
    glp_write_lp( lp, NULL, pNameIt );

    if (newCols)
    {
        printf("\nPress Enter to proceed to the next iteration\n");
        getchar();

        printf("New columns added. Starting a new iteration.\n\n");
        goto SOLVE_RELAXATION;
    }

    /// solution
    double objVal = glp_get_obj_val( lp );
    printf("solution found with cost %g\n", objVal );
    int nPattern = 0, j;

    for ( j=1 ; (j<=glp_get_num_cols(lp)) ; ++j )
        if ( ( glp_get_col_prim( lp, j ) >= 0.01 ) && ( strstr( glp_get_col_name( lp, j ), "lbda" ) ) )
        {
            printf(" %s  %g \n", glp_get_col_name(lp, j), glp_get_col_prim( lp, j ) );
            int nzCol = glp_get_mat_col( lp, j, colIdx-1, colCoef-1 );
            int k;
            printf("\t");
            for ( k=0 ; k<nzCol ; ++k )
                printf( "%d(%d) ", w[colIdx[k]-1], (int)colCoef[k] );
            printf("\n");

            nPattern++;
        }

    printf("\n");

    glp_free( lp );

    glp_free_env();

TERMINATE:
    if (pi)
        free( pi );
    if (dpWork)
        free(dpWork);
    if (selected)
        free( selected );
    if (ipi)
        free( ipi );

    if (colIdx)
        free(colIdx);
    if (colCoef)
        free(colCoef);

    return exit_code;
}

void fill_initial_LP()
{
    glp_set_prob_name( lp, "cutting stock" );

    int i;

    /// artificial variables, one for each
    /// constraint which will be created after
    glp_add_cols( lp, m );
    for ( i=1 ; (i<=m) ; ++i )
    {
        char name[64];
        sprintf( name, "art_%d", i );
        glp_set_col_name( lp, i, name );
        glp_set_col_bnds( lp, i, GLP_DB, 0.0, 1.0 );
        glp_set_obj_coef( lp, i, WEIGHT_SLACK_VARIABLE );
    }

    /// creating one constraint for each size
    glp_add_rows( lp, m );
    for ( i=1 ; (i<=m) ; ++i )
    {
        int idx[] = { 0, i };
        double coef[] = { 0.0, b[i-1] };
        glp_set_mat_row( lp, i, 1, idx, coef );
        glp_set_row_bnds( lp, i, GLP_LO, b[i-1], DBL_MAX );
    }
}

int solve_knapsack( int cap, int n, int weight[], int profit[],
                    int selected[], int work[] )
{
    int bestProfit = -1;
    /// best profit within a given capacity
    int *bestProfitWC = work;
    int *bestItemToIncludeWC =  work + (cap+1);
    memset( bestProfitWC, 0, sizeof(int)*(cap+1) );
    memset( selected, 0, sizeof(int)*n );

    int i;
    for ( i=0 ; (i<=cap) ; ++i )
        bestItemToIncludeWC[i] = -1;

    int lastUsefullCapacity = -1;

    int curCap;
    for ( curCap=1 ; (curCap<=cap) ; ++curCap )
    {
        for ( i=0 ; (i<n) ; ++i )
        {
            int profitWith = -1;
            int rCap = curCap - weight[ i ];
            if ( rCap >= 0 )
                profitWith = profit[i] + bestProfitWC[rCap];

            if ( profitWith > bestProfitWC[curCap] )
            {
                bestProfitWC[curCap] = profitWith;
                bestItemToIncludeWC[curCap] = i;
                lastUsefullCapacity = curCap;
            }
        }
    }

    if ( lastUsefullCapacity == -1 )
        return 0;

    bestProfit = bestProfitWC[lastUsefullCapacity];

    curCap = lastUsefullCapacity;
    while ( curCap >0 )
    {
        int curItem = bestItemToIncludeWC[curCap];
        if ( curItem == -1 )
            return bestProfit;
        selected[ curItem ]++;
        curCap -= weight[ curItem ];
    }

    return bestProfit;
}