#include <glpk.h>
#include "ei_simplex.h"

EISimplex* ei_simplex_new() {

    EISimplex* sx = (EISimplex*)malloc(sizeof(EISimplex));
    assert(sx!=NULL);
    
    sx->lp_problem=glp_create_prob();
    assert(sx->lp_problem!=NULL);

    return sx;

}

void ei_simplex_delete(EISimplex* sx) {

    if (sx!=NULL) {
        if (sx->lp_problem!=NULL) glp_delete_prob(sx->lp_problem);
    }
    free(sx);
}


int ei_simplex_solve_scp(EISimplex* sx,SCPView* view,int itrLimit) {

    int result;

    glp_smcp smcp;
    

    /*Check parameters and conditions*/
    assert(sx!=NULL && view!=NULL);

    /*Simplex method parameters*/
    glp_init_smcp(&smcp);
    smcp.msg_lev=GLP_MSG_OFF;
    if (itrLimit>0) smcp.it_lim=itrLimit;

    /*Just plain minimalization scp*/
    ei_simplex_formulate_scp(sx,view,1);

    /*1st stage real numbers*/
    if(glp_simplex(sx->lp_problem,&smcp)==GLP_EITLIM) return EI_SIMPLEX_OVERRUN;
    result=glp_get_status(sx->lp_problem);

    /*2nd stage rational numbers*/
    /*if (result!=GLP_OPT) {*/
        if(glp_exact(sx->lp_problem,&smcp)==GLP_EITLIM) return EI_SIMPLEX_OVERRUN;
    /*}*/

    /*What is the quality of the solution. If no optimum found,
      we report infeasible state*/
    result=glp_get_status(sx->lp_problem);
    if (result!=GLP_OPT) return EI_SIMPLEX_INFEAS;
    return EI_SIMPLEX_OPTIMUM;
    
}

int ei_simplex_solve_lexmin(EISimplex* sx,SCPView* view,int maxCost,int auxCols) {

    int i,j,primalState;
    glp_smcp smcp;
    int retVal;
    

    double* strip=NULL;
    int* indexer=NULL;

    /*Check*/
    assert(sx!=NULL);
    assert(view!=NULL);
    assert(auxCols<=view->columns);

    /*NOTE: GLPK uses indices style 1..n instead of 0..n-1
     *
     *      If LP problem formulation is
     *      confusing, export formulated LP problem in some
     *      standard form and examine it. Export all problems
     *      in the loop below before glp_simplex() call
     */

    strip=(double*)malloc((view->columns+1)*sizeof(double));
    indexer=(int*)malloc((view->columns+1)*sizeof(int));

    assert(strip!=NULL && indexer!=NULL);
    

    glp_init_smcp(&smcp);
    
    

    /*Formulate lp_problem for SCP, this is a base*/
    ei_simplex_formulate_scp(sx,view,0);

       /*Indexer*/
        for(j=1;j<view->columns+1;j++) indexer[j]=j;
        /*Strip of costs*/
        for(j=0;j<view->columns;j++) strip[j+1]=scp_view_get_cost(view,j);
    
        /*There is also requirement for cost, so we add a new constraint
        *row that will limit maximum cost
        */
        glp_add_rows(sx->lp_problem,1);
        glp_set_mat_row(sx->lp_problem,glp_get_num_rows(sx->lp_problem),view->columns,indexer,strip);
        glp_set_row_bnds(sx->lp_problem,glp_get_num_rows(sx->lp_problem),GLP_UP,0,maxCost);
   

  
    /*If there is a necessity for auxiliary columns, add row representing
     *auxiliary constraint. This is used for group testing heuristic inside
     *the LCE algorithm. This ensures that at least one of first auxCols
     *has to be part of the solution
     */
    if (auxCols>0) {
        glp_add_rows(sx->lp_problem,1);
        for(j=0;j<auxCols;j++) {
            strip[j+1]=1;
        }
        glp_set_mat_row(sx->lp_problem,glp_get_num_rows(sx->lp_problem),auxCols,indexer,strip);
        glp_set_row_bnds(sx->lp_problem,glp_get_num_rows(sx->lp_problem),GLP_LO,1,0);
    }



    /*Clear objective function*/
    for(j=0;j<view->columns;j++) glp_set_obj_coef(sx->lp_problem,j+1,0);


    smcp.meth=GLP_PRIMAL;
    smcp.presolve=GLP_OFF;
    smcp.msg_lev=GLP_MSG_OFF;
    smcp.tol_bnd=0.000000000000001;
    smcp.tol_piv=0.000000000000001;
    smcp.tol_dj= 0.000000000000001;
    
    /*Performing succesive problem solving, while increasing number
      of fixed variables*/
    for (i=0;i<view->columns;i++) {

        

        /*In i-th iteration, we are trying to minimize i-th column*/
        glp_set_obj_coef(sx->lp_problem,i+1,1);

        /*NOTE: this is a place for occassional LP problem export*/
        
        /*First stage, 120 iterations in real numbers*/
        smcp.it_lim=160;
        glp_simplex(sx->lp_problem,&smcp);
        glp_scale_prob(sx->lp_problem,GLP_SF_AUTO);
        glp_simplex(sx->lp_problem,&smcp);

        /*Second stage, rational numbers*/
        smcp.it_lim=PS_INFINITY;
        glp_scale_prob(sx->lp_problem,GLP_SF_AUTO);
        glp_exact(sx->lp_problem,&smcp);
        
        
        
        /*What was solved, if some*/
        primalState=glp_get_status(sx->lp_problem);
        

        /*Not possible to solve with current set of variables, cancelling*/
        if (primalState!=GLP_OPT) {
            retVal=EI_SIMPLEX_INFEAS;
            goto f_end;
        }

        /*In other case, we contine with iteration
         *Moving to next column
         */
        glp_set_obj_coef(sx->lp_problem,i+1,0);

        /*I-th variable is now fixed*/
        glp_set_col_bnds(sx->lp_problem,i+1,GLP_FX,glp_get_col_prim(sx->lp_problem,i+1),0);

        /*Continuing to next iteration*/
    }
    retVal=EI_SIMPLEX_OPTIMUM;

    f_end:
    free(strip);
    free(indexer);

    return retVal;
    
    
}



