! Program hyperplane_main: driver for hyperplane-newton routine
!
! Author: Robert K. Moniot, Fordham University
! Date:   October 2005
!
! Performs Deming least-squares fit of input data to a straight line in
! 2 dimensions, plane in 3 dimensions, or in general a hyperplane of
! n-1 dimensions in n dimensions.
!
! Input (stdin) consists of 3 lines:
!   First line has 5 numbers separated by space and/or comma:
!     max_iters (integer) = limit on number of iterations (default 10)
!     conv_crit (real) = convergence criterion, precision of result (default 1e-12)
!     v (integer) = verbosity level, 0 for just the facts, 1 for details (default 0)
!     wait_iters (integer) = limit on number of fixed-point iterations to
!             use before starting Newton iteration.  (default 2)
!     wait_crit (real) = convergence criterion
!   Second line is name of method used (newton), not used except for information.
!   Third line is name of file containing data for the fit.
!
! Data file (specified on line 3 of stdin) consists of the following:
!   First line has mdim = number of dimensions of data space (integer),
!      followed by any text to be used as a title.
!   Second line has a series of mdim text items to be used as column
!      headings for the list of data points.
!   Following lines contain the data points.  Each data point starts on
!      a new line, and may span any number of lines.  The data values are
!      given first, then the sigmas of the data values.  Values are separated
!      by blanks, newlines, and/or commas.
!
! Copyright (c) 2005 by Robert K. Moniot.
!
! Permission is hereby granted, free of charge, to any person
! obtaining a copy of this software and associated documentation
! files (the "Software"), to deal in the Software without
! restriction, including without limitation the rights to use,
! copy, modify, merge, publish, distribute, sublicense, and/or
! sell copies of the Software, and to permit persons to whom the
! Software is furnished to do so, subject to the following
! conditions:

! The above copyright notice and this permission notice shall be
! included in all copies or substantial portions of the
! Software.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY
! KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
! WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
! PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
! COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
! OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
! SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
!
! Acknowledgement: the above permission notice is known as the X
! Consortium license.

program hyperplane_main

  use precision_module
  use hyperplane_fit

  implicit none

  real(kind=DP), allocatable, dimension(:,:) :: Y, sigma_Y ! data, +/-
  real(kind=DP), allocatable, dimension(:) :: a, sigma_a  ! params of fit, +/-
  real(kind=DP), allocatable, dimension(:,:) :: cov_a  ! covariance matrix of a
  real(kind=DP) :: chisq
  real(kind=DP), allocatable, dimension(:) :: num_sigma_a ! numerical result
  real(kind=DP), allocatable, dimension(:,:) :: num_cov_a ! numerical result
  character(len=8), allocatable, dimension(:) :: Ylabels ! coord identifiers
  real(kind=DP), allocatable, dimension(:) :: Ypoint, sigma_Ypoint ! temp for input
  real(kind=DP), allocatable, dimension(:,:) :: Ytemp, sigma_Ytemp ! temp for input
  integer :: i, j, npts, mdims
  integer :: stat
  character(len=*), parameter :: datafmt="(a2,(6(1x,f12.5)))"
  character(len=80) :: title

! optional parameters, read from  stdin.
  real(kind=DP) :: conv_crit, wait_crit
  integer :: max_iters, wait_iters, v
  character(len=10) :: method

  character (len=80) :: datfilename
  integer, parameter :: datunit = 11

! Set defaults for input quantities
  max_iters = 10
  conv_crit = 1.0e-14
  v = 0
  wait_iters = 2
  wait_crit = 1e-4
