#include <stdio.h>

#include "ampl/ampl_c.h"

int main(int argc, char **argv) {
  // Create an AMPL instance
  AMPL *ampl;
  AMPL_CALL(AMPL_Create(&ampl));

  /*
  // If the AMPL installation directory is not in the system search path:
  AMPL_ENVIRONMENT *env;
  AMPL_EnvironmentCreate(&env, "full path to the AMPL installation directory", "");
  AMPL_CALL(AMPL_CreateWithEnv(&ampl, env));
  */

  // Number of steps of the efficient frontier
  const int steps = 10;

  if (argc > 1) AMPL_CALL(AMPL_SetOption(ampl, "solver", argv[1]));

  AMPL_CALL(AMPL_SetBoolOption(ampl, "reset_initial_guesses", true));
  AMPL_CALL(AMPL_SetBoolOption(ampl, "send_statuses", false));
  AMPL_CALL(AMPL_SetOption(ampl, "solver", "cplex"));

  // Load the AMPL model from file
  const char *modelDirectory = argc == 3 ? argv[2] : "../models";
  char dir[256];
  char mod[256];
  char run[256];
  sprintf(dir, "%s%s", modelDirectory, "/qpmv");
  sprintf(mod, "%s%s", modelDirectory, "/qpmv/qpmv.mod");
  sprintf(run, "%s%s", modelDirectory, "/qpmv/qpmvbit.run");

  AMPL_CALL(AMPL_Read(ampl, mod));
  AMPL_CALL(AMPL_Read(ampl, run));

  // Set tables directory (parameter used in the script above)
  AMPL_CALL(AMPL_ParameterSetString(ampl, "data_dir", dir));

  // Read tables
  AMPL_CALL(AMPL_ReadTable(ampl, "assetstable"));
  AMPL_CALL(AMPL_ReadTable(ampl, "astrets"));

  // Relax the integrality
  AMPL_CALL(AMPL_SetBoolOption(ampl, "relax_integrality", true));

  // Solve the problem
  AMPL_CALL(AMPL_Solve(ampl, "", ""));

  // Calibrate the efficient frontier range
  AMPL_DATAFRAME *values;
  double minret;
  size_t index;
  size_t nrows;

  AMPL_CALL(AMPL_VariableGetValue(ampl, "portret", &minret));
  AMPL_CALL(AMPL_EntityGetValues(ampl, "averret", NULL, 0, &values));
  AMPL_CALL(AMPL_DataFrameGetColumnIndex(values, "averret", &index));
  AMPL_CALL(AMPL_DataFrameGetNumRows(values, &nrows));

  double maxret = 0;

  for (size_t i = 0; i < nrows; i++) {
    AMPL_VARIANT *value;
    double value_dbl;

    AMPL_CALL(AMPL_DataFrameElement(values, i, index, &value));
    AMPL_VariantGetNumericValue(value, &value_dbl);
    if (value_dbl > maxret)
      maxret = value_dbl;
  }
  AMPL_DataFrameFree(&values);

  double stepsize = (maxret - minret) / steps;
  double returns[steps];
  double variances[steps];

  for (int i = 0; i < steps; i++) {
    printf("Solving for return = %g\n", maxret - (i - 1) * stepsize);

    // Set target return to the desired point
    AMPL_CALL(AMPL_ParameterSetNumeric(ampl, "targetret", maxret - (i - 1) * stepsize));
    AMPL_CALL(AMPL_Eval(ampl, "let stockopall:={};let stockrun:=stockall;"));

    // Relax integrality
    AMPL_CALL(AMPL_SetDblOption(ampl, "relax_integrality", 1));
    AMPL_CALL(AMPL_Solve(ampl, "", ""));
    double objval;
    AMPL_CALL(AMPL_GetValueNumeric(ampl, "cst", &objval));
    printf("QP result = %g\n", objval);

    // Adjust included stocks
    AMPL_CALL(AMPL_Eval(ampl, "let stockrun:={i in stockrun:weights[i]>0};"));
    AMPL_CALL(AMPL_Eval(ampl, "let stockopall:={i in stockrun:weights[i]>0.5};"));

    // Set integrality back
    AMPL_CALL(AMPL_SetBoolOption(ampl, "relax_integrality", false));
    AMPL_CALL(AMPL_Solve(ampl, "", ""));
    AMPL_CALL(AMPL_GetValueNumeric(ampl, "cst", &objval));
    printf("QMIP result = %g\n", objval);

    // Store data of corrent frontier point
    returns[i] = maxret - (i - 1) * stepsize;
    variances[i] = objval;
  }

  // Display efficient frontier points
  printf("RETURN    VARIANCE\n");
  for (int i = 0; i < steps; i++)
    printf("%-6f  %-6f\n", returns[i], variances[i]);

  AMPL_Free(&ampl);
  
  return 0;
}