// qmisim.c   this program executes a  QMI  Quantum Machine Instruction
//            Given the h and J values, print the s values at the minimum of
//            Ising Hamiltonian E(s) = {sum(h_i*s_i) + sum(J_ik*s_i*s_k)}
//            compile   gcc -O qmisim qmisim.c -lm
//            run       qmisim xxx.qmi
//
// xxx.qmi  file format: lines starting with  #  are comments
// number_of_h  number_of_J
// h_index  h_value               # number_of_h  of these lines
// s_i_index  s_k_index  J_value  # number_of_J  of these lines
//
// constraints for this version: -1 <= h <= +1, -1 <= J <= +1, s = -1 or +1
// 0 <= i <= 511, 0 <= k <= 511, maximum of 2048 J_ik
// warnings if input would violate Dwave hardware

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>

#undef max
#define max(a,b) ((a)>(b)?(a):(b))

static int nexts(int ns, double h[], double s[]);
static double f(int ns, int nJ, double h[], double J[],
                int Ji[], int Jk[], double s[]);
static int legalc(int s, int legal[]);
static int checkik(int i, int k);
void randomize(int ns, double h[], double s[], double sr[]);

int main(int argc, char* argv[])
{
  int ns = 0;
  int is = 0;
  double s[512];     // start with 0, max 511
  double sr[512];    // randomized values
  double sbest[512];
  double smin = -1.01;
  double smax =  1.01;
  double y, ebest = 1.0e+10;
  int nh = 0;
  double h[512];
  double hmin = -1.01; // just for checking input
  double hmax =  1.01;
  int nJ = 0;
  int iJ = 0;
  int Ji[2048];
  int Jk[2048];
  double J[2048];
  double Jmin = -1.01; // just for checking input
  double Jmax =  1.01;
  int i, j, k, farg;
  double weight;
  FILE * inp;
  char line[128];
  int done = 0;
  int err = 0;
  int ifrand = 0;
  int debug = 0;

  printf("qmisim running \n");
  if(argc<2) {printf("need file name  xxx.qmi \n"); exit(1);}
  farg = 1;
  if(argc>2) // check for option, just 1 allowed
  {
    farg++;
    if(!strcmp(argv[1],"-r")) ifrand = 1;
    if(!strcmp(argv[1],"-d")) debug = 1;
  }
  if(argc<1+farg) {printf("need file name  xxx.qmi \n"); exit(1);}
  printf("reading QMI from %s \n", argv[farg]);
 
  inp = fopen(argv[farg], "r");
  if(inp == NULL) {printf("can not open for read %s \n", argv[farg]); exit(1);}

  while(feof(inp)==0)
  {
    fgets(line, 128, inp);
    printf("%s", line);
    if(line[0] == '#') continue;
    if(nh==0)
    {
      sscanf(line, "%d %d", &nh, &nJ);
      for(i=0; i<512; i++)
      {
        s[i] = 0.0;
        h[i] = 0.0;
        sbest[i] = 0.0;
      }
    }
    else if(is<nh) // count h lines
    {
      sscanf(line, "%d %lf", &i, &weight);
      h[i] = weight;
      s[i] = -1.0;
      is++;
      ns = max(ns,i); // allow gaps, yet they waste time
      if(debug) printf("h[%d]=%f \n", i, h[i]);
      if(weight<hmin || weight>hmax) printf("warning h out of range\n");
    }
    else  // read Jik
    {
      sscanf(line, "%d %d %lf", &i, &k, &weight);
      Ji[iJ] = i;
      Jk[iJ] = k;
      J[iJ] = weight;
      if(debug) printf("J %d %d = %f \n", Ji[iJ], Jk[iJ], J[iJ]);
      iJ++;
      if(weight<Jmin || weight>Jmax) printf("warning J out of range\n");
      err = checkik(i, k);
    }
    if(iJ==nJ) break;
  }
  printf("qmi read, nh=%d, nJ=%d, largest s,h index=%d\n", nh, nJ, ns);

  done = 0;
  while(done==0) // main computation
  {
    if(ifrand)
    { 
      randomize(ns, h, s, sr);
      y = f(ns, nJ, h, J, Ji, Jk, sr); // evaluate Ising Hamiltonian
    }
    else // non random case
    {
      y = f(ns, nJ, h, J, Ji, Jk, s); // evaluate Ising Hamiltonian
    }
    if(debug>0)
    {
      printf("energy=%g \n", y); // debug print
      for(i=1; i<=ns; i++)
      {
        if(h[i]==0.0) continue;
        if(ifrand==0) printf("s[%d] = %f \n", i, s[i]);
        if(ifrand==1) printf("sr[%d] = %f \n", i, sr[i]);
      }
    }

    if(y<ebest)
    {
      ebest = y;
      for(i=0; i<=ns; i++)
      {
        if(h[i]==0.0) continue;
        if(ifrand==0) sbest[i] = s[i];
        if(ifrand==1) sbest[i] = sr[i];
      }
    }
    done = nexts(ns, h, s); // skip if h==0.0
  } // end while, all combinations evaluated

  printf("\nminimum energy = %g at: \n", ebest);
  for(i=1; i<=ns; i++)
  {
    if(h[i]==0.0) continue;
    if(ifrand==0) printf("s[%d] = %f \n", i, sbest[i]);
    if(ifrand==1) printf("sr[%d] = %f \n", i, sbest[i]);
  }
  printf("\nqmisim finished \n");
  return 0;
} // end main

