module defgpdata
    implicit none
    real(8),parameter::pi=3.14159265358979323846D0
    !Structure for drawing
    type gpdata
        integer::ndiv=1000
        character::scom*100,sxlabel*100,sylabel*100
        character::slogx*3,slogy*3
        character::sxfmt*20,syfmt*20
        real(4)::xmin,xmax,dx,ymin,ymax,dy
    end type gpdata
    type(gpdata)::gp
end module defgpdata


program f90_SIMPLEX
    !Curve fitting by Simplex method
    use defgpdata
    implicit none
    integer::i,j
    integer::ndata,mpara
    real(8),allocatable::datax(:),datay(:)
    real(8),allocatable::param(:),paraini(:)
    integer::iteration,iflag
    real(8)::work

    character::fnameR*100,fnameW*100,fnameF*100
    character::linebuf*1000,dummy*200

    mpara=mpara_set()    !Number of regression coefficients
    allocate(param(0:mpara-1))
    allocate(paraini(0:mpara-1))

    call getarg(1,fnameR) !Input file name
    call getarg(2,fnameW) !Output file name
    call getarg(3,fnameF) !Figure file name by gnuplot (eps)
    do i=0,mpara-1
        call getarg(4+i,dummy) !Initial values of regrssion coefficients
        read(dummy,*) paraini(i)   !To convert character 'dummy' to real 'paraini'
    end do

    !Data input
    open(11,file=fnameR,status='old')
    read(11,'(a)') gp%scom
    read(11,*) gp%sxlabel,gp%sylabel
    read(11,*) gp%slogx,gp%xmin,gp%xmax,gp%dx,gp%sxfmt
    read(11,*) gp%slogy,gp%ymin,gp%ymax,gp%dy,gp%syfmt
    read(11,*) ndata
    ndata=ndata-1
    allocate(datax(0:ndata))
    allocate(datay(0:ndata))
    do i=0,ndata
        read(11,*),datax(i),datay(i)
    end do
    close(11)

    call SIMPLEX(ndata,datax,datay,mpara,param,paraini,iteration,iflag)

    !Print out
    open(12, file=fnameW, status='replace')
    write(12,'("iflag,",i0)') ,iflag
    write(12,'("iteration,",i0)') ,iteration
    do i=0,mpara-1
        write(linebuf,'("param(",i0,"),",e15.7)') i,param(i)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    write (12,'(a)') 'datax,datay,y_estimate,residual'
    do i=0,ndata
        work=Vfunc(datax(i),param,mpara)
        write(linebuf,*) datax(i),',',datay(i),',',work,',',datay(i)-work
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    close(12)

    !Make file for drawing
    open(12, file='_temp0.prn', status='replace')
    do i=0,ndata
        write(12,*) datax(i),datay(i)
    end do
    close(12)
    open(12, file='_temp1.prn', status='replace')
    do i=0,gp%ndiv
        work=gp%xmin+(gp%xmax-gp%xmin)/dble(gp%ndiv)*dble(i)
        write(12,*) work,Vfunc(work,param,mpara)
    end do
    close(12)

    call GPLOT(fnameF)

    call system('gnuplot "gpl_temp.txt"')

    stop

