""" Interpos Python package
Antoine Merle and Olivier Sauter
https://crppwww.epfl.ch/~sauter/interpos/

This package uses f2py to allow access in python to subroutines of the interpos fortran library.
For now it is not possible to use the fortran generic interface.
A python function has been used to handle the different types of calls.
This function processes the arguments and calls the corresponding fortran subroutines 
(possible choices are intlinear, intquadratic, interpos_def or interpos_defper).
In the future, the function could be modified to allow for optional arguments
(but remember that in python, you cannot use a positional argument after a keyword argument).

The interpos python package requires the installation of numpy.

"""

__author__ = "Antoine Merle (antoine.merle@epfl.ch)"
__version__ = "8.3"

__all__ = ['interpos']


import interpos_fun
import numpy as np

def interpos(kopt,xin,yin,xout,taus,nbc,ybc,sigmain):

    """Interpolation routine based on a python wrapper of the fortran library interpos.

    = documentation based on the one of the matlab mex file =

    Typical use:

    yout,youtp,youtpp,youtint = interpos(kopt,xin,yin,xout,taus,nbc,ybc,sigmain)

    For now, all the arguments are mandatory.
    There should be exactly 1 or 4 output arguments.
    If there is only 1 output argument, then it will be a tuple of size 4 with each index corresponding respectively to yout, youtp, youtpp and youtint.
 
    The arguments have the following meaning with respect to the function y(x)
    that one wants to interpolate/extrapolate or compute its derivatives and integrals:
 
    kopt: option for interpolation and extrapolation method (default is kopt=13):
    .     = 1, 11 or 21 : linear interpolation
    .     = 2, 12 or 22 : quadratic "
    .     = 3, 13 or 23 : cubic spline interpolation, with tension=taus
    .     < 10 : send warning if need to extrapolate
    .     in [11,19] : extrapolate using a lower order than for interpolation
    .     in [21,29] : extrapolate using same order as interpolation
 
    xin : array giving the input values of x
    yin : array giving the input values of y(x=xin(i))
 
    xout: array of values of x at which the function and/or its derivatives
    .     are to be computed. If xout is not given, assumes xout=xin
    yout: interpolated values of y at x=xout(i).
    youtp: 1st derivative of y at x=xout(i)
    youtpp: 2nd derivative of y at x=xout(i)
    youtint: Integral of (y dx) from x=xin(1) to x=xout(i)
 
    taus  : tension value for cubic spline interpolation. If not given, uses taus=0
            if taus < 0 (typically -1) uses a default taus value: taus=abs(taus)*default_taus
            (the default_taus is typically min(delta_x)**3)
 
    nbc(2): [NBCLFT NBCRGT] (default: [0 0])
      BOUNDARY CONDITIONS, 4 TYPES DETERMINED BY THE VALUE OF (nbc(1)=NBCLFT
      and nbc(2) for left and right-hand side BC.):
 
      0) VALUE OF SECOND DERIVATIVE AT XBCLFT OR RGT IS GIVEN (0 OR 10)
      1) VALUE OF 1ST        "       "   "     "  "   "   "   (1 OR 11)
      2) VALUE OF FUNCTION AT XBCLFT OR RGT IS GIVEN          (2 OR 12)
      if nbc(1)=-1, then assumes periodic boundary condition and uses ybc(1) for the period
 
      The value of nbc(1 or 2) should be >= 10 if the BC is not at an end point:
      XBCLFT~=xin(1) or XBCRGT~=xin(end)
 
      Examples:
        A good value for radial profiles in rho between [0,1] is to specify
        the first derivative = 0 at left and second derivative=0 at right:
        => nbc = [1 0] and ybc=[0. 0.]
 
        or to obtained the first derivative at right-hand side from a
        lagrangian interpolation of the last points:
        => nbc=[1 1] and ybc=[0. 1e32]
 
        Note: The B.C. can only be given anywhere but within the interval [xin(1),xin(end)]
        Therefore, say rho=[0.1 ... 1] and one wants to impose zero
        derivative at rho=0, one should add a point in the input before. One uses sigmain to avoid forcing the value
        calling interpos:
              rho_eff(1)=0.;
              rho_eff(2:length(rho)+1) = rho;
        then  [yout]=interpos([0.;rho],[yin(1);yin],..,taus,[1 0],[0. 0.],[1000;ones(size(rho))]);
 
    ybc(2 or 4 or 6) (default: [0 0]): [YBCLFT YBCRGT] or [YBCLFT YBCRGT XBCLFT XBCRGT] with
      THE VALUE IS GIVEN BY YBCLFT OR YBCRGT RESPECTIVELY.
 
      FOR nbc type 1: IF (YBCLFT OR YBCRGT > 1E31 THEN DER. FROM LAGRANGIAN INTERP.
      FOR nbc type 1: IF (YBCLFT OR YBCRGT <-1E31 THEN DER. FROM LINEAR     INTERP.
 
      IF NBCLFT OR NBCRGT IS < 10, PXIN(1) OR PXIN(KNIN) IS USED INSTEAD
      OF XBCLFT OR XBCRGT, RESPECTIVELY => XBCLFT OR XBCRGT NOT USED
 
    ybc(6): [YBCLFT YBCRGT XBCLFT XBCRGT PXEXP0 PXEXPDL]
         Enables one to specify a gaussian wight to taus with PXEXP0 and PXEXPDL:
               taus_eff =  taus * EXP(-((XIN-PXEXP0)/PXEXPDL)**2)
         IF PXEXP0 NOT IN [PXIN(1),PXIN(KNIN)], EXP() IGNORED AND taus_eff=cst= taus
         if ybc(5:6) not given, then  PXEXP0=xin(1)-1. and pxexpdl=1. to get constant taus
 
    ybc(1)=period for periodic boundary condition, if nbc(1)=-1. Defines the condition that y(x+period)=y(x)
           In this case xin(end) should not be equal to xin(1)+period, since it is redundant. However interpos just
           does not use the end point in this case.
         if period=-1, assumes period=xin(end)-xin(1). Useful if xin goes from 0 to 2pi for example
 
     sigmain: error_bar at each (x,y) point used for the fit. The effective taus value will then
            be taus(i)=taus .* sigmain(i) ./ min(sigmain). So uses the relative values of sigmain.
    """

    # adapted from interpos_libs/interpos_matlab_top.f90

    def fun_sigma(xx): ztaueff*exp(-zcofexp*(xx-zxexp0)**2/zxexpdl**2)
    #
    #.......................................................................
    #
    #   1. defaults values
    #
    #
    #   2. Input arguments
    #
    #
    #   2.1 kopt
    #
    zkopt = kopt
    kopt = np.fabs(kopt)
    kopt_sign = np.sign(zkopt)
    #
    #   define interpolation type and extrapolation type
    #
    nel = 10
    inttype = np.mod(kopt,10)
    iextrapo = kopt/nel
    #
    #   2.2 get xin
    #
    nin = np.size(xin)
    if ((inttype==3) and (nin <= 3)):
        if (nin <= 2):
            inttype = 1
            print 'nin = ',nin,' <3 , cannot compute spline, uses linear interpolation'
        else:
            inttype = 2
            print 'nin = ',nin,' =3 , cannot compute spline, uses quadratic interpolation'
    #
    #   2.3 get yin
    #
    #
    #   2.4 get xout
    #
    nout = np.size(xout)
    #
    #   2.5 get taus
    #
    #
    #   2.6 get nbc,ybc
    #
    znbc = nbc
    zybc = ybc
    nbc = np.array([0,0])
    ybc = np.array([0.,0.,0.,0.,0.,0.])
    idoexp_error=0
    if np.size(znbc) > 0:
        if np.size(znbc) == 1:
            nbc[0] = znbc[0]
            nbc[1] = znbc[0]
            ybc[0] = -1.
        else:
            nbc[0] = znbc[0]
            nbc[1] = znbc[1]
    if np.size(zybc) == 0:
        ybc[0] = 0.
        ybc[1] = 0.
        idoexp_error = 0
    else:
        ybc[:np.size(zybc)] = zybc
        if np.size(zybc) <= 4:
            ybc[4] = xin[0]-xin[-1]
            ybc[5] = 1.
        else:
            idoexp_error = 1
    #
    #   2.6 get simain
    #
    if sigmain is None:
        ztaueff = 1.
        zxexp0 = ybc[4]
        zxexpdl = ybc[5]
        zcofexp = 0.
        if idoexp_error == 1:
            zcofexp = 1.
        sigmain = fun_simain(xin)
    #
    #   3. Allocate return arguments arrays (in same geometry as xin)
    #
    yout = np.ones(nout)
    youtp = np.ones(nout)
    youtpp = np.ones(nout)
    youtint = np.ones(nout)
    iflag = 0
    #   pointers for yout, youtp and youtpp and youtint respectively
    # ioptder = 0
    # ioptder = 1
    # ioptder = 2
    # ioptder = 3
    if (nin == 1):
        yout = yin
        youtp = 0.
        youtpp = 0.
        youtint = 0.
        return yout,youtp,youtpp,youtint
    #
    #   4. compute interpolated function and its derivatives
    #    
    if (inttype == 1):
        if (iextrapo == 0):
            iextrapo = 0
        elif (iextrapo == 1):
            iextrapo = 1
        elif (iextrapo == 2):
            iextrapo = 10
        else:
            iextrapo = 1
        #
        iextrapo = kopt_sign * iextrapo
        interpos_fun.intlinear(xin,yin,nin,xout,yout,youtp,youtpp,youtint,nout,ioptder,iextrapo)
    #
    if (inttype == 2):
        if (iextrapo == 0):
            iextrapo = 0
        elif (iextrapo == 1):
            iextrapo = 21
        elif (iextrapo == 2):
            iextrapo = 1
        elif (iextrapo == 3):
            iextrapo = 2
        elif (iextrapo == 4):
            iextrapo = 10
        else:
            iextrapo = 2
        #
        iextrapo = kopt_sign * iextrapo
        interpos_fun.intquadratic(xin,yin,nin,xout,yout,youtp,youtpp,youtint,nout,ioptder,iextrapo,nbc[0])
    #
    if (inttype == 3):
        if (iextrapo == 0):
            iextrapo = 0
        elif (iextrapo == 1):
            iextrapo = 32
        elif (iextrapo == 2):
            iextrapo = 1
        elif (iextrapo == 3):
            iextrapo = 2
        elif (iextrapo == 4):
            iextrapo = 3
        elif (iextrapo == 5):
            iextrapo = 31
        elif (iextrapo == 6):
            iextrapo = 10
        else:
            iextrapo = 32
        #
        if (nbc[0] < 0):
            # periodic boundary conditions
            interpos_fun.interpos_module.interpos_defper(xin=xin,yin=yin,nin=nin,nout=nout,tension=taus,xout=xout,yout=yout,youtp=youtp,youtpp=youtpp,youtint=youtint,nbc=nbc[0],ybc=ybc[0],sigma=sigmain,info=iflag)
        else:
            # Standard case
            interpos_fun.interpos_module.interpos_def(xin=xin,yin=yin,nin=nin,nout=nout,tension=taus,xout=xout,yout=yout,youtp=youtp,youtpp=youtpp,youtint=youtint,nbc=nbc,ybc=ybc,sigma=sigmain,option=zkopt,info=iflag)
        #
        if (iflag != 0):
            print ' '
            print ' problem in interpos'
            print ' iflag = ',iflag
            return
    
    return yout,youtp,youtpp,youtint