! 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