static int nexts(int ns, double h[], double s[])
{
  int i, j, k;

  for(i=0; i<=ns; i++)
  {
    if(h[i]==0.0) continue; // skip unused 
    j = i;
    if(s[i]!=1.0) break;    
    if(j==ns && s[j]==1.0) return 1;
    s[i] = -1.0;
  }
  // printf("j=%d, s[j]=%f \n", j, s[j]);
  s[j] = 1.0;
  return 0;
} // end nexts

static double f(int ns, int nJ, double h[], double J[],
                int Ji[], int Jk[], double s[])
{
  double e = 0.0;
  int i, j, k;

  for(i=0; i<=ns; i++) e += h[i]*s[i]; // some h may be zero
  for(i=0; i<nJ; i++) e += J[i]*s[Ji[i]]*s[Jk[i]];
  return e;
} // end f 

static int legalc(int s, int legal[])
{
  int localrow;
  int localcol;
  int s4, s8, s16, s64;
  int base;
  int n = 6;

  // local
  s4 = s%4;
  s8 = s%8;
  base = s - s8;
  s16 = s%16;
  s64 = s%64;
  // printf("s=%3d, s4=%3d, s8=%3d, s16=%3d, s64=%3d\n", s, s4, s8, s16, s64);
  localrow = s4;
  localcol = 1;
  if(s8>=4) localcol = 2;
  if(localcol==1)
  {
    legal[0] = base+4;
    legal[1] = base+5;
    legal[2] = base+6;
    legal[3] = base+7;
    legal[4] = base-64+s4;
    if(legal[4]<0) {legal[4] = legal[3]; n--;}
    legal[5] = base+64+s4;
    if(legal[5]>511) {legal[5] = legal[3]; n--;}
  }
  else
  {
    legal[0] = base+0;
    legal[1] = base+1;
    legal[2] = base+2;
    legal[3] = base+3;
    legal[4] = base-8+s8;
    if(s64<=7) {legal[4] = legal[3]; n--;}
    legal[5] = base+8+s8;
    if(s64>=56) {legal[5] = legal[3]; n--;}
  }  
  return n;
} // end legalc

static int checkik(int i, int k)
{
  int legal[6];
  int n;
  int j;

  n = legalc(i, legal);
  for(j=0; j<6; j++)
  {
    if(k==legal[j]) return 0;
  }
  printf("coupling %d %d not valid\n", i, k);
  return 1;
} // end checkik

void randomize(int ns, double h[], double s[], double sr[])
{
  static int first = 1;
  static double two34;
  long int ran;
  int i;

  if(first)
  {
    srandom(time(NULL)); // delete if same sequence wanted
    first = 0;
    two34 = pow(2.0, 34.0); // range 0 to 1/8 for random() 
  }
  for(i=1; i<=ns; i++)
  {
    if(h[i]==0.0) continue;
    sr[i] = s[i] + (random()/two34 - 0.0625); // add range -1/16 to +1/16
  }
} // end randomize

// end qmisim.c