contains

    integer function mpara_set()
        !*****************************************************
        !Number of parameters
        !*****************************************************
        !Lorentz type+linear curve (mpara=5)   input file: inp_SIMP0.csv
        !mpara_set=5

        !Damped oscillation (mpara=4)   input file: inp_SIMP1.csv
        mpara_set=4
    end function mpara_set


    real(8) function Vfunc(x,param,mpara)
        !**********************************************
        !Calculation of value of function: Vfunc
        !param(i): regression coefficient (i=0 to mpara-1)
        !**********************************************
        integer::mpara
        real(8),intent(in)::x
        real(8),intent(in)::param(0:mpara-1)

        !Lorentz type+linear curve (mpara=5)   input file: inp_SIMP0.csv
        !Vfunc=param(2)*param(1)**2.0D0/((x-param(0))**2.0D0+param(1)**2.0D0)+param(3)+param(4)*x

        !Damped oscillation (mpara=4)   input file: inp_SIMP1.csv
        Vfunc=param(0)*exp(-param(1)*x)*sin(param(2)*x+param(3))
    end function Vfunc


    subroutine INIVAL(mpara,vertex,stepsize,tolerance)
        !*******************************************
        !Set initial values
        !vertex   : initial values of parameters
        !stepsize : initial values of searce steps
        !tolerance: values for judgement of convergence
        !*******************************************
        integer,intent(in)::mpara
        real(8),intent(out)::vertex(0:mpara,0:mpara-1)
        real(8),intent(out)::stepsize(0:mpara-1)
        real(8),intent(out)::tolerance(0:mpara-1)

        integer::j

        do j=0,mpara-1
            vertex(0,j)=1.0D0
            stepsize(j)=0.1D0*vertex(0,j)
            tolerance(j)=1.0D-6
        end do
    end subroutine INIVAL


    subroutine SIMPLEX(ndata,datax,datay,mpara,param,paraini,iteration,iflag)
        integer,intent(in)::ndata,mpara
        integer,intent(out)::iteration,iflag
        real(8),intent(in)::datax(0:ndata),datay(0:ndata),paraini(0:mpara-1)
        real(8),intent(inout)::param(0:mpara-1)

        real(8)::vertex(0:mpara,0:mpara-1)
        real(8)::criterion(0:mpara)
        real(8)::stepsize(0:mpara-1)
        real(8)::tolerance(0:mpara-1)
        real(8)::s,v1,v2,cr,v,vmax,vmin,vsum,p,x,y
        integer::i,j,k,kworst,kbest,ks

        call INIVAL(mpara,vertex,stepsize,tolerance)

        do i=0,mpara-1
            vertex(0,i)=paraini(i)
        end do

        s=sqrt(0.5)
        v1=s*(sqrt(dble(mpara+1))-1.0D0)/dble(mpara)
        do j=0,mpara-1
            v2=vertex(0,j)+v1*stepsize(j)
            do k=1,mpara
                vertex(k,j)=v2
            end do
            vertex(j+1,j)=v2+s*stepsize(j)
        end do
        do k=0,mpara
            do j=0,mpara-1
                param(j)=vertex(k,j)
            end do
            cr=0.0D0    !Vcrit-Start
            do i=0,ndata
                x=datax(i)
                y=Vfunc(x,param,mpara)
                cr=cr+(y-datay(i))*(y-datay(i))
            end do    !Vcrit-End
            criterion(k)=cr
        end do
        iteration=0
        do iteration=1,10000
            kworst=0
            kbest=0
            do k=1,mpara
                if(criterion(k)>criterion(kworst))kworst=k
                if(criterion(k)<criterion(kbest))kbest=k
            end do
            iflag=0
            do j=0,mpara-1
                vmax=vertex(0,j)
                vmin=vmax
                do k=1,mpara
                    v=vertex(k,j)
                    if(v>=vmax)vmax=v
                    if(v<vmin)vmin=v
                end do
                if(vmax-vmin>tolerance(j))iflag=1
            end do
            if(iflag==0)exit
            do j=0,mpara-1
                vsum=0
                do k=0,mpara
                    if(k/=kworst)vsum=vsum+vertex(k,j)
                end do
                param(j)=2.0D0*vsum/dble(mpara)-vertex(kworst,j)
            end do
            cr=0.0D0    !Vcrit-Start
            do i=0,ndata
                x=datax(i)
                y=Vfunc(x,param,mpara)
                cr=cr+(y-datay(i))*(y-datay(i))
            end do    !Vcrit-End
            if(cr<criterion(kbest))ks=0
            if(criterion(kbest)<=cr.and.cr<=criterion(kworst))ks=1
            if(criterion(kworst)<cr)ks=2
            select case(ks)
                case(0)
                    do j=0,mpara-1
                        p=1.5D0*param(j)-0.5D0*vertex(kworst,j)
                        vertex(kworst,j)=param(j)
                        param(j)=p
                    end do
                    criterion(kworst)=cr
                    cr=0.0D0    !Vcrit-Start
                    do i=0,ndata
                        x=datax(i)
                        y=Vfunc(x,param,mpara)
                        cr=cr+(y-datay(i))*(y-datay(i))
                    end do    !Vcrit-End
                    if(cr<criterion(kworst))then
                        do j=0,mpara-1    !REPLACE-Start
                            vertex(kworst,j)=param(j)
                        end do
                        criterion(kworst)=cr    !REPLACE-End
                    end if
                case(1)
                    do j=0,mpara-1    !REPLACE-Start
                        vertex(kworst,j)=param(j)
                    end do
                    criterion(kworst)=cr    !REPLACE-End
                case(2)
                    do j=0,mpara-1
                        param(j)=0.75D0*vertex(kworst,j)+0.25D0*param(j)
                    end do
                    cr=0.0D0    !Vcrit-Start
                    do i=0,ndata
                        x=datax(i)
                        y=Vfunc(x,param,mpara)
                        cr=cr+(y-datay(i))*(y-datay(i))
                    end do    !Vcrit-End
                    if(cr<criterion(kworst))then
                        do j=0,mpara-1    !REPLACE-Start
                            vertex(kworst,j)=param(j)
                        end do
                        criterion(kworst)=cr    !REPLACE-End
                    else
                        do k=0,mpara
                            if(k/=kbest)then
                                do j=0,mpara-1
                                    param(j)=0.5D0*(vertex(k,j)+vertex(kbest,j))
                                    vertex(k,j)=param(j)
                                end do
                                cr=0.0D0    !Vcrit-Start
                                do i=0,ndata
                                    x=datax(i)
                                    y=Vfunc(x,param,mpara)
                                    cr=cr+(y-datay(i))*(y-datay(i))
                                end do    !Vcrit-End
                                criterion(k)=cr
                            end if
                        end do
                    end if
            end select
        end do
        do j=0,mpara-1
            param(j)=vertex(kbest,j)
        end do
    end subroutine SIMPLEX


    subroutine GPLOT(fnameF)
        character::fnameF*100
        character::swork1*50,swork2*50
        character::linebuf*1000
        real(4)::work
        integer::i

        !/****************************/
        !/*gnuplot制御用ファイル作成*/
        !/****************************/
        open(12, file="gpl_temp.txt", status='replace')
        linebuf='reset' ;write(12,'(a)') trim(linebuf)
        linebuf='set terminal png'                    ;write(12,'(a)') trim(linebuf)
        linebuf='set output "'//trim(fnameF)//'"'     ;write(12,'(a)') trim(linebuf)
        linebuf='set xlabel "'//trim(gp%sxlabel)//'"' ;write(12,'(a)') trim(linebuf)
        linebuf='set ylabel "'//trim(gp%sylabel)//'"' ;write(12,'(a)') trim(linebuf)
        if(index(gp%slogx,"L")/=0)then
            linebuf='set logscale x'                  ;write(12,'(a)') trim(linebuf)
            linebuf='set mxtics 10'                   ;write(12,'(a)') trim(linebuf)
        end if
        if(index(gp%slogx,"N")/=0)then
            write(swork1,*) gp%dx
            linebuf='set xtics '//trim(swork1)            ;write(12,'(a)') trim(linebuf)
            linebuf='set format x "'//trim(gp%sxfmt)//'"' ;write(12,'(a)') trim(linebuf)
        end if
        if(index(gp%slogy,"L")/=0)then
            linebuf='set logscale y'                  ;write(12,'(a)') trim(linebuf)
            linebuf='set mytics 10'                   ;write(12,'(a)') trim(linebuf)
        end if
        if(index(gp%slogx,"N")/=0)then
            write(swork1,*) gp%dy
            linebuf='set ytics '//trim(swork1)            ;write(12,'(a)') trim(linebuf)
            linebuf='set format y "'//trim(gp%sxfmt)//'"' ;write(12,'(a)') trim(linebuf)
        end if
        write(swork1,*) gp%xmin
        write(swork2,*) gp%xmax
        linebuf='set xrange ['//trim(swork1)//':'//trim(swork2)//']' ;write(12,'(a)') trim(linebuf)
        write(swork1,*) gp%ymin
        write(swork2,*) gp%ymax
        linebuf='set yrange ['//trim(swork1)//':'//trim(swork2)//']' ;write(12,'(a)') trim(linebuf)
        linebuf='set grid xtics ytics mxtics mytics'     ;write(12,'(a)') trim(linebuf)
        linebuf='plot \'                                 ;write(12,'(a)') trim(linebuf)
        linebuf='"_temp0.prn" with points pointtype 1 pointsize 1.5 title "Input data", \' ;write(12,'(a)') trim(linebuf)
        linebuf='"_temp1.prn" with lines linetype 3 linewidth 2 title "Estimated"'       ;write(12,'(a)') trim(linebuf)
        linebuf='quit' ;write(12,'(a)') trim(linebuf)
        close(12)
        !/****************************/
        !/*gnuplot制御*/
        !/****************************/
        call system('gnuplot gpl_temp.txt')

    end subroutine GPLOT


    subroutine del_spaces(s)
        character (*), intent (inout) :: s
        character (len=len(s)) tmp
        integer i, j
        j = 1
        do i = 1, len(s)
            if (s(i:i)==' ') cycle
            tmp(j:j) = s(i:i)
            j = j + 1
        end do
        s = tmp(1:j-1)
    end subroutine del_spaces


end program f90_SIMPLEX