Hi, guys,
I have implemented the least squares using orthogonal multinomials basis. The result is ok in 1D problems, but not good in 2D.
The structure of the codes is:
a module file - orthogonal_multinomials_least_squares.f90
a test main file - test_OMLS_1d.f90 & test_OMLS_2d.f90
Can someone check the codes and find out the problems?
Thanks!![[bigears]](/data/assets/smilies/bigears.gif "[bigears] [bigears]")
I have implemented the least squares using orthogonal multinomials basis. The result is ok in 1D problems, but not good in 2D.
The structure of the codes is:
a module file - orthogonal_multinomials_least_squares.f90
a test main file - test_OMLS_1d.f90 & test_OMLS_2d.f90
Code:
program main
use orthogonal_multinomials_least_squares
implicit none
integer, parameter :: n = 10
real, parameter :: dx = 0.1
real :: x(1,n), f(n), f0(n)
integer :: i, j, k
real :: r
do i = 1, n
x(1,i) = i*dx
call random_number(r)
x(1,i) = x(1,i)+r*dx/10.0
end do
! f = exp(-100.0*(x(1,:)-0.5)**2)
f = (x(1,:)-0.5)**2
call OMLS_init(1, 3)
call OMLS_fit(n, x, f, n, x, f0)
call OMLS_finalize
open(10, file="input_1d.dat", form="unformatted", &
access="direct", recl=sizeof(r)/4*n)
write(10, rec=1) x(1,:)
write(10, rec=2) f
close(10)
open(11, file="output_1d.dat", form="unformatted", &
access="direct", recl=sizeof(r)/4*n)
write(11, rec=1) x(1,:)
write(11, rec=2) f0
close(11)
stop
end program main
Code:
program main
use orthogonal_multinomials_least_squares
implicit none
integer, parameter :: n = 100, nx = 10, ny = 10
real, parameter :: dx = 0.1, dy = 0.1
real :: x(2,n), f(n), f0(n)
integer :: i, j, k
real :: r
do j = 1, ny
do i = 1, nx
k = (j-1)*nx+i
x(1,k) = i*dx
x(2,k) = j*dy
call random_number(r)
x(1,k) = x(1,k)+r*dx/10.0
x(2,k) = x(2,k)+r*dy/10.0
end do
end do
! f = exp(-100.0*((x(1,:)-0.5)**2+(x(2,:)-0.5)**2))
f = (x(1,:)-0.5)**2+(x(2,:)-0.5)**2
call OMLS_init(2, 2)
call OMLS_fit(n, x, f, n, x, f0)
call OMLS_finalize
open(10, file="input_2d.dat", form="unformatted", access="direct", recl=sizeof(r)/4*n)
write(10, rec=1) x(1,:)
write(10, rec=2) x(2,:)
write(10, rec=3) f
close(10)
open(11, file="output_2d.dat", form="unformatted", access="direct", recl=sizeof(r)/4*n)
write(11, rec=1) x(1,:)
write(11, rec=2) x(2,:)
write(11, rec=3) f0
close(11)
stop
end program main
Code:
! orthogonal_multinomials_least_squares module
! Purpose:
! (1) construct orthogonal multinomials according to discrete data
! (2) use the generated orthogonal basis to reconstruct the continuous
! distribution of the data
! Author:
! Name: DONG Li
! Email: dongli@lasg.iap.ac.cn
! Mathematics supporter:
! Name: SONG Xiang
! Email: ?
module orthogonal_multinomials_least_squares
implicit none
integer :: num_dim
integer :: max_exponent
integer :: num_row
integer :: num_monomial
integer, allocatable :: ranges(:,:)
integer, allocatable :: exponents(:,:)
real, allocatable :: coef(:)
real, allocatable :: psi_buf1(:,:)
real, allocatable :: psi_buf2(:)
interface psi
module procedure psi_internal_point
module procedure psi_external_point
end interface psi
contains
subroutine OMLS_init(nd, me)
integer, intent(in) :: nd, me
integer :: i, j, k, m, n
num_dim = nd
max_exponent = me
num_row = me+1
! generate the structure of monomials
! number of monomials in each range
allocate(ranges(2:num_row,num_dim))
ranges(2,:) = 1
do i = 3, num_row
do j = 1, num_dim
ranges(i,j) = sum(ranges(i-1,j:num_dim))
end do
end do
! total number of monomials
num_monomial = sum(ranges)+1
allocate(exponents(num_monomial,num_dim))
exponents(1,:) = 0
do i = 2, num_dim+1
exponents(i,:) = 0
exponents(i,i-1) = 1
end do
m = num_dim+2 ! offset for exponents
do i = 3, num_row
n = 1 ! offset base for exponents in previous row
do k = 2, i-2
n = n+sum(ranges(k,:))
end do
do j = 1, num_dim ! loop for ranges in current row
do k = 1, sum(ranges(i-1,j:num_dim)) ! loop for monomial
exponents(m,:) = exponents(n+k,:)
exponents(m,j) = exponents(m,j)+1
m = m+1
end do
n = n+sum(ranges(i-1,j:j))
end do
end do
allocate(coef(num_monomial))
write(*, "('ranges:')")
do i = 2, num_row
write(*, "('row ', i2, ': ')", advance="no") i
write(*, *)ranges(i,:)
end do
write(*, *)
write(*, "('exponents:')")
do i = 1, num_monomial
write(*, *) exponents(i,:)
end do
return
end subroutine OMLS_init
subroutine OMLS_fit(num_point, x, f, num_point0, x0, f0)
integer, intent(in) :: num_point, num_point0
real, intent(in) :: x(num_dim,num_point), f(num_point)
real, intent(in) :: x0(num_dim,num_point0)
real, intent(out) :: f0(num_point0)
real :: diag_A(num_monomial), B(num_monomial)
real :: tmp
real :: A(num_monomial,num_monomial) ! for debug
integer :: k ! for debug
integer :: i, j
allocate(psi_buf1(num_point,num_monomial))
allocate(psi_buf2(num_monomial))
psi_buf2 = 0.0
! set up the buffers
do j = 1, num_monomial
do i = 1, num_point
psi_buf1(i,j) = psi(num_point, x, j, i)
end do
do i = 1, num_point
psi_buf2(j) = psi_buf2(j)+psi_buf1(i,j)**2
end do
end do
B = 0.0
do j = 1, num_monomial
diag_A(j) = psi_buf2(j)
do i = 1, num_point
B(j) = B(j)+psi_buf1(i,j)*f(i)
end do
end do
coef = B/diag_A
f0 = 0.0
do j = 1, num_point0
do i = 1, num_monomial
f0(j) = f0(j)+coef(i)*psi(num_point, x, i, x0(:,j))
end do
end do
! for debug
A = 0.0
do k = 1, num_monomial
do j = 1, num_monomial
do i = 1, num_point
A(j,k) = A(j,k)+psi_buf1(i,j)*psi_buf1(i,k)
end do
end do
write(*,*) A(:,k)
end do
write(*, "('The coefficients are:')")
do i = 1, num_monomial
write(*, *) coef(i)
end do
return
end subroutine OMLS_fit
real function psi_internal_point(num_point, x, mono_indx, point_indx) result(ans)
integer, intent(in) :: num_point
real, intent(in), target :: x(num_dim,num_point)
integer, intent(in) :: mono_indx
integer, intent(in) :: point_indx
real :: tmp
real, pointer :: xk(:)
integer :: i, j
if(mono_indx == 1) then
ans = 1
return
end if
do i = num_dim, 1, -1
if(exponents(mono_indx,i) /= 0) then
xk => x(i,:)
ans = x(i,point_indx)*psi_buf1(point_indx,mono_indx-1)
exit
end if
end do
do i = 1, mono_indx-1
tmp = 0.0
do j = 1, num_point
tmp = tmp+xk(j)*psi_buf1(j,mono_indx-1)*psi_buf1(j,i)
end do
ans = ans-tmp/psi_buf2(i)*psi_buf1(point_indx,i)
end do
return
end function psi_internal_point
recursive real function psi_external_point(num_point, x, mono_indx, x0) result(ans)
integer, intent(in) :: num_point
real, intent(in), target :: x(num_dim,num_point)
integer, intent(in) :: mono_indx ! current monomial index (recursive)
real, intent(in) :: x0(num_dim)
real :: tmp
real, pointer :: xk(:)
integer :: i, j
if(mono_indx == 1) then
ans = 1
return
end if
do i = num_dim, 1, -1
if(exponents(mono_indx,i) /= 0) then
xk => x(i,:)
ans = x0(i)*psi_external_point(num_point, x, mono_indx-1, x0)
exit
end if
end do
do i = 1, mono_indx-1
tmp = 0.0
do j = 1, num_point
tmp = tmp+xk(j)*psi_buf1(j,mono_indx-1)*psi_buf1(j,i)
end do
ans = ans-tmp/psi_buf2(i)*psi_external_point(num_point, x, i, x0)
end do
return
end function psi_external_point
subroutine OMLS_finalize
if(allocated(ranges)) deallocate(ranges)
if(allocated(exponents)) deallocate(exponents)
if(allocated(coef)) deallocate(coef)
if(allocated(psi_buf1)) deallocate(psi_buf1)
if(allocated(psi_buf2)) deallocate(psi_buf2)
end subroutine OMLS_finalize
end module orthogonal_multinomials_least_squaresCan someone check the codes and find out the problems?
Thanks!
![[dazed]](/data/assets/smilies/dazed.gif "[dazed] [dazed]")
You can plot the result of 2D-case. When the data points are sampled from parabolic surface, the approximated surface should fit the original one very well (I suppose so). Since the basis function is multinomial, the information contained in this approximation system should be enough to capture parabolic surface. The 1D-case is fine.