void ei_simplex_formulate_scp(EISimplex* sx,SCPView* view,int withObjective) {
    
    int i=0;
    int j=0;

    double* strip=NULL;
    int* indexer=NULL;

    strip=(double*)malloc((view->columns+1)*sizeof(double));
    indexer=( int*)malloc((view->columns+1)*sizeof(int));

    assert(strip!=NULL && indexer!=NULL);

    for(i=0;i<view->columns;i++) {
        indexer[i+1]=i+1;
    }

    /*Set number of rows and columns in the  matrix*/
    glp_add_cols(sx->lp_problem,view->columns);
    glp_add_rows(sx->lp_problem,view->rows);

    /*Set objective function coeficients according to costs of the columns*/
    if (withObjective==1) {
        for(j=0;j<view->columns;j++) glp_set_obj_coef(sx->lp_problem,j+1,scp_view_get_cost(view,j));
    }

    /*Fill the main matrix*/
    for (i=0;i<view->rows;i++) {
        /*Fill the row strip with 0 and 1*/
        for(j=0;j<view->columns;j++) {
            strip[j+1]=scp_view_get(view,i,j);
        }
        /*Set the row*/
        glp_set_mat_row(sx->lp_problem,i+1,view->columns,indexer,strip);
        /*This row has to be >=1*/
        glp_set_row_bnds(sx->lp_problem,i+1,GLP_LO,1,0);
    }
    /*All columns have to be >=0*/
    for(j=0;j<view->columns;j++) glp_set_col_bnds(sx->lp_problem,j+1,GLP_LO,0,0);

    /*Minimization is always performed*/
    glp_set_obj_dir(sx->lp_problem,GLP_MIN);

    free(strip);
    free(indexer);

}


void ei_simplex_fill_result(EISimplex* sx,double* retVector) {

    int j=0;
    int cols=0;

    /*Check parameter*/
    assert(retVector!=NULL && sx!=NULL);

    cols=glp_get_num_cols(sx->lp_problem);

    for(j=0;j<cols;j++) {
        retVector[j]=glp_get_col_prim(sx->lp_problem,j+1);
    }
}

void ei_simplex_get_rhs(EISimplex* sx, int column, double* val) {
    assert(sx!=NULL && val!=NULL);
    *val=glp_get_col_prim(sx->lp_problem,column+1);
}

double ei_simplex_get_cost(EISimplex* sx) {
    assert(sx!=NULL);
    return glp_get_obj_val(sx->lp_problem);
}

