! ch13_polynomial_models.f95
! Companion code for Chapter 13: Optimization and Data Fitting.
!
! Compare polynomial degrees using training and independent validation data.
! Fit by Householder QR with centered, scaled x values. No external libraries.
! Edit data and degrees in the driver. Console output only.
! Smaller training error does not necessarily imply better predictions.
! Reserve fresh test data for an unbiased assessment after model selection.
!
! Compile (Fortran 95):
! gfortran -std=f95 -Wall -Wextra -O2 ch13_polynomial_models.f95 -o ch13_polynomial_models
! Run: ./ch13_polynomial_models (Windows: ch13_polynomial_models.exe)

module polynomial_methods
    implicit none
    private
    integer, parameter :: dp=selected_real_kind(15,307)
    public :: dp, polynomial_fit, polynomial_value, initialize_random, normal_random
contains
    subroutine polynomial_fit(x,y,degree,coeff,mu)
        real(dp), intent(in) :: x(:),y(:)
        integer, intent(in) :: degree
        real(dp), intent(out) :: coeff(:),mu(2)
        real(dp) :: A(size(x),degree+1),rhs(size(x)),v(size(x)),z(size(x))
        real(dp) :: length,scale,diagonal,projection,rank_tol
        integer :: n,m,i,j,k

        n=size(x)
        m=degree+1
        if (degree<0 .or. n<m .or. n<2 .or. size(y)/=n .or. size(coeff)/=m) &
            stop 'Invalid polynomial degree or data dimensions.'
        if (.not. all(abs(x)<=huge(1.0_dp)) .or. .not. all(abs(y)<=huge(1.0_dp))) &
            stop 'Data must be finite.'
        mu(1)=sum(x)/real(n,dp)
        mu(2)=sqrt(sum((x-mu(1))**2)/real(n-1,dp))
        if (mu(2)<=0.0_dp .or. .not. (mu(2)<=huge(1.0_dp))) stop 'Require a usable spread in x.'
        z=(x-mu(1))/mu(2)
        ! Columns are increasing powers of z: 1,z,z^2,... .
        A(:,1)=1.0_dp
        do j=2,m
            A(:,j)=A(:,j-1)*z
        end do
        if (.not. all(abs(A)<=huge(1.0_dp))) stop 'Polynomial basis overflow: lower the degree.'
        rhs=y
        rank_tol=epsilon(1.0_dp)*real(max(n,m),dp)*maxval(abs(A))
        do k=1,m
            ! Construct a normalized Householder reflector for column k.
            v=0.0_dp
            v(k:n)=A(k:n,k)
            scale=maxval(abs(v(k:n)))
            if (scale<=rank_tol) stop 'Polynomial design matrix is numerically rank deficient.'
            length=scale*sqrt(sum((v(k:n)/scale)**2))
            diagonal=-sign(length,v(k))
            v(k)=v(k)-diagonal
            scale=maxval(abs(v(k:n)))
            v(k:n)=v(k:n)/scale
            length=sqrt(sum(v(k:n)**2))
            v(k:n)=v(k:n)/length
            do j=k,m
                projection=2.0_dp*sum(v(k:n)*A(k:n,j))
                A(k:n,j)=A(k:n,j)-projection*v(k:n)
            end do
            projection=2.0_dp*sum(v(k:n)*rhs(k:n))
            rhs(k:n)=rhs(k:n)-projection*v(k:n)
            A(k,k)=diagonal
            if (k<n) A(k+1:n,k)=0.0_dp
        end do
        ! Back substitution in the triangular QR system.
        coeff=0.0_dp
        do i=m,1,-1
            if (abs(A(i,i))<=rank_tol) stop 'Polynomial design matrix is numerically rank deficient.'
            coeff(i)=(rhs(i)-sum(A(i,i+1:m)*coeff(i+1:m)))/A(i,i)
        end do
    end subroutine polynomial_fit

    function polynomial_value(coeff,mu,x) result(y)
        real(dp), intent(in) :: coeff(:),mu(2),x
        real(dp) :: y,z
        integer :: j
        z=(x-mu(1))/mu(2)
        y=coeff(size(coeff))
        do j=size(coeff)-1,1,-1
            y=y*z+coeff(j)
        end do
    end function polynomial_value
    subroutine initialize_random(seed)
        integer, intent(in) :: seed
        integer, allocatable :: seeds(:)
        integer :: n,i
        call random_seed(size=n)
        allocate(seeds(n))
        do i=1,n
            seeds(i)=1+modulo(modulo(seed,1000000)+37*modulo(i,1000000),1000000)
        end do
        call random_seed(put=seeds)
        deallocate(seeds)
        ! Repeatable on the same compiler/runtime; sequences depend on compiler.
    end subroutine initialize_random

    function normal_random() result(z)
        real(dp) :: z,u1,u2
        do
            call random_number(u1)
            if (u1>0.0_dp) exit
        end do
        call random_number(u2)
        z=sqrt(-2.0_dp*log(u1))*cos(2.0_dp*acos(-1.0_dp)*u2)
    end function normal_random

end module polynomial_methods

program ch13_polynomial_models
    use polynomial_methods, only: dp, polynomial_fit, polynomial_value, initialize_random, normal_random
    implicit none
    integer, parameter :: n=18,n_validation=101,n_models=3,seed=7
    integer, parameter :: degrees(n_models)=(/1,5,14/)
    integer, parameter :: max_degree=14 ! Set to at least the largest selected degree.
    real(dp), parameter :: noise_sigma=0.18_dp
    real(dp) :: x(n),y(n),x_validation(n_validation),y_validation(n_validation)
    real(dp) :: coefficients(max_degree+1),mu(2),train_sse,validation_sse,r,pi
    integer :: i,k,m
    ! integer :: clock_seed  ! Uncomment for clock-based seeding.

    if (any(degrees<0) .or. any(degrees>max_degree)) stop 'Adjust max_degree for the selected degrees.'
    call initialize_random(seed)
    ! Alternative: comment out the line above, uncomment the clock_seed
    ! declaration above, and uncomment both lines below.
    ! call system_clock(count=clock_seed)
    ! call initialize_random(clock_seed)

    pi=acos(-1.0_dp)
    do i=1,n
        x(i)=real(i-1,dp)/real(n-1,dp)
        y(i)=sin(2.0_dp*pi*x(i))+0.5_dp*x(i)+noise_sigma*normal_random()
    end do
    ! Independent measurements, never used to determine the coefficients.
    do i=1,n_validation
        x_validation(i)=real(i-1,dp)/real(n_validation-1,dp)
        y_validation(i)=sin(2.0_dp*pi*x_validation(i))+0.5_dp*x_validation(i)+noise_sigma*normal_random()
    end do
    ! Replace both loops with measured training and validation data as needed.
    ! Keep validation data out of the coefficient fit. Fits are unweighted.

    write(*,'(/,a)') 'Polynomial model comparison'
    write(*,'(a)') 'Errors below are RMSE, using N in the denominator.'
    do k=1,n_models
        m=degrees(k)+1
        call polynomial_fit(x,y,degrees(k),coefficients(1:m),mu)
        train_sse=0.0_dp
        do i=1,n
            r=y(i)-polynomial_value(coefficients(1:m),mu,x(i))
            train_sse=train_sse+r*r
        end do
        validation_sse=0.0_dp
        do i=1,n_validation
            r=y_validation(i)-polynomial_value(coefficients(1:m),mu,x_validation(i))
            validation_sse=validation_sse+r*r
        end do
        write(*,'(/,a,i0)') 'Degree:          ',degrees(k)
        write(*,'(a,es16.8)') 'Training RMSE:   ',sqrt(train_sse/real(n,dp))
        write(*,'(a,es16.8)') 'Validation RMSE: ',sqrt(validation_sse/real(n_validation,dp))
        write(*,'(a,2es20.12)') 'Center and scale: ',mu
        write(*,'(a)') 'Coefficients of increasing powers of z = (x-center)/scale:'
        do i=1,m
            write(*,'(a,i0,a,es23.15)') 'c(',i-1,') = ',coefficients(i)
        end do
    end do
    write(*,'(/,a)') 'A single validation set does not establish the best model in general.'
end program ch13_polynomial_models