!
!   Input section
!
  read(unit=*,fmt=*) max_iters, conv_crit, v, wait_iters, wait_crit
  read(unit=*,fmt="(a10)") method
  read(unit=*,fmt="(a80)") datfilename

  write(unit=*,fmt=*) "hyperplane program"
  write(unit=*,fmt=*) "max_iters=", max_iters,  &
       " conv_crit=", real(conv_crit,kind=SP), &
       " verbosity=", v
  write(unit=*,fmt=*) "Method=", trim(method)
     write(unit=*,fmt=*) "wait_iters=", wait_iters, &
          " wait_crit=", real(wait_crit,kind=SP)
  write(unit=*,fmt=*) "Reading data from ", datfilename


  open(unit=datunit,file=datfilename,status="old",action="read",iostat=stat)
  if( stat /= 0 ) then
     write(unit=*,fmt=*) "Error opening ", trim(datfilename), &
          " for reading: iostat=", stat
     stop
  end if

  read(unit=datunit,fmt="(a80)") title
  ! read the number of dimensions from 1st part of title
  read(unit=title,fmt=*) mdims
  if( mdims < 2 ) then
     print *, "Error: number of dimensions =",mdims," must be at least 2"
     stop
  end if

  allocate(Ylabels(1:mdims))
  read(unit=datunit,fmt=*) Ylabels

! Ypoint and sigma_Ypoint hold one input data vector & its error vector
! while growing data array to hold it.
  allocate(Ypoint(1:mdims), sigma_Ypoint(1:mdims))

  ! read the data points up to EOF
  ! The arrays Y and sigma_Y grow to hold each new point, so no max on input
  npts = 0
  do
     read(unit=datunit,fmt=*,iostat=stat) Ypoint(1:mdims), sigma_Ypoint(1:mdims)
     if( stat < 0 ) then
        exit
     end if
     if( npts > 0 ) then        ! not the first point
        ! wish we had a REALLOCATE statement
        allocate(Ytemp(1:npts,1:mdims),sigma_Ytemp(1:npts,1:mdims))
        Ytemp = Y               ! hold the previously read data
        sigma_Ytemp = sigma_Y
        deallocate(Y, sigma_Y)
     end if
     npts = npts+1
     allocate(Y(1:npts,1:mdims), sigma_Y(1:npts,1:mdims))
     if( npts > 1 ) then        ! not the first point
        Y(1:npts-1,1:mdims) = Ytemp
        sigma_Y(1:npts-1,1:mdims) = sigma_Ytemp
        deallocate(Ytemp, sigma_Ytemp)
     end if
     Y(npts,1:mdims) = Ypoint
     sigma_Y(npts,1:mdims) = sigma_Ypoint
  end do

  close(datunit)
  write(unit=*,fmt=*) npts, " data points read"
  write(unit=*,fmt="(1x,a85)")  "dims=" // title
  write(unit=*,fmt="(5x,(10a9))") Ylabels
  do i=1, npts
     write(unit=*, fmt=*) "  ", real(Y(i,1:mdims),SP)
     write(unit=*, fmt=*) "+-", real(sigma_Y(i,1:mdims),SP)
  end do

  allocate(a(1:mdims), sigma_a(1:mdims))
  allocate(cov_a(1:mdims,1:mdims))

  call hyperplane(npts, mdims, Y, sigma_Y, a, chisq, stat, V=cov_a, &
       max_iters=max_iters, conv_crit=conv_crit, &
       fixedpt_conv_crit=wait_crit, fixedpt_iters=wait_iters, verbosity=v)

  if( stat /= 0 ) then
     if( stat > 1 ) then
        print *, "Hyperplane fit failed"
        stop
     else
        print *, "Warning: solution did not converge"
     end if
  end if

  if( v >= 1 ) then
     write(unit=*,fmt=*) "cov(a) ="
     do j=1,mdims
        write(unit=*,fmt=*) real(cov_a(j,1:mdims),SP)
     end do
  end if

  do j=1, mdims
          ! Guard against exception in case of bugs
     if( cov_a(j,j) < 0 ) then
        print *, "ERROR: negative variance of a(",j,")=",cov_a(j,j)
        sigma_a(j) = -1.0
     else
        sigma_a(j) = sqrt(cov_a(j,j))
     end if
  end do


  write(unit=*,fmt=*)
  write(unit=*, fmt=*) " a=", real(a(1:mdims),SP)
  write(unit=*, fmt=*) " +-", real(sigma_a(1:mdims),SP)

  write(unit=*, fmt="(a,f10.3)") " chi squared=", chisq


end program hyperplane_main