/*
Copyright 2016 Ian Jauslin

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

#include <stdio.h>
#include <stdarg.h>
// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears
#define MPFR_USE_VA_LIST
#include <mpfr.h>
#include <math.h>
#include <libinum.h>

#include "types.h"
#include "hh_integral.h"
#include "hh_root.h"
#include "hh_integral_double.h"
#include "hh_root_double.h"
#include "ss_integral_double.h"
#include "zz_integral.h"
#include "zz_integral_double.h"
#include "parser.h"
#include "definitions.h"

// usage message
int print_usage();
// read arguments
int read_args(int argc, char* argv[], hh_params* params, mpfr_t tolerance, unsigned int* maxiter, unsigned int* order, unsigned int* computation_nr, unsigned int* threads, unsigned int* use_double);

// compute the first loop correction to the topological phase diagram
int compute_hh(hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order);
// using doubles instead of mpfr
int compute_hh_double(hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order);

// compute the sunrise diagram
int compute_ss(TYPE_I, hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order, unsigned int threads);
// using doubles instead of mpfr
int compute_ss_double(TYPE_I_DOUBLE, hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order, unsigned int threads);

// codes for possible computations
#define COMPUTATION_PHASE 1
#define COMPUTATION_ZZ 2
#define COMPUTATION_ZZZZ 3

int main(int argc, char* argv[]){
  hh_params params;
  mpfr_t tolerance;
  unsigned int maxiter;
  unsigned int order;
  int ret;
  unsigned int computation_nr;
  unsigned int threads;
  unsigned int use_double;

  // default computation: phase diagram
  computation_nr=COMPUTATION_PHASE;

  mpfr_inits(params.W, params.sinphi, params.phi, params.t1, params.t2, params.lambda, tolerance, NULL);
  // read command line arguments
  ret=read_args(argc, argv, &params, tolerance, &maxiter, &order, &computation_nr, &threads, &use_double);
  if(ret<0){
    mpfr_clears(params.W, params.sinphi, params.phi, params.t1, params.t2, params.lambda, tolerance, NULL);
    return(-1);
  }
  if(ret>0){
    mpfr_clears(params.W, params.sinphi, params.phi, params.t1, params.t2, params.lambda, tolerance, NULL);
    return(0);
  }

  // phase diagram
  if(computation_nr==COMPUTATION_PHASE){
    // compute the first-loop correction to the phase diagram
    if(use_double==0){
      compute_hh(params, tolerance, maxiter, order);
    }
    else{
      compute_hh_double(params, tolerance, maxiter, order);
    }
  }
  else if(computation_nr==COMPUTATION_ZZ){
    // compute second-order correction to z1-z2
    if(use_double==0){
      compute_ss(&zz_I, params, tolerance, maxiter, order, threads);
    }
    else{
      compute_ss_double(&zz_I_double, params, tolerance, maxiter, order, threads);
    }
  }
  else if(computation_nr==COMPUTATION_ZZZZ){
    // compute second-order correction to z1+z2
    if(use_double==0){
      compute_ss(&ZZ_I, params, tolerance, maxiter, order, threads);
    }
    else{
      compute_ss_double(&ZZ_I_double, params, tolerance, maxiter, order, threads);
    }
  }

  mpfr_clears(params.W, params.sinphi, params.phi, params.t1, params.t2, params.lambda, tolerance, NULL);

  return(0);
}

// usage message
int print_usage(){
  fprintf(stderr, "usage:\n       hhtop phase [-p params] [-v] [-O order] [-t tolerance] [-N maxiter] [-P precision] [-E emax]\n       hhtop z1-z2 [-p params] [-v] [-O order] [-t tolerance] [-N maxiter] [-P precision] [-E emax] [-C threads]\n       hhtop z1+z2 [-p params] [-v] [-O order] [-t tolerance] [-N maxiter] [-P precision] [-E emax] [-C threads]\n\n       hhtop -D phase [-p params] [-v] [-O order] [-t tolerance] [-N maxiter]\n       hhtop -D z1-z2 [-p params] [-v] [-O order] [-t tolerance] [-N maxiter] [-C threads]\n\n       hhtop -V [-v]\n\n");
  return(0);
}

// read command line arguments
#define CP_FLAG_PARAMS 1
#define CP_FLAG_ORDER 2
#define CP_FLAG_TOLERANCE 3
#define CP_FLAG_MAXITER 4
#define CP_FLAG_MPFR_PREC 5
#define CP_FLAG_MPFR_EXP 6
#define CP_FLAG_THREADS 7
int read_args(int argc, char* argv[], hh_params* params, mpfr_t tolerance, unsigned int* maxiter, unsigned int* order, unsigned int* computation_nr, unsigned int* threads, unsigned int* use_double){
  int i;
  int ret;
  // temporary long int
  long int tmp_lint;
  // temporary unsigned int
  unsigned int tmp_uint;
  // pointers
  char* ptr;
  // flag that indicates what argument is being read
  int flag=0;
  // pointer to various arguments
  char* tolerance_str=NULL;
  char* params_str=NULL;
  // whether to print the variables after they are read
  int print_vars=0;
  // keep track of which flags were used (to check for incompatibilities)
  unsigned char pflag=0;
  unsigned char Oflag=0;
  unsigned char tflag=0;
  unsigned char Nflag=0;
  unsigned char Pflag=0;
  unsigned char Eflag=0;
  unsigned char vflag=0;
  unsigned char Cflag=0;
  unsigned char Dflag=0;
  unsigned char Vflag=0;

  // defaults
  mpfr_set_d(params->t1, 1., MPFR_RNDN);
  mpfr_set_d(params->t2, .1, MPFR_RNDN);
  mpfr_set_d(params->lambda, .01, MPFR_RNDN);
  mpfr_set_d(params->sinphi, 1., MPFR_RNDN);
  mpfr_set_d(tolerance, 1e-11, MPFR_RNDN);
  *maxiter=1000000;
  *order=10;
  params->omega=+1;
  *threads=1;
  *use_double=0;

  mpfr_const_pi(params->phi, MPFR_RNDN);
  mpfr_div_ui(params->phi, params->phi, 2, MPFR_RNDN);

  // default W=-3*sqrt(3)*t2*sin(phi)
  mpfr_sqrt_ui(params->W, 3, MPFR_RNDN);
  mpfr_mul_si(params->W, params->W, -3, MPFR_RNDN);
  mpfr_mul(params->W, params->W, params->sinphi, MPFR_RNDN);
  mpfr_mul(params->W, params->W, params->t2, MPFR_RNDN);

  // loop over arguments
  for(i=1;i<argc;i++){
    // flag
    if(argv[i][0]=='-'){
      for(ptr=((char*)argv[i])+1;*ptr!='\0';ptr++){
	switch(*ptr){
	// parameters
	case 'p':
	  flag=CP_FLAG_PARAMS;
	  pflag=1;
	  break;
	// order of the integration
	case 'O':
	  flag=CP_FLAG_ORDER;
	  Oflag=1;
	  break;
	// tolerance
	case 't':
	  flag=CP_FLAG_TOLERANCE;
	  tflag=1;
	  break;
	// maximal number of Newton steps
	case 'N':
	  flag=CP_FLAG_MAXITER;
	  Nflag=1;
	  break;
	// mpfr precision
	case 'P':
	  flag=CP_FLAG_MPFR_PREC;
	  Pflag=1;
	  break;
	// mpfr emax
	case 'E':
	  flag=CP_FLAG_MPFR_EXP;
	  Eflag=1;
	  break;
	// print value of variables
	case 'v':
	  print_vars=1;
	  vflag=1;
	  break;
	// number of threads
	case 'C':
	  flag=CP_FLAG_THREADS;
	  Cflag=1;
	  break;
	// use doubles instead of mpfr
	case 'D':
	  *use_double=1;
	  Dflag=1;
	  // set prec to that of long doubles (for consistency when reading cli arguments with many digits)
	  mpfr_set_default_prec(64);
	  break;
	// print version
	case 'V':
	  Vflag=1;
	  break;
	default:
	  fprintf(stderr, "unrecognized option '-%c'\n", *ptr);
	  print_usage();
	  return(-1);
	  break;
	}
      }
    }
    // parameters
    else if(flag==CP_FLAG_PARAMS){
      // read str later (after having set the MPFR precision and emax)
      params_str=argv[i];
      flag=0;
    }
    // order of the integration
    else if(flag==CP_FLAG_ORDER){
      ret=sscanf(argv[i],"%u",&tmp_uint);
      if(ret!=1){
	fprintf(stderr, "error: '-O' should be followed by an unsigned int\n       got '%s'\n",argv[i]);
	return(-1);
      }
      *order=tmp_uint;
      flag=0;
    }
    // tolerance
    else if(flag==CP_FLAG_TOLERANCE){
      // read str later (after having set the MPFR precision and emax)
      tolerance_str=argv[i];
      flag=0;
    }
    // maximal number of Newton steps
    else if(flag==CP_FLAG_MAXITER){
      ret=sscanf(argv[i],"%u",maxiter);
      if(ret!=1){
	fprintf(stderr, "error: '-N' should be followed by a positive int\n       got '%s'\n",argv[i]);
	return(-1);
      }
      flag=0;
    }
    // mpfr precision
    else if(flag==CP_FLAG_MPFR_PREC){
      ret=sscanf(argv[i],"%ld",&tmp_lint);
      if(ret!=1){
	fprintf(stderr, "error: '-P' should be followed by a long int\n       got '%s'\n",argv[i]);
	return(-1);
      }
      mpfr_set_default_prec(tmp_lint);
      flag=0;
    }
    // mpfr emax
    else if(flag==CP_FLAG_MPFR_EXP){
      ret=sscanf(argv[i],"%ld",&tmp_lint);
      if(ret!=1){
	fprintf(stderr, "error: '-E' should be followed by a long int\n       got '%s'\n",argv[i]);
	return(-1);
      }
      mpfr_set_emax(tmp_lint);
      flag=0;
    }
    // number of threads to use for the computation
    else if(flag==CP_FLAG_THREADS){
      ret=sscanf(argv[i],"%u",threads);
      if(ret!=1){
	fprintf(stderr, "error: '-C' should be followed by a positive int\n       got '%s'\n",argv[i]);
	return(-1);
      }
      flag=0;
    }
    // computation to run
    else{
      if(str_cmp(argv[i], "phase")==1){
	*computation_nr=COMPUTATION_PHASE;
      }
      else if(str_cmp(argv[i], "z1-z2")==1){
	*computation_nr=COMPUTATION_ZZ;
      }
      else if(str_cmp(argv[i], "z1+z2")==1){
	*computation_nr=COMPUTATION_ZZZZ;
      }
      else{
	fprintf(stderr, "error: unrecognized computation: '%s'\n",argv[i]);
	print_usage();
	return(-1);
      }
      flag=0;
    }
  }
  if(tolerance_str!=NULL){
    ret=mpfr_set_str(tolerance, tolerance_str, 10, MPFR_RNDN);
    if(ret<0){
      fprintf(stderr, "error: '-t' should be followed by an MPFR floating point number\n       got '%s'\n", tolerance_str);
      return(-1);
    }
  }
  if(params_str!=NULL){
    ret=read_params(params, params_str);
    if(ret<0){
      return(ret);
    }
  }

  // check for incompatible flags
  if((Vflag==1 && (pflag!=0 || Oflag!=0 || tflag!=0 || Nflag!=0 || Pflag!=0 || Eflag!=0 || Cflag!=0 || Dflag!=0)) || \
     (*computation_nr!=COMPUTATION_ZZ && *computation_nr!=COMPUTATION_ZZZZ && Cflag==1) || \
     (Dflag==1 && (Pflag==1 || Eflag==1)) \
     ){
    print_usage();
    return(-1);
  }

  // print version and exit
  if(Vflag==1){
    printf("hhtop " VERSION "\n");
    printf("libinum " LIBINUM_VERSION "\n");
    if(vflag==1){
      // print datatype information
      printf("\n\n");
      print_datatype_info(stdout);
    }
    return(1);
  }

  // print variables
  if(print_vars==1){
    fprintf(stderr, "t1=");
    fprint_mpfr(stderr,params->t1);
    fprintf(stderr, "\n");

    fprintf(stderr, "t2=");
    fprint_mpfr(stderr,params->t2);
    fprintf(stderr, "\n");

    fprintf(stderr, "lambda=");
    fprint_mpfr(stderr,params->lambda);
    fprintf(stderr, "\n");

    fprintf(stderr, "phi=");
    fprint_mpfr(stderr,params->phi);
    fprintf(stderr, "\n");

    fprintf(stderr, "sinphi=");
    fprint_mpfr(stderr,params->sinphi);
    fprintf(stderr, "\n");

    fprintf(stderr, "W=");
    fprint_mpfr(stderr,params->W);
    fprintf(stderr, "\n");

    fprintf(stderr, "omega=%d\n",params->omega);

    fprintf(stderr, "\ntolerance=");
    fprint_mpfr(stderr,tolerance);
    fprintf(stderr, "\n");

    fprintf(stderr, "order of integration=%d\n", *order);

    fprintf(stderr, "\nMPFR precision=%ld\n", mpfr_get_default_prec());
    fprintf(stderr, "MPFR emax=%ld\n", mpfr_get_emax());

    fprintf(stderr, "\n");
  }

  return(0);
}

// compute the first loop correcion to the topological phase diagram
int compute_hh(hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order){
  mpfr_t val;
  int ret;
  args_integration args_int;

  // compute weights
  ret=gauss_legendre_weights_mpfr(order, tolerance, maxiter, &(args_int.abcissa), &(args_int.weights));
  // return codes
  if(ret==LIBINUM_ERROR_MAXITER){
    fprintf(stderr, "error: maximum number of iterations reached when computing the integration abcissa\n       try increasing the precision or the tolerance\n");
    return(ret);
  }
  else if(ret==LIBINUM_ERROR_NAN){
    fprintf(stderr, "error: infinity encountered when computing the integration abcissa\n");
    return(ret);
  }

  // set args
  args_int.params=params;
  mpfr_init(args_int.params.W);

  mpfr_init(val);
  // initial value
  mpfr_sqrt_ui(val, 3, MPFR_RNDN);
  mpfr_mul_ui(val, val, 3, MPFR_RNDN);
  mpfr_mul(val, val, params.sinphi, MPFR_RNDN);
  mpfr_mul(val, val, params.t2, MPFR_RNDN);
  if(params.omega==1){
    mpfr_neg(val, val, MPFR_RNDN);
  }

  // compute root
  ret=root_newton_inplace_mpfr(&val, &integration_wrapper, &d_integration_wrapper, tolerance, maxiter, &args_int);
  // return codes
  if(ret==LIBINUM_ERROR_MAXITER){
    fprintf(stderr, "error: maximum number of iterations reached when computing the solution to m=0\n       try increasing the precision or the tolerance\n");
    mpfr_clear(val);
    mpfr_clear(args_int.params.W);
    array_mpfr_free(args_int.abcissa);
    array_mpfr_free(args_int.weights);
    return(ret);
  }
  else if(ret==LIBINUM_ERROR_NAN){
    fprintf(stderr, "error: infinity encountered: either the integrand is singular or the derivative of the integral vanishes at a point of the Newton iteration\n");
    mpfr_clear(val);
    mpfr_clear(args_int.params.W);
    array_mpfr_free(args_int.abcissa);
    array_mpfr_free(args_int.weights);
    return(ret);
  }

  fprint_mpfr(stdout, val);
  printf("\n");

  mpfr_clear(val);
  mpfr_clear(args_int.params.W);
  array_mpfr_free(args_int.abcissa);
  array_mpfr_free(args_int.weights);

  return(0);
}
// using double instead of mpfr
int compute_hh_double(hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order){
  long double tolerance_d;
  hh_args_integration_double args_int;
  long double val;
  int ret;

  // convert mpfr to double
  tolerance_d=mpfr_get_ld(tolerance, MPFR_RNDN);
  hh_params_todouble(&(args_int.params), params);

  // compute weights
  ret=gauss_legendre_weights_ldouble(order, tolerance_d, maxiter, &(args_int.abcissa), &(args_int.weights));

  // return codes
  if(ret==LIBINUM_ERROR_MAXITER){
    fprintf(stderr, "error: maximum number of iterations reached when computing the integration abcissa\n");
    return(ret);
  }
  else if(ret==LIBINUM_ERROR_NAN){
    fprintf(stderr, "error: infinity encountered when computing the integration abcissa\n");
    return(ret);
  }


  // initial value
  val=-args_int.params.omega*3*sqrtl(3)*args_int.params.t2*args_int.params.sinphi;

  ret=root_newton_inplace_ldouble(&val, &hh_integration_wrapper_double, &hh_d_integration_wrapper_double, tolerance_d, maxiter, &args_int);

  // return codes
  if(ret==LIBINUM_ERROR_MAXITER){
    fprintf(stderr, "error: maximum number of iterations reached when computing the solution to m=0\n");
    array_ldouble_free(args_int.abcissa);
    array_ldouble_free(args_int.weights);
    return(ret);
  }
  else if(ret==LIBINUM_ERROR_NAN){
    fprintf(stderr, "error: infinity encountered: either the integrand is singular or the derivative of the integral vanishes at a point of the Newton iteration\n");
    array_ldouble_free(args_int.abcissa);
    array_ldouble_free(args_int.weights);
    return(ret);
  }

  printf("% .19Le\n",val);

  array_ldouble_free(args_int.abcissa);
  array_ldouble_free(args_int.weights);

  return(0);
}

// compute the sunrise diagram
int compute_ss(TYPE_I, hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order, unsigned int threads){
  mpfr_t val;
  int ret;
  array_mpfr abcissa;
  array_mpfr weights;

  // compute weights
  ret=gauss_legendre_weights_mpfr(order, tolerance, maxiter, &abcissa, &weights);
  // return codes
  if(ret==LIBINUM_ERROR_MAXITER){
    fprintf(stderr, "error: maximum number of iterations reached when computing the integration abcissa\n       try increasing the precision, the tolerance, or the maximal number of Newton steps\n");
    return(ret);
  }
  else if(ret==LIBINUM_ERROR_NAN){
    fprintf(stderr, "error: infinity encountered when computing the integration abcissa\n");
    return(ret);
  }

  mpfr_init(val);

  // compute integral
  ret=ss_integrate(&val, I, params, abcissa, weights, threads);
  // return codes
  if(ret==LIBINUM_ERROR_NAN){
    fprintf(stderr, "error: infinity encountered: the integrand is singular\n");
    mpfr_clear(val);
    array_mpfr_free(abcissa);
    array_mpfr_free(weights);
    return(ret);
  }

  fprint_mpfr(stdout, val);
  printf("\n");

  mpfr_clear(val);
  array_mpfr_free(abcissa);
  array_mpfr_free(weights);

  return(0);
}
// using doubles instead of mpfr
int compute_ss_double(TYPE_I_DOUBLE, hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order, unsigned int threads){
  long double tolerance_d;
  hh_params_double params_d;
  long double val;
  int ret;
  array_ldouble abcissa;
  array_ldouble weights;

  // convert mpfr to double
  tolerance_d=mpfr_get_ld(tolerance, MPFR_RNDN);
  hh_params_todouble(&params_d, params);

  // compute weights
  ret=gauss_legendre_weights_ldouble(order, tolerance_d, maxiter, &abcissa, &weights);
  // return codes
  if(ret==LIBINUM_ERROR_MAXITER){
    fprintf(stderr, "error: maximum number of iterations reached when computing the integration abcissa\n       try increasing the tolerance, or the maximal number of Newton steps\n");
    return(ret);
  }
  else if(ret==LIBINUM_ERROR_NAN){
    fprintf(stderr, "error: infinity encountered when computing the integration abcissa\n");
    return(ret);
  }

  // compute integral
  ret=ss_integrate_double(&val, I, params_d, abcissa, weights, threads);
  // return codes
  if(ret==LIBINUM_ERROR_NAN){
    fprintf(stderr, "error: infinity encountered: the integrand is singular\n");
    array_ldouble_free(abcissa);
    array_ldouble_free(weights);
    return(ret);
  }

  printf("% .19Le\n",val);

  array_ldouble_free(abcissa);
  array_ldouble_free(weights);

  return(0);
}