## Copyright 2021-2023 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.

using Distributed
@everywhere using FastGaussQuadrature
@everywhere using LinearAlgebra
@everywhere using Polynomials
@everywhere using Printf
@everywhere using SpecialFunctions
include("chebyshev.jl")
include("integration.jl")
include("interpolation.jl")
include("tools.jl")
include("multithread.jl")
include("optimization.jl")
include("potentials.jl")
include("print.jl")
include("easyeq.jl")
include("simpleq-Kv.jl")
include("simpleq-iteration.jl")
include("simpleq-hardcore.jl")
include("anyeq.jl")
function main()
  ## defaults

  # values of rho,e, when needed
  rho=1e-6
  e=1e-4

  # potential
  v=k->v_exp(k,1.)
  a0=a0_exp(1.)
  # arguments of the potential
  v_param_a=1.
  v_param_b=1.
  v_param_c=1.
  v_param_e=1.

  # plot range when plotting in rho
  # linear
  minrho=1e-6
  maxrho=1e2
  nrho=0
  # logarithmic
  minlrho=-6.
  maxlrho=2.
  nlrho=100
  # list
  rhos=Array{Float64,1}(undef,0)
  # plot range when plotting in e
  minle=-6.
  maxle=2.
  nle=100
  es=Array{Float64,1}(undef,0)
  # plot range when plotting in x
  xmin=0.
  xmax=100.
  nx=100
  # plot range when plotting in k
  kmin=0.
  kmax=10.
  nk=100

  # cutoffs
  tolerance=1e-11
  order=100
  maxiter=21

  # for anyeq
  # P
  P=11
  # N
  N=12
  # number of splines
  J=10

  # starting rho from which to incrementally initialize Newton algorithm
  minlrho_init=-6.
  nlrho_init=0

  # Hann window for Fourier transforms
  windowL=1e3

  # approximation for easyeq
  # bK,bL
  easyeq_simpleq_approx=Easyeq_approx(0.,0.)
  easyeq_medeq_approx=Easyeq_approx(1.,1.)
  easyeq_approx=easyeq_simpleq_approx

  # approximation for anyeq
  # aK,bK,gK,aL1,bL1,aL2,bL2,gL2,aL3,bL3,gl3
  anyeq_simpleq_approx=Anyeq_approx(0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.)
  anyeq_bigeq_approx=Anyeq_approx(1.,1.,1.,1.,1.,1.,1.,1.,0.,0.,0.)
  anyeq_fulleq_approx=Anyeq_approx(1.,1.,1.,1.,1.,1.,1.,1.,1.,0.,1.)
  anyeq_medeq_approx=Anyeq_approx(0.,1.,1.,0.,1.,0.,0.,0.,0.,0.,0.)
  anyeq_compleq_approx=Anyeq_approx(1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.)
  anyeq_approx=anyeq_bigeq_approx

  # numerical approximations of derivatives
  dx=1e-7
  dk=1e-7

  # initial guess for 2pt_max
  x0=1.
  k0=1.
  # maximum step in 2pt_max
  maxstep=Inf
  # tolerance for max search
  tolerance_max=Inf

  # read cli arguments
  (params,potential,method,savefile,command)=read_args(ARGS)

  # read params
  if params!=""
    for param in split(params,";")
      terms=split(param,"=")
      if length(terms)!=2
	print(stderr,"error: could not read parameter '",param,"'.\n")
	exit(-1)
      end
      lhs=string(terms[1])
      rhs=string(terms[2])
      if lhs=="rho"
	rho=parse(Float64,rhs)
      elseif lhs=="minlrho_init"
	minlrho_init=parse(Float64,rhs)
      elseif lhs=="nlrho_init"
	nlrho_init=parse(Int64,rhs)
      elseif lhs=="e"
	e=parse(Float64,rhs)
      elseif lhs=="tolerance"
	tolerance=parse(Float64,rhs)
      elseif lhs=="order"
	order=parse(Int64,rhs)
      elseif lhs=="maxiter"
	maxiter=parse(Int64,rhs)
      elseif lhs=="rhos"
	rhos=parse_list(rhs)
      elseif lhs=="minlrho"
	minlrho=parse(Float64,rhs)
      elseif lhs=="maxlrho"
	maxlrho=parse(Float64,rhs)
      elseif lhs=="nlrho"
	nlrho=parse(Int64,rhs)
      elseif lhs=="minrho"
	minrho=parse(Float64,rhs)
      elseif lhs=="maxrho"
	maxrho=parse(Float64,rhs)
      elseif lhs=="nrho"
	nrho=parse(Int64,rhs)
      elseif lhs=="es"
	es=parse_list(rhs)
      elseif lhs=="minle"
	minle=parse(Float64,rhs)
      elseif lhs=="maxle"
	maxle=parse(Float64,rhs)
      elseif lhs=="nle"
	nle=parse(Int64,rhs)
      elseif lhs=="xmin"
	xmin=parse(Float64,rhs)
      elseif lhs=="xmax"
	xmax=parse(Float64,rhs)
      elseif lhs=="nx"
	nx=parse(Int64,rhs)
      elseif lhs=="kmin"
	kmin=parse(Float64,rhs)
      elseif lhs=="kmax"
	kmax=parse(Float64,rhs)
      elseif lhs=="nk"
	nk=parse(Int64,rhs)
      elseif lhs=="P"
	P=parse(Int64,rhs)
      elseif lhs=="N"
	N=parse(Int64,rhs)
      elseif lhs=="J"
	J=parse(Int64,rhs)
      elseif lhs=="window_L"
	windowL=parse(Float64,rhs)
      elseif lhs=="dx"
	dx=parse(Float64,rhs)
      elseif lhs=="x0"
	x0=parse(Float64,rhs)
      elseif lhs=="dk"
	dk=parse(Float64,rhs)
      elseif lhs=="k0"
	k0=parse(Float64,rhs)
      elseif lhs=="maxstep"
	maxstep=parse(Float64,rhs)
      elseif lhs=="tolerance_max"
	tolerance_max=parse(Float64,rhs)
      elseif lhs=="aK"
	anyeq_approx.aK=parse(Float64,rhs)
      elseif lhs=="bK"
	anyeq_approx.bK=parse(Float64,rhs)
	easyeq_approx.bK=parse(Float64,rhs)
      elseif lhs=="gK"
	anyeq_approx.gK=parse(Float64,rhs)
      elseif lhs=="aL1"
	anyeq_approx.aL1=parse(Float64,rhs)
      elseif lhs=="bL"
	easyeq_approx.bL=parse(Float64,rhs)
      elseif lhs=="bL1"
	anyeq_approx.bL1=parse(Float64,rhs)
      elseif lhs=="aL2"
	anyeq_approx.aL2=parse(Float64,rhs)
      elseif lhs=="bL2"
	anyeq_approx.bL2=parse(Float64,rhs)
      elseif lhs=="gL2"
	anyeq_approx.gL2=parse(Float64,rhs)
      elseif lhs=="aL3"
	anyeq_approx.aL3=parse(Float64,rhs)
      elseif lhs=="bL3"
	anyeq_approx.bL3=parse(Float64,rhs)
      elseif lhs=="gL3"
	anyeq_approx.gL3=parse(Float64,rhs)
      elseif lhs=="v_a"
	v_param_a=parse(Float64,rhs)
      elseif lhs=="v_b"
	v_param_b=parse(Float64,rhs)
      elseif lhs=="v_c"
	v_param_c=parse(Float64,rhs)
      elseif lhs=="v_e"
	v_param_e=parse(Float64,rhs)
      elseif lhs=="eq"
	if rhs=="simpleq"
	  easyeq_approx=easyeq_simpleq_approx
	  anyeq_approx=anyeq_simpleq_approx
	elseif rhs=="medeq"
	  easyeq_approx=easyeq_medeq_approx
	  anyeq_approx=anyeq_medeq_approx
	elseif rhs=="bigeq"
	  anyeq_approx=anyeq_bigeq_approx
	elseif rhs=="fulleq"
	  anyeq_approx=anyeq_fulleq_approx
	elseif rhs=="compleq"
	  anyeq_approx=anyeq_compleq_approx
	else
	    print(stderr,"error: unrecognized equation: ",rhs,"\n")
	    exit(-1)
	end
      else
	print(stderr,"unrecognized parameter '",lhs,"'.\n")
	exit(-1)
      end
    end
  end

  ## read potential
  if potential=="exp"
    v=k->v_exp(k,v_param_a)
    a0=a0_exp(v_param_a)
  elseif potential=="expcry"
    v=k->v_expcry(k,v_param_a,v_param_b)
    a0=a0_expcry(v_param_a,v_param_b)
  elseif potential=="npt"
    v=k->v_npt(k)
    a0=a0_npt()
  elseif potential=="alg"
    v=v_alg
    a0=a0_alg
  elseif potential=="algwell"
    v=v_algwell
    a0=a0_algwell
  elseif potential=="exact"
    v=k->v_exact(k,v_param_b,v_param_c,v_param_e)
    a0=a0_exact(v_param_b,v_param_c,v_param_e)
  elseif potential=="tent"
    v=k->v_tent(k,v_param_a,v_param_b)
    a0=a0_tent(v_param_a,v_param_b)
  else
    print(stderr,"unrecognized potential '",potential,"'.\n'")
    exit(-1)
  end

  ## set parameters
  # rhos
  if length(rhos)==0
    # linear only if nrho is specified
    if nrho>0
      rhos=Array{Float64,1}(undef,nrho)
      for j in 0:nrho-1
	rhos[j+1]=(nrho==1 ? minrho : minrho+(maxrho-minrho)/(nrho-1)*j)
      end
    else
      rhos=Array{Float64,1}(undef,nlrho)
      for j in 0:nlrho-1
	rhos[j+1]=(nlrho==1 ? 10^minlrho : 10^(minlrho+(maxlrho-minlrho)/(nlrho-1)*j))
      end
    end
  end

  # es
  if length(es)==0
    es=Array{Float64,1}(undef,nle)
    for j in 0:nle-1
      es[j+1]=(nle==1 ? 10^minle : 10^(minle+(maxle-minle)/(nle-1)*j))
    end
  end

  # splines
  taus=Array{Float64,1}(undef,J+1)
  for j in 0:J
    taus[j+1]=-1+2*j/J
  end

  # tolerance_max
  if tolerance_max==Inf
    tolerance_max=tolerance
  end

  ## run command
  if command=="scattering_length"
    @printf("% .15e\n",a0)
  elseif method=="easyeq"
    if command=="energy"
      easyeq_energy(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx)
    # e(rho)
    elseif command=="energy_rho"
      easyeq_energy_rho(rhos,minlrho_init,nlrho_init,order,a0,v,maxiter,tolerance,easyeq_approx)
    # u(k)
    elseif command=="uk"
      easyeq_uk(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx)
    # u(x)
    elseif command=="ux"
      easyeq_ux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,easyeq_approx)
    elseif command=="uux"
      easyeq_uux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,easyeq_approx)
    # condensate fraction
    elseif command=="condensate_fraction"
      easyeq_condensate_fraction(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx)
    elseif command=="condensate_fraction_rho"
      easyeq_condensate_fraction_rho(rhos,minlrho_init,nlrho_init,order,a0,v,maxiter,tolerance,easyeq_approx)
    # 2pt correlation
    elseif command=="2pt"
      easyeq_2pt(xmin,xmax,nx,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxiter,tolerance,easyeq_approx)
    # max of 2pt correlation
    elseif command=="2pt_max"
      easyeq_2pt_max(dx,x0,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx)
    # max of 2pt correlation as a function of rho
    elseif command=="2pt_max_rho"
      easyeq_2pt_max_rho(rhos,dx,x0,minlrho_init,nlrho_init,order,windowL,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx)
    # fourier transform of 2pt correlation
    elseif command=="2pt_fourier"
      easyeq_2pt_fourier(kmin,kmax,nk,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxiter,tolerance,easyeq_approx)
    # max of fourier transform of 2pt correlation
    elseif command=="2pt_fourier_max"
      easyeq_2pt_fourier_max(dk,k0,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx)
    # max of 2pt correlation as a function of rho
    elseif command=="2pt_fourier_max_rho"
      easyeq_2pt_fourier_max_rho(rhos,dk,k0,minlrho_init,nlrho_init,order,windowL,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx)
    # momentum distribution
    elseif command=="momentum_distribution"
      easyeq_momentum_distribution(kmin,kmax,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxiter,tolerance,easyeq_approx)
    else
      print(stderr,"unrecognized command '",command,"'.\n")
      exit(-1)
    end
  elseif method=="simpleq-hardcore"
    if command=="energy_rho"
      simpleq_hardcore_energy_rho(rhos,taus,P,N,J,maxiter,tolerance)
    elseif command=="ux"
      simpleq_hardcore_ux(rho,taus,P,N,J,maxiter,tolerance)
    elseif command=="condensate_fraction_rho"
      simpleq_hardcore_condensate_fraction_rho(rhos,taus,P,N,J,maxiter,tolerance)
    end
  elseif method=="simpleq-iteration"
    # u_n(x) using iteration
    if command=="u"
      simpleq_iteration_ux(order,e,v,maxiter,xmin,xmax,nx)
    # rho(e) using iteration
    elseif command=="rho_e"
      simpleq_iteration_rho_e(es,order,v,maxiter)
    else
      print(stderr,"unrecognized command '",command,"'.\n")
      exit(-1)
    end
  elseif method=="anyeq"
    if command=="save_Abar"
      anyeq_save_Abar(taus,P,N,J,v,anyeq_approx)
    elseif command=="energy"
      anyeq_energy(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
    # e(rho)
    elseif command=="energy_rho"
      anyeq_energy_rho(rhos,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
    elseif command=="energy_rho_init_prevrho"
      anyeq_energy_rho_init_prevrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
    elseif command=="energy_rho_init_nextrho"
      anyeq_energy_rho_init_nextrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
    # u(j)
    elseif command=="uk"
      anyeq_uk(minlrho_init,nlrho_init,taus,P,N,J,rho,a0,v,maxiter,tolerance,anyeq_approx,savefile)
    # u(j)
    elseif command=="ux"
      anyeq_ux(minlrho_init,nlrho_init,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,anyeq_approx,savefile)
    # condensate fraction
    elseif command=="condensate_fraction"
      anyeq_condensate_fraction(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
    elseif command=="condensate_fraction_rho"
      anyeq_condensate_fraction_rho(rhos,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
    # momentum distribution
    elseif command=="momentum_distribution"
      anyeq_momentum_distribution(kmin,kmax,rho,minlrho_init,nlrho_init,taus,P,N,J,windowL,a0,v,maxiter,tolerance,anyeq_approx,savefile)
    elseif command=="2pt"
      anyeq_2pt_correlation(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,anyeq_approx,savefile)
    elseif command=="2pt_max"
      anyeq_2pt_correlation_max(rho,minlrho_init,nlrho_init,dx,x0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile)
    elseif command=="2pt_max_rho"
      anyeq_2pt_correlation_max_rho(rhos,minlrho_init,nlrho_init,dx,x0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile)
    elseif command=="2pt_fourier"
      anyeq_2pt_correlation_fourier(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,kmin,kmax,nk,anyeq_approx,savefile)
    elseif command=="2pt_fourier_test"
      anyeq_2pt_correlation_fourier_test(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmax,kmin,kmax,nk,anyeq_approx,savefile)
    elseif command=="2pt_fourier_max"
      anyeq_2pt_correlation_fourier_max(rho,minlrho_init,nlrho_init,dk,k0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile)
    elseif command=="2pt_fourier_max_rho"
      anyeq_2pt_correlation_fourier_max_rho(rhos,minlrho_init,nlrho_init,dk,k0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile)
    elseif command=="uncondensed_2pt"
      anyeq_uncondensed_2pt_correlation(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,anyeq_approx,savefile)
    # compressibility
    elseif command=="compressibility_rho"
      anyeq_compressibility_rho(rhos,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
    else
      print(stderr,"unrecognized command: '",command,"'\n")
      exit(-1)
    end
  elseif method=="simpleq-Kv"
    if command=="2pt"
      simpleq_Kv_2pt(minlrho_init,nlrho_init,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx)
    elseif command=="condensate_fraction_rho"
      simpleq_Kv_condensate_fraction(rhos,taus,P,N,J,a0,v,maxiter,tolerance)
    elseif command=="Kv"
      simpleq_Kv_Kv(minlrho_init,nlrho_init,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx)
    else
      print(stderr,"unrecognized command: '",command,"'\n")
      exit(-1)
    end
  else
    print(stderr,"unrecognized method: '",method,"'\n")
    exit(-1)
  end

end

# parse a comma separated list as an array of Float64
function parse_list(
  str::String
)
  elems=split(str,",")
  out=Array{Float64,1}(undef,length(elems))
  for i in 1:length(elems)
    out[i]=parse(Float64,elems[i])
  end
  return out
end

# read cli arguments
function read_args(
  ARGS
)
  # flag
  flag=""

  # output strings
  params=""
  # default potential
  potential="exp"
  # default method
  method="easyeq"
  
  savefile=""
  command=""

  # loop over arguments
  for arg in ARGS
    # argument that starts with a dash
    if arg[1]=='-'
      # go through characters after dash
      for char in arg[2:length(arg)]
	
	# set params
	if char=='p'
	  # raise flag
	  flag="params"
	elseif char=='U'
	  # raise flag
	  flag="potential"
	elseif char=='M'
	  # raise flag
	  flag="method"
	elseif char=='s'
	  # raise flag
	  flag="savefile"
	else
	  print_usage()
	  exit(-1)
	end
      end
    # arguments that do not start with a dash
    else
      if flag=="params"
	params=arg
      elseif flag=="potential"
	potential=arg
      elseif flag=="method"
	method=arg
      elseif flag=="savefile"
	savefile=arg
      else
	command=arg
      end
      # reset flag
      flag=""
    end
  end

  if command==""
    print_usage()
    exit(-1)
  end

  return (params,potential,method,savefile,command)
end

# usage
function print_usage()
  print(stderr,"usage:  simplesolv [-p params] [-U potential] [-M method] [-s savefile] <command>\n")
end

main()