/*
 * A program to perform Data Envelopment Analysis. The problem
 * formulation is read from a GMPL model in the file dea.mod.
 * This program solves the LP for each DMU and collects the resulting
 * information
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <glpk.h>

/* given a string of the form "inputs[space]", return a string of the form "space" */
char *pretty(char *name) {
  char *beg, *end;
  beg = strchr(name, '[');
  end = strchr(name, ']');
  if (beg == NULL || end == NULL) 
    return NULL;
  else {
    *end = '\0';
    return beg+1;
  }
}

int main(void) { 
  glp_prob *lp;
  glp_tran *tran;
  int ret, ret1;
  int i, j, colw, n, m;
  int *ind, *ind1;
  double *val, *val1;
  int num_outputs;
  char symbol[50]; 
  char *sep = ",";
  FILE *of;

  lp = glp_create_prob();
  tran = glp_mpl_alloc_wksp();
  ret = glp_mpl_read_model(tran, "dea.mod", 1);
  if (ret != 0) {
    fprintf(stderr, "Error on translating model\n");
    goto skip;
  }
  ret = glp_mpl_read_data(tran, "totpt.dat");
  if (ret != 0) {
    fprintf(stderr, "Error on translating data\n");
    goto skip;
  }
  ret = glp_mpl_generate(tran, NULL);
  if (ret != 0) {
    fprintf(stderr, "Error on generating model\n");
    goto skip;
  }
  glp_mpl_build_prob(tran, lp);

  ret = glp_write_prob(lp, 0, "dea.txt");  /* write the problem in GLPK LP/MIP format */
  if (ret != 0)
    fprintf(stderr, "Error when writing problem\n");

  n  = glp_get_num_rows(lp);
  m  = glp_get_num_cols(lp);

  ind  = calloc(n, sizeof(int));
  val  = calloc(n, sizeof(double));
  ind1 = calloc(n, sizeof(int));
  val1 = calloc(n, sizeof(double));

  glp_create_index(lp);
  colw = glp_find_col(lp, "w");                  /* column number for column w */
  ret1 = glp_get_mat_col(lp, colw, ind1, val1);  /* store the column for w */
  num_outputs = ret1 - 1;                        /* there is one nonzero entry for w itself */


  /* open the output file and print the header */
  of = fopen("totpt.csv", "w");
  fprintf(of, "dmu"); fprintf(of, sep);
  fprintf(of, "efficiency"); fprintf(of, sep);
  for (i = 1; i <= n; ++i) {  
    if (glp_get_row_type(lp, i)==GLP_LO || glp_get_row_type(lp, i)==GLP_UP) {
      strcpy(symbol, glp_get_row_name(lp, i));
      fprintf(of, "%s", pretty(symbol));
      if (i < n) fprintf(of, sep);
    }
  }
  fprintf(of, "\n");

  /* solve the LP for each DMU */
  for (j = 1; j <= m; ++j) {
    if (strcmp(glp_get_col_name(lp, j), "w")==0) continue;

    ret = glp_get_mat_col(lp, j, ind, val);       /* store the column for the dmu represented by column j */
    for (i = 1; i <= num_outputs; ++i)
      val1[i] = -val[i];
    glp_set_mat_col(lp, colw, ret1, ind1, val1);  /* set the new output coefficients in the column for w */

    /* set the new RHS, i.e. the new input coefficients */
    for (i = 1; i <= n; ++i) {
      if (glp_get_row_type(lp, i) == GLP_UP)                /* only want rows corresponding to inputs */
	glp_set_row_bnds(lp, i, GLP_UP, 0.0, val[ind[i]]);  /* set the new RHS for the input constraints */
    }

    glp_simplex(lp, NULL);

    strcpy(symbol,glp_get_col_name(lp, j));
    fprintf(of, "%s", pretty(symbol));         /* DMU name */
    fprintf(of, sep);
    fprintf(of, "%f", 1/glp_get_obj_val(lp));  /* efficiency number */
    fprintf(of, sep);
    for (i = 1; i <= n; ++i) {       /* weights obtained as the dual values */
      if (glp_get_row_type(lp, i)==GLP_UP) { 
	fprintf(of, "%f", glp_get_row_dual(lp, i));
	if (i < n) fprintf(of, sep);
      }
      if (glp_get_row_type(lp, i)==GLP_LO) {  /* a nonzero dual value will be negative here */
	fprintf(of, "%f", -glp_get_row_dual(lp, i));
	if (i < n) fprintf(of, sep);
      }
    }
    fprintf(of, "\n");
  }
  fclose(of);

 skip: glp_mpl_free_wksp(tran);
  glp_delete_prob(lp);
  return 0;
}