From ea6b295776efaac4bd4666c87206ad056214b0f1 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Wed, 4 Feb 2026 12:43:08 +0530 Subject: [PATCH 01/29] kabsch first commit --- example/CMakeLists.txt | 1 + example/spatial/CMakeLists.txt | 1 + example/spatial/example_kabsch.f90 | 70 +++++++++++ src/CMakeLists.txt | 2 + src/spatial/CMakeLists.txt | 14 +++ src/spatial/stdlib_spatial.fypp | 30 +++++ src/spatial/stdlib_spatial_kabsch.fypp | 164 +++++++++++++++++++++++++ 7 files changed, 282 insertions(+) create mode 100644 example/spatial/CMakeLists.txt create mode 100644 example/spatial/example_kabsch.f90 create mode 100644 src/spatial/CMakeLists.txt create mode 100644 src/spatial/stdlib_spatial.fypp create mode 100644 src/spatial/stdlib_spatial_kabsch.fypp diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index 8b330c6e5..038d0796a 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -45,6 +45,7 @@ if (STDLIB_QUADRATURE) endif() add_subdirectory(selection) add_subdirectory(sorting) +add_subdirectory(spatial) add_subdirectory(specialfunctions_gamma) if (STDLIB_SPECIALMATRICES) add_subdirectory(specialmatrices) diff --git a/example/spatial/CMakeLists.txt b/example/spatial/CMakeLists.txt new file mode 100644 index 000000000..54540eb25 --- /dev/null +++ b/example/spatial/CMakeLists.txt @@ -0,0 +1 @@ +ADD_EXAMPLE(kabsch) \ No newline at end of file diff --git a/example/spatial/example_kabsch.f90 b/example/spatial/example_kabsch.f90 new file mode 100644 index 000000000..de364b170 --- /dev/null +++ b/example/spatial/example_kabsch.f90 @@ -0,0 +1,70 @@ +program example_kabsch + use stdlib_linalg_constants, only: dp + use stdlib_spatial + implicit none + + integer, parameter :: d = 3, N = 4 + real(dp) :: P(d, N), Q(d, N), Q_1(d, N) + real(dp) :: R(d, d), t(d) + real(dp) :: c, rmsd + logical :: scale + + integer :: i + + ! Reference point set P (columns are points) + P(:,1) = [0.0_dp, 0.0_dp, 0.0_dp] + P(:,2) = [1.0_dp, 0.0_dp, 0.0_dp] + P(:,3) = [0.0_dp, 1.0_dp, 0.0_dp] + P(:,4) = [0.0_dp, 0.0_dp, 1.0_dp] + + ! Target point set Q + ! Here Q is a rotated + translated + scaled version of P + ! Target point set Q = 2 * Rz * P + [1, 2, 3] + +Q(:,1) = [1.0_dp, 2.0_dp, 3.0_dp] ! P1 = (0,0,0) +Q(:,2) = [1.0_dp, 4.0_dp, 3.0_dp] ! P2 = (1,0,0) -> (0,1,0) +Q(:,3) = [-1.0_dp, 2.0_dp, 3.0_dp] ! P3 = (0,1,0) -> (-1,0,0) +Q(:,4) = [1.0_dp, 2.0_dp, 5.0_dp] ! P4 = (0,0,1) + + + scale = .true. + + call kabsch(P, Q, scale, R, t, c, rmsd) + + print *, "" + print *, "Matrix P:" + do i = 1, d + print "(4F10.5)", P(i,:) + end do + + print *, "" + print *, "Matrix Q:" + do i = 1, d + print "(4F10.5)", Q(i,:) + end do + + print *, "Rotation matrix R:" + do i = 1, d + print "(3F10.5)", R(i,:) + end do + + print *, "" + print *, "Translation vector t:" + print "(3F10.5)", t + + ! Apply the transformation: Q_1 = c * R * Q + t + do i = 1, N + Q_1(:,i) = c * matmul(R, Q(:,i)) + t + end do + + print *, "" + print *, "Aligned Q (should match P):" + do i = 1, d + print "(4F10.5)", Q_1(i,:) + end do + + print *, "" + print *, "Scale factor c:", c + print *, "RMSD:", rmsd + +end program example_kabsch diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f955ed352..a4f6dde9e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -44,6 +44,7 @@ ADD_SUBDIR(system) ADD_SUBDIR(stats) add_subdirectory(sparse) +add_subdirectory(spatial) set(fppFiles stdlib_version.fypp @@ -70,5 +71,6 @@ target_link_libraries(${PROJECT_NAME} PUBLIC ${PROJECT_NAME}_strings ${PROJECT_NAME}_blas ${PROJECT_NAME}_lapack ${PROJECT_NAME}_lapack_extended ${PROJECT_NAME}_sparse + ${PROJECT_NAME}_spatial ${OPTIONAL_LIB} ) diff --git a/src/spatial/CMakeLists.txt b/src/spatial/CMakeLists.txt new file mode 100644 index 000000000..c78edd60e --- /dev/null +++ b/src/spatial/CMakeLists.txt @@ -0,0 +1,14 @@ +set(spatial_fppFiles + stdlib_spatial.fypp + stdlib_spatial_kabsch.fypp + ) + +set(spatial_cppFiles + ) + +set(spatial_f90Files + ) + +configure_stdlib_target(${PROJECT_NAME}_spatial spatial_f90Files spatial_fppFiles spatial_cppFiles) + +target_link_libraries(${PROJECT_NAME}_spatial PUBLIC ${PROJECT_NAME}_constants ${PROJECT_NAME}_linalg_core ${PROJECT_NAME}_linalg ${PROJECT_NAME}_intrinsics) \ No newline at end of file diff --git a/src/spatial/stdlib_spatial.fypp b/src/spatial/stdlib_spatial.fypp new file mode 100644 index 000000000..eb72f7f56 --- /dev/null +++ b/src/spatial/stdlib_spatial.fypp @@ -0,0 +1,30 @@ +#:include "common.fypp" +#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) +#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) +#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES +module stdlib_spatial + use stdlib_linalg_constants + use stdlib_constants + use stdlib_linalg, only: svd, det, eye, trace + use stdlib_intrinsics, only: stdlib_sum + implicit none + private + public :: kabsch + + interface kabsch + #:for k, t, s in (KINDS_TYPES) + #:if k!="xdp" + module subroutine kabsch_${s}$(P, Q, scale, R, t, c, rmsd, W) + ${t}$, intent(in) :: P(:, :) + ${t}$, intent(in) :: Q(:, :) + logical, intent(in) :: scale + ${t}$, intent(out) :: R(:,:) + ${t}$, intent(out) :: t(:) + real(${k}$), intent(out) :: c + real(${k}$), intent(out) :: rmsd + real(${k}$), intent(in), optional :: W(:) + end subroutine + #:endif + #:endfor + end interface +end module stdlib_spatial diff --git a/src/spatial/stdlib_spatial_kabsch.fypp b/src/spatial/stdlib_spatial_kabsch.fypp new file mode 100644 index 000000000..1e2a40c73 --- /dev/null +++ b/src/spatial/stdlib_spatial_kabsch.fypp @@ -0,0 +1,164 @@ +#:include "common.fypp" +#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) +#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) +#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES +submodule(stdlib_spatial) stdlib_spatial_kabsch + +contains + #:for k, t, s in (KINDS_TYPES) + #:if k!="xdp" + module subroutine kabsch_${s}$(P, Q, scale, R, t, c, rmsd, W) + ${t}$, intent(in) :: P(:, :) + ${t}$, intent(in) :: Q(:, :) + logical, intent(in) :: scale + ${t}$, intent(out) :: R(:,:) + ${t}$, intent(out) :: t(:) + real(${k}$), intent(out) :: c + real(${k}$), intent(out) :: rmsd + real(${k}$), intent(in), optional :: W(:) + + ! Internal variables. + integer(ilp) :: i, j, k, d, N + ${t}$, allocatable :: c_P(:), c_Q(:), covariance(:,:), U(:,:), Vt(:,:), B(:,:), vec(:) + real(${k}$) :: sum_w, vp, vq, variance_p + real(${k}$), allocatable :: S(:) + #:if t.startswith('complex') + ${t}$ :: detR + #:endif + + if(size(P,dim=1)/=size(Q,dim=1) .or. size(P,dim=1)/=size(R,dim=1) .or. size(P,dim=1)/=size(R,dim=2) & + .or. size(P,dim=1)/=size(t)) then + error stop "array sizes do not match" + end if + if(size(P,dim=2)/=size(Q,dim=2)) then + error stop "array sizes do not match" + end if + if (present(W)) then + if (size(W) /= size(P,dim=2)) then + error stop "array sizes do not match" + end if + end if + d = size(P,dim=1) + N = size(P,dim=2) + + ! Compute centroid + sum_w = 1.0_${k}$ / N + if(present(W)) sum_w = 1.0_${k}$ / stdlib_sum(W) + + #:if t.startswith('complex') + allocate(c_P(d), source=zero_c${k}$) + allocate(c_Q(d), source=zero_c${k}$) + #:else + allocate(c_P(d), source=zero_${k}$) + allocate(c_Q(d), source=zero_${k}$) + #:endif + + if(present(W)) then + do k = 1, N + c_P(1:d) = c_P(1:d) + P(1:d,k)*W(k) + c_Q(1:d) = c_Q(1:d) + Q(1:d,k)*W(k) + end do + else + do k = 1, N + c_P(1:d) = c_P(1:d) + P(1:d,k) + c_Q(1:d) = c_Q(1:d) + Q(1:d,k) + end do + end if + c_P = c_P * sum_w + c_Q = c_Q * sum_w + + ! Compute covariance matrix. + #:if t.startswith('complex') + allocate(covariance(d,d), source=zero_c${k}$) + #:else + allocate(covariance(d,d), source=zero_${k}$) + #:endif + variance_p = zero_${k}$ + + if (present(W)) then + do i = 1, N + do j = 1, d + vq = Q(j,i) - c_Q(j) + do k = 1, d + vp = P(k,i) - c_P(k) + covariance(k,j) = covariance(k,j) + vp * vq * W(i) + end do + end do + do k = 1, d + vp = P(k,i) - c_P(k) + variance_p = variance_p + vp*vp * W(i) + end do + end do + else + do i = 1, N + do j = 1, d + vq = Q(j,i) - c_Q(j) + do k = 1, d + vp = P(k,i) - c_P(k) + covariance(k,j) = covariance(k,j) + vp * vq + end do + end do + do k = 1, d + vp = P(k,i) - c_P(k) + variance_p = variance_p + vp*vp + end do + end do + end if + covariance = covariance * sum_w + variance_p = variance_p * sum_w + + #:if t.startswith('complex') + allocate(U(d,d), source=zero_c${k}$) + allocate(Vt(d,d), source=zero_c${k}$) + #:else + allocate(U(d,d), source=zero_${k}$) + allocate(Vt(d,d), source=zero_${k}$) + #:endif + allocate(S(d), source=zero_${k}$) + + ! SVD + call svd(covariance, S, U, Vt) + + #:if t.startswith('complex') + allocate(B(d,d), source=zero_c${k}$) + #:else + allocate(B(d,d), source=zero_${k}$) + #:endif + + ! Validate right-handed coordinate system + R = matmul(U, Vt) + B = eye(d,d) + #:if t.startswith('complex') + detR = det(R) + B(d,d) = detR / abs(detR) + #:else + B(d,d) = sign(one_${k}$, det(R)) + #:endif + + ! Optimal rotation matrix. + R = matmul(U,matmul(B,Vt)) + + ! Scaling factor + c = (sum(S(1:d-1))+B(d,d)*S(d)) /variance_p + if (.not. scale) c = 1.0_${k}$ + + ! Translation vector + t = c_P - c*matmul(R , c_Q ) + + #:if t.startswith('complex') + allocate(vec(d), source=zero_c${k}$) + #:else + allocate(vec(d), source=zero_${k}$) + #:endif + rmsd = 0.0_${k}$ + do i = 1, N + vec(1:d) = c * matmul(R , Q(1:d,i) ) + vec(1:d) = vec(1:d) + t(1:d) + vec(1:d) = vec(1:d) - P(1:d,i) + rmsd = rmsd + dot_product(vec,vec) + end do + rmsd = sqrt( rmsd * sum_w ) + end subroutine + #:endif + #:endfor +end submodule stdlib_spatial_kabsch From 655af91ab7efe8484a9b0b27cf4fc490c6dbfe14 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Thu, 5 Feb 2026 23:28:22 +0530 Subject: [PATCH 02/29] fixed the logic issue --- example/spatial/example_kabsch.f90 | 23 +++++---- src/spatial/stdlib_spatial_kabsch.fypp | 65 +++++++++++++------------- 2 files changed, 44 insertions(+), 44 deletions(-) diff --git a/example/spatial/example_kabsch.f90 b/example/spatial/example_kabsch.f90 index de364b170..dcbd3acb5 100644 --- a/example/spatial/example_kabsch.f90 +++ b/example/spatial/example_kabsch.f90 @@ -1,6 +1,7 @@ program example_kabsch use stdlib_linalg_constants, only: dp use stdlib_spatial + use stdlib_math, only: all_close implicit none integer, parameter :: d = 3, N = 4 @@ -12,19 +13,15 @@ program example_kabsch integer :: i ! Reference point set P (columns are points) - P(:,1) = [0.0_dp, 0.0_dp, 0.0_dp] - P(:,2) = [1.0_dp, 0.0_dp, 0.0_dp] - P(:,3) = [0.0_dp, 1.0_dp, 0.0_dp] - P(:,4) = [0.0_dp, 0.0_dp, 1.0_dp] +Q(:,1) = [0.0_dp, 0.0_dp, 0.0_dp] +Q(:,2) = [1.0_dp, 0.0_dp, 0.0_dp] +Q(:,3) = [0.0_dp, 1.0_dp, 0.0_dp] +Q(:,4) = [0.0_dp, 0.0_dp, 1.0_dp] - ! Target point set Q - ! Here Q is a rotated + translated + scaled version of P - ! Target point set Q = 2 * Rz * P + [1, 2, 3] - -Q(:,1) = [1.0_dp, 2.0_dp, 3.0_dp] ! P1 = (0,0,0) -Q(:,2) = [1.0_dp, 4.0_dp, 3.0_dp] ! P2 = (1,0,0) -> (0,1,0) -Q(:,3) = [-1.0_dp, 2.0_dp, 3.0_dp] ! P3 = (0,1,0) -> (-1,0,0) -Q(:,4) = [1.0_dp, 2.0_dp, 5.0_dp] ! P4 = (0,0,1) +P(:,1) = [2.0_dp, 0.0_dp, 1.0_dp] +P(:,2) = [2.0_dp, 1.0_dp, 1.0_dp] +P(:,3) = [1.0_dp, 0.0_dp, 2.0_dp] +P(:,4) = [3.0_dp, 0.0_dp, 1.0_dp] scale = .true. @@ -63,6 +60,8 @@ program example_kabsch print "(4F10.5)", Q_1(i,:) end do + print*, "Check P and RQ + t: ", all_close(P, Q_1, rel_tol=1.0e-6_dp, abs_tol=1.0e-12_dp) + print *, "" print *, "Scale factor c:", c print *, "RMSD:", rmsd diff --git a/src/spatial/stdlib_spatial_kabsch.fypp b/src/spatial/stdlib_spatial_kabsch.fypp index 1e2a40c73..63d5149c5 100644 --- a/src/spatial/stdlib_spatial_kabsch.fypp +++ b/src/spatial/stdlib_spatial_kabsch.fypp @@ -18,7 +18,7 @@ contains real(${k}$), intent(in), optional :: W(:) ! Internal variables. - integer(ilp) :: i, j, k, d, N + integer(ilp) :: i, j, point, d, N ${t}$, allocatable :: c_P(:), c_Q(:), covariance(:,:), U(:,:), Vt(:,:), B(:,:), vec(:) real(${k}$) :: sum_w, vp, vq, variance_p real(${k}$), allocatable :: S(:) @@ -41,9 +41,8 @@ contains d = size(P,dim=1) N = size(P,dim=2) - ! Compute centroid - sum_w = 1.0_${k}$ / N - if(present(W)) sum_w = 1.0_${k}$ / stdlib_sum(W) + sum_w = one_${k}$ / N + if(present(W)) sum_w = one_${k}$ / stdlib_sum(W) #:if t.startswith('complex') allocate(c_P(d), source=zero_c${k}$) @@ -53,15 +52,16 @@ contains allocate(c_Q(d), source=zero_${k}$) #:endif + ! Compute centroid if(present(W)) then - do k = 1, N - c_P(1:d) = c_P(1:d) + P(1:d,k)*W(k) - c_Q(1:d) = c_Q(1:d) + Q(1:d,k)*W(k) + do point = 1, N + c_P(1:d) = c_P(1:d) + P(1:d,point)*W(point) + c_Q(1:d) = c_Q(1:d) + Q(1:d,point)*W(point) end do else - do k = 1, N - c_P(1:d) = c_P(1:d) + P(1:d,k) - c_Q(1:d) = c_Q(1:d) + Q(1:d,k) + do point = 1, N + c_P(1:d) = c_P(1:d) + P(1:d,point) + c_Q(1:d) = c_Q(1:d) + Q(1:d,point) end do end if c_P = c_P * sum_w @@ -76,37 +76,38 @@ contains variance_p = zero_${k}$ if (present(W)) then - do i = 1, N + do point = 1, N + do i = 1, d + vp = P(i,point) - c_P(i) + variance_p = variance_p + vp*vp * W(point) + end do do j = 1, d - vq = Q(j,i) - c_Q(j) - do k = 1, d - vp = P(k,i) - c_P(k) - covariance(k,j) = covariance(k,j) + vp * vq * W(i) + vq = Q(j,point) - c_Q(j) + do i = 1, d + vp = P(i,point) - c_P(i) + covariance(i,j) = covariance(i,j) + vp*vq * W(point) end do end do - do k = 1, d - vp = P(k,i) - c_P(k) - variance_p = variance_p + vp*vp * W(i) - end do end do else - do i = 1, N + do point = 1, N + do i = 1, d + vp = P(i,point) - c_P(i) + variance_p = variance_p + vp*vp + end do do j = 1, d - vq = Q(j,i) - c_Q(j) - do k = 1, d - vp = P(k,i) - c_P(k) - covariance(k,j) = covariance(k,j) + vp * vq + vq = Q(j,point) - c_Q(j) + do i = 1, d + vp = P(i,point) - c_P(i) + covariance(i,j) = covariance(i,j) + vp*vq end do end do - do k = 1, d - vp = P(k,i) - c_P(k) - variance_p = variance_p + vp*vp - end do end do end if + covariance = covariance * sum_w variance_p = variance_p * sum_w - + #:if t.startswith('complex') allocate(U(d,d), source=zero_c${k}$) allocate(Vt(d,d), source=zero_c${k}$) @@ -139,8 +140,8 @@ contains R = matmul(U,matmul(B,Vt)) ! Scaling factor - c = (sum(S(1:d-1))+B(d,d)*S(d)) /variance_p - if (.not. scale) c = 1.0_${k}$ + c = variance_p / (sum(S(1:d-1))+B(d,d)*S(d)) + if (.not. scale) c = one_${k}$ ! Translation vector t = c_P - c*matmul(R , c_Q ) @@ -150,7 +151,7 @@ contains #:else allocate(vec(d), source=zero_${k}$) #:endif - rmsd = 0.0_${k}$ + rmsd = zero_${k}$ do i = 1, N vec(1:d) = c * matmul(R , Q(1:d,i) ) vec(1:d) = vec(1:d) + t(1:d) From 60818bb235a157697f894324d017e9d988eb80ab Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Fri, 6 Feb 2026 19:17:48 +0530 Subject: [PATCH 03/29] add: comments and complex support --- src/spatial/stdlib_spatial.fypp | 27 +++++++++-- src/spatial/stdlib_spatial_kabsch.fypp | 62 +++++++++++++++++++------- 2 files changed, 70 insertions(+), 19 deletions(-) diff --git a/src/spatial/stdlib_spatial.fypp b/src/spatial/stdlib_spatial.fypp index eb72f7f56..3f743c53b 100644 --- a/src/spatial/stdlib_spatial.fypp +++ b/src/spatial/stdlib_spatial.fypp @@ -12,19 +12,38 @@ module stdlib_spatial public :: kabsch interface kabsch + !----------------------------------------------------------------------- + !> Compute the optimal similarity transform (Kabsch–Umeyama): + !> + !> P ≈ c * R * Q + t + !> + !> where: + !> - R is an orthogonal rotation matrix + !> - c is an optional scale factor + !> - t is a translation vector + !> + !> The transformation minimizes the RMSD between corresponding columns + !> of P and Q, optionally using weights. + !----------------------------------------------------------------------- #:for k, t, s in (KINDS_TYPES) - #:if k!="xdp" - module subroutine kabsch_${s}$(P, Q, scale, R, t, c, rmsd, W) + module subroutine kabsch_${s}$(P, Q, R, t, c, rmsd, W, scale) + !> Reference point set (d × N) ${t}$, intent(in) :: P(:, :) + !> Target point set (d × N) ${t}$, intent(in) :: Q(:, :) - logical, intent(in) :: scale + !> Optimal rotation matrix (d × d) ${t}$, intent(out) :: R(:,:) + !> Translation vector (d) ${t}$, intent(out) :: t(:) + !> Scale factor real(${k}$), intent(out) :: c + !> Root-mean-square deviation real(${k}$), intent(out) :: rmsd + !> Optional weights real(${k}$), intent(in), optional :: W(:) + !> Enable scaling + logical, intent(in), optional :: scale end subroutine - #:endif #:endfor end interface end module stdlib_spatial diff --git a/src/spatial/stdlib_spatial_kabsch.fypp b/src/spatial/stdlib_spatial_kabsch.fypp index 63d5149c5..d5df39661 100644 --- a/src/spatial/stdlib_spatial_kabsch.fypp +++ b/src/spatial/stdlib_spatial_kabsch.fypp @@ -6,26 +6,36 @@ submodule(stdlib_spatial) stdlib_spatial_kabsch contains #:for k, t, s in (KINDS_TYPES) - #:if k!="xdp" - module subroutine kabsch_${s}$(P, Q, scale, R, t, c, rmsd, W) + module subroutine kabsch_${s}$(P, Q, R, t, c, rmsd, W, scale) + !> Reference point set (d × N) ${t}$, intent(in) :: P(:, :) + !> Target point set (d × N) ${t}$, intent(in) :: Q(:, :) - logical, intent(in) :: scale + !> Optimal rotation matrix (d × d) ${t}$, intent(out) :: R(:,:) + !> Translation vector (d) ${t}$, intent(out) :: t(:) + !> Scale factor real(${k}$), intent(out) :: c + !> Root-mean-square deviation real(${k}$), intent(out) :: rmsd + !> Optional weights real(${k}$), intent(in), optional :: W(:) + !> Enable scaling + logical, intent(in), optional :: scale ! Internal variables. integer(ilp) :: i, j, point, d, N ${t}$, allocatable :: c_P(:), c_Q(:), covariance(:,:), U(:,:), Vt(:,:), B(:,:), vec(:) - real(${k}$) :: sum_w, vp, vq, variance_p + ${t}$ :: vp, vq + real(${k}$) :: sum_w, variance_p real(${k}$), allocatable :: S(:) + logical :: scale_ #:if t.startswith('complex') ${t}$ :: detR #:endif - + + ! Dimension checks if(size(P,dim=1)/=size(Q,dim=1) .or. size(P,dim=1)/=size(R,dim=1) .or. size(P,dim=1)/=size(R,dim=2) & .or. size(P,dim=1)/=size(t)) then error stop "array sizes do not match" @@ -40,6 +50,8 @@ contains end if d = size(P,dim=1) N = size(P,dim=2) + scale_ = .true. + if(present(scale)) scale_ = scale sum_w = one_${k}$ / N if(present(W)) sum_w = one_${k}$ / stdlib_sum(W) @@ -52,7 +64,7 @@ contains allocate(c_Q(d), source=zero_${k}$) #:endif - ! Compute centroid + ! Compute centroids of P and Q if(present(W)) then do point = 1, N c_P(1:d) = c_P(1:d) + P(1:d,point)*W(point) @@ -67,7 +79,7 @@ contains c_P = c_P * sum_w c_Q = c_Q * sum_w - ! Compute covariance matrix. + ! Compute covariance matrix H = (P - c_P)^T (Q - c_Q) and variance of P #:if t.startswith('complex') allocate(covariance(d,d), source=zero_c${k}$) #:else @@ -79,13 +91,21 @@ contains do point = 1, N do i = 1, d vp = P(i,point) - c_P(i) + #:if t.startswith('complex') + variance_p = variance_p + conjg(vp)*vp * W(point) + #:else variance_p = variance_p + vp*vp * W(point) + #:endif end do do j = 1, d vq = Q(j,point) - c_Q(j) do i = 1, d vp = P(i,point) - c_P(i) + #:if t.startswith('complex') + covariance(i,j) = covariance(i,j) + vp*conjg(vq) * W(point) + #:else covariance(i,j) = covariance(i,j) + vp*vq * W(point) + #:endif end do end do end do @@ -93,13 +113,21 @@ contains do point = 1, N do i = 1, d vp = P(i,point) - c_P(i) + #:if t.startswith('complex') + variance_p = variance_p + conjg(vp)*vp + #:else variance_p = variance_p + vp*vp + #:endif end do do j = 1, d vq = Q(j,point) - c_Q(j) do i = 1, d vp = P(i,point) - c_P(i) + #:if t.startswith('complex') + covariance(i,j) = covariance(i,j) + vp*conjg(vq) + #:else covariance(i,j) = covariance(i,j) + vp*vq + #:endif end do end do end do @@ -116,10 +144,10 @@ contains allocate(Vt(d,d), source=zero_${k}$) #:endif allocate(S(d), source=zero_${k}$) - - ! SVD + + ! SVD of covariance matrix H -> H = U * S * Vt call svd(covariance, S, U, Vt) - + #:if t.startswith('complex') allocate(B(d,d), source=zero_c${k}$) #:else @@ -141,11 +169,12 @@ contains ! Scaling factor c = variance_p / (sum(S(1:d-1))+B(d,d)*S(d)) - if (.not. scale) c = one_${k}$ + if (.not. scale_) c = one_${k}$ ! Translation vector t = c_P - c*matmul(R , c_Q ) + ! Compute RMSD #:if t.startswith('complex') allocate(vec(d), source=zero_c${k}$) #:else @@ -153,13 +182,16 @@ contains #:endif rmsd = zero_${k}$ do i = 1, N - vec(1:d) = c * matmul(R , Q(1:d,i) ) + vec(1:d) = c * matmul(R , Q(1:d,i)) vec(1:d) = vec(1:d) + t(1:d) vec(1:d) = vec(1:d) - P(1:d,i) - rmsd = rmsd + dot_product(vec,vec) + #:if t.startswith('complex') + rmsd = rmsd + real(dot_product(vec, vec), kind=${k}$) + #:else + rmsd = rmsd + dot_product(vec, vec) + #:endif end do - rmsd = sqrt( rmsd * sum_w ) + rmsd = sqrt(rmsd * sum_w) end subroutine - #:endif #:endfor end submodule stdlib_spatial_kabsch From 400b85d2b5e6c48cf2dd94ca7f786fca43e47781 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Sat, 7 Feb 2026 11:26:30 +0530 Subject: [PATCH 04/29] removed complex kinds --- example/spatial/example_kabsch.f90 | 123 ++++++++++++++++++------- src/spatial/stdlib_spatial.fypp | 5 +- src/spatial/stdlib_spatial_kabsch.fypp | 61 +----------- 3 files changed, 95 insertions(+), 94 deletions(-) diff --git a/example/spatial/example_kabsch.f90 b/example/spatial/example_kabsch.f90 index dcbd3acb5..5e6ede868 100644 --- a/example/spatial/example_kabsch.f90 +++ b/example/spatial/example_kabsch.f90 @@ -1,69 +1,122 @@ -program example_kabsch +program example_kabsch_real use stdlib_linalg_constants, only: dp use stdlib_spatial use stdlib_math, only: all_close + use stdlib_linalg, only: svd, det implicit none integer, parameter :: d = 3, N = 4 real(dp) :: P(d, N), Q(d, N), Q_1(d, N) - real(dp) :: R(d, d), t(d) - real(dp) :: c, rmsd - logical :: scale + real(dp) :: R(d, d), R_original(d, d) + real(dp) :: t(d), t_original(d) + real(dp) :: c, c_original + real(dp) :: rmsd + + integer :: i, j + real(dp) :: r1 + + call random_seed() + + ! ---------------------------- + ! Random reference points Q + ! ---------------------------- + do j = 1, N + do i = 1, d + call random_number(r1) + Q(i,j) = r1 + end do + end do - integer :: i + ! ------------------------------------------------ + ! Random proper rotation matrix R_original + ! Constructed via SVD: R = U * V^T + ! ------------------------------------------------ + do i = 1, d + do j = 1, d + call random_number(r1) + R_original(i,j) = r1 + end do + end do - ! Reference point set P (columns are points) -Q(:,1) = [0.0_dp, 0.0_dp, 0.0_dp] -Q(:,2) = [1.0_dp, 0.0_dp, 0.0_dp] -Q(:,3) = [0.0_dp, 1.0_dp, 0.0_dp] -Q(:,4) = [0.0_dp, 0.0_dp, 1.0_dp] + block + real(dp) :: U(d,d), Vt(d,d), S(d) -P(:,1) = [2.0_dp, 0.0_dp, 1.0_dp] -P(:,2) = [2.0_dp, 1.0_dp, 1.0_dp] -P(:,3) = [1.0_dp, 0.0_dp, 2.0_dp] -P(:,4) = [3.0_dp, 0.0_dp, 1.0_dp] + call svd(R_original, S, U, Vt) + R_original = matmul(U, Vt) + ! Enforce det = +1 (no reflection) + if (det(R_original) < 0.0_dp) then + U(:,d) = -U(:,d) + R_original = matmul(U, Vt) + end if + end block + R_original(:,1) = -R_original(:,1) - scale = .true. + print *, "det(R_original):", det(R_original) - call kabsch(P, Q, scale, R, t, c, rmsd) + ! ---------------------------- + ! Random scale and translation + ! ---------------------------- + call random_number(r1) + c_original = 0.5_dp + 2.0_dp * r1 - print *, "" - print *, "Matrix P:" do i = 1, d - print "(4F10.5)", P(i,:) + call random_number(r1) + t_original(i) = r1 end do + ! ---------------------------- + ! Construct P = c*R*Q + t + ! ---------------------------- + do j = 1, N + P(:,j) = c_original * matmul(R_original, Q(:,j)) + t_original + end do + + ! ---------------------------- + ! Call Kabsch–Umeyama + ! ---------------------------- + call kabsch(P, Q, R, t, c, rmsd) + print *, "" - print *, "Matrix Q:" + print *, "Original rotation R_original:" do i = 1, d - print "(4F10.5)", Q(i,:) + print *, R_original(i,:) end do - print *, "Rotation matrix R:" + print *, "" + print *, "Recovered rotation R:" do i = 1, d - print "(3F10.5)", R(i,:) + print *, R(i,:) end do print *, "" - print *, "Translation vector t:" - print "(3F10.5)", t + print *, "Original translation t_original:" + print *, t_original - ! Apply the transformation: Q_1 = c * R * Q + t - do i = 1, N - Q_1(:,i) = c * matmul(R, Q(:,i)) + t + print *, "" + print *, "Recovered translation t:" + print *, t + + ! ---------------------------- + ! Apply recovered transform + ! ---------------------------- + do j = 1, N + Q_1(:,j) = c * matmul(R, Q(:,j)) + t end do print *, "" - print *, "Aligned Q (should match P):" - do i = 1, d - print "(4F10.5)", Q_1(i,:) - end do + print *, "Check P ≈ c*R*Q + t: ", & + all_close(P, Q_1, rel_tol=1.0e-10_dp, abs_tol=1.0e-12_dp) - print*, "Check P and RQ + t: ", all_close(P, Q_1, rel_tol=1.0e-6_dp, abs_tol=1.0e-12_dp) + print *, "" + print *, "Original scale c_original:", c_original + print *, "Recovered scale c:", c print *, "" - print *, "Scale factor c:", c print *, "RMSD:", rmsd -end program example_kabsch + print *, "" + print *, "Check R ≈ R_original: ", & + all_close(R, R_original, rel_tol=1.0e-10_dp, abs_tol=1.0e-12_dp) + +end program example_kabsch_real diff --git a/src/spatial/stdlib_spatial.fypp b/src/spatial/stdlib_spatial.fypp index 3f743c53b..93051a122 100644 --- a/src/spatial/stdlib_spatial.fypp +++ b/src/spatial/stdlib_spatial.fypp @@ -1,7 +1,6 @@ #:include "common.fypp" #:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) -#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) -#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES +#:set KINDS_TYPES = R_KINDS_TYPES module stdlib_spatial use stdlib_linalg_constants use stdlib_constants @@ -36,7 +35,7 @@ module stdlib_spatial !> Translation vector (d) ${t}$, intent(out) :: t(:) !> Scale factor - real(${k}$), intent(out) :: c + ${t}$, intent(out) :: c !> Root-mean-square deviation real(${k}$), intent(out) :: rmsd !> Optional weights diff --git a/src/spatial/stdlib_spatial_kabsch.fypp b/src/spatial/stdlib_spatial_kabsch.fypp index d5df39661..5525a913c 100644 --- a/src/spatial/stdlib_spatial_kabsch.fypp +++ b/src/spatial/stdlib_spatial_kabsch.fypp @@ -1,7 +1,6 @@ #:include "common.fypp" #:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) -#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) -#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES +#:set KINDS_TYPES = R_KINDS_TYPES submodule(stdlib_spatial) stdlib_spatial_kabsch contains @@ -16,7 +15,7 @@ contains !> Translation vector (d) ${t}$, intent(out) :: t(:) !> Scale factor - real(${k}$), intent(out) :: c + ${t}$, intent(out) :: c !> Root-mean-square deviation real(${k}$), intent(out) :: rmsd !> Optional weights @@ -31,9 +30,6 @@ contains real(${k}$) :: sum_w, variance_p real(${k}$), allocatable :: S(:) logical :: scale_ - #:if t.startswith('complex') - ${t}$ :: detR - #:endif ! Dimension checks if(size(P,dim=1)/=size(Q,dim=1) .or. size(P,dim=1)/=size(R,dim=1) .or. size(P,dim=1)/=size(R,dim=2) & @@ -56,13 +52,8 @@ contains sum_w = one_${k}$ / N if(present(W)) sum_w = one_${k}$ / stdlib_sum(W) - #:if t.startswith('complex') - allocate(c_P(d), source=zero_c${k}$) - allocate(c_Q(d), source=zero_c${k}$) - #:else allocate(c_P(d), source=zero_${k}$) allocate(c_Q(d), source=zero_${k}$) - #:endif ! Compute centroids of P and Q if(present(W)) then @@ -79,33 +70,21 @@ contains c_P = c_P * sum_w c_Q = c_Q * sum_w - ! Compute covariance matrix H = (P - c_P)^T (Q - c_Q) and variance of P - #:if t.startswith('complex') - allocate(covariance(d,d), source=zero_c${k}$) - #:else + ! Compute covariance matrix H = (P - c_P) * (Q - c_Q)^T and variance of P allocate(covariance(d,d), source=zero_${k}$) - #:endif variance_p = zero_${k}$ if (present(W)) then do point = 1, N do i = 1, d vp = P(i,point) - c_P(i) - #:if t.startswith('complex') - variance_p = variance_p + conjg(vp)*vp * W(point) - #:else variance_p = variance_p + vp*vp * W(point) - #:endif end do do j = 1, d vq = Q(j,point) - c_Q(j) do i = 1, d vp = P(i,point) - c_P(i) - #:if t.startswith('complex') - covariance(i,j) = covariance(i,j) + vp*conjg(vq) * W(point) - #:else covariance(i,j) = covariance(i,j) + vp*vq * W(point) - #:endif end do end do end do @@ -113,21 +92,13 @@ contains do point = 1, N do i = 1, d vp = P(i,point) - c_P(i) - #:if t.startswith('complex') - variance_p = variance_p + conjg(vp)*vp - #:else variance_p = variance_p + vp*vp - #:endif end do do j = 1, d vq = Q(j,point) - c_Q(j) do i = 1, d vp = P(i,point) - c_P(i) - #:if t.startswith('complex') - covariance(i,j) = covariance(i,j) + vp*conjg(vq) - #:else covariance(i,j) = covariance(i,j) + vp*vq - #:endif end do end do end do @@ -136,33 +107,19 @@ contains covariance = covariance * sum_w variance_p = variance_p * sum_w - #:if t.startswith('complex') - allocate(U(d,d), source=zero_c${k}$) - allocate(Vt(d,d), source=zero_c${k}$) - #:else allocate(U(d,d), source=zero_${k}$) allocate(Vt(d,d), source=zero_${k}$) - #:endif allocate(S(d), source=zero_${k}$) ! SVD of covariance matrix H -> H = U * S * Vt call svd(covariance, S, U, Vt) - #:if t.startswith('complex') - allocate(B(d,d), source=zero_c${k}$) - #:else allocate(B(d,d), source=zero_${k}$) - #:endif ! Validate right-handed coordinate system R = matmul(U, Vt) B = eye(d,d) - #:if t.startswith('complex') - detR = det(R) - B(d,d) = detR / abs(detR) - #:else B(d,d) = sign(one_${k}$, det(R)) - #:endif ! Optimal rotation matrix. R = matmul(U,matmul(B,Vt)) @@ -172,24 +129,16 @@ contains if (.not. scale_) c = one_${k}$ ! Translation vector - t = c_P - c*matmul(R , c_Q ) + t = c_P - c*matmul(R , c_Q) ! Compute RMSD - #:if t.startswith('complex') - allocate(vec(d), source=zero_c${k}$) - #:else allocate(vec(d), source=zero_${k}$) - #:endif rmsd = zero_${k}$ do i = 1, N - vec(1:d) = c * matmul(R , Q(1:d,i)) + vec(1:d) = c * matmul(R, Q(1:d,i)) vec(1:d) = vec(1:d) + t(1:d) vec(1:d) = vec(1:d) - P(1:d,i) - #:if t.startswith('complex') - rmsd = rmsd + real(dot_product(vec, vec), kind=${k}$) - #:else rmsd = rmsd + dot_product(vec, vec) - #:endif end do rmsd = sqrt(rmsd * sum_w) end subroutine From 5208695dd031c55708aeb89b31856306b39bf896 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Sat, 7 Feb 2026 12:41:37 +0530 Subject: [PATCH 05/29] add complex types back, correct reflection --- example/spatial/example_kabsch.f90 | 72 +++++++++++++------------- src/spatial/stdlib_spatial.fypp | 5 +- src/spatial/stdlib_spatial_kabsch.fypp | 63 ++++++++++++++++++---- 3 files changed, 92 insertions(+), 48 deletions(-) diff --git a/example/spatial/example_kabsch.f90 b/example/spatial/example_kabsch.f90 index 5e6ede868..b8deea26f 100644 --- a/example/spatial/example_kabsch.f90 +++ b/example/spatial/example_kabsch.f90 @@ -1,4 +1,4 @@ -program example_kabsch_real +program example_kabsch use stdlib_linalg_constants, only: dp use stdlib_spatial use stdlib_math, only: all_close @@ -6,76 +6,78 @@ program example_kabsch_real implicit none integer, parameter :: d = 3, N = 4 - real(dp) :: P(d, N), Q(d, N), Q_1(d, N) - real(dp) :: R(d, d), R_original(d, d) - real(dp) :: t(d), t_original(d) - real(dp) :: c, c_original - real(dp) :: rmsd + complex(dp) :: P_original(d, N), Q_original(d, N), P_recovered(d, N) + complex(dp) :: R_recovered(d, d), R_original(d, d) + complex(dp) :: t_recovered(d), t_original(d) + complex(dp) :: c_recovered, c_original + real(dp) :: rmsd integer :: i, j - real(dp) :: r1 + real(dp) :: r1, r2 call random_seed() ! ---------------------------- - ! Random reference points Q + ! Random complex reference points Q ! ---------------------------- do j = 1, N do i = 1, d call random_number(r1) - Q(i,j) = r1 + call random_number(r2) + Q_original(i,j) = cmplx(r1, r2, kind=dp) end do end do ! ------------------------------------------------ - ! Random proper rotation matrix R_original - ! Constructed via SVD: R = U * V^T + ! Random complex unitary matrix R_original + ! Constructed via SVD: R = U * V^H ! ------------------------------------------------ do i = 1, d do j = 1, d call random_number(r1) - R_original(i,j) = r1 + call random_number(r2) + R_original(i,j) = cmplx(r1, r2, kind=dp) end do end do block - real(dp) :: U(d,d), Vt(d,d), S(d) + complex(dp) :: U(d,d), Vt(d,d) + real(dp) :: S(d) call svd(R_original, S, U, Vt) - R_original = matmul(U, Vt) + R_original = matmul(U, Vt) ! unitary - ! Enforce det = +1 (no reflection) - if (det(R_original) < 0.0_dp) then - U(:,d) = -U(:,d) - R_original = matmul(U, Vt) - end if end block - R_original(:,1) = -R_original(:,1) + + ! Complex reflection of pi phase flip + R_original(:, 1) = -R_original(:,1) - print *, "det(R_original):", det(R_original) + print *, "abs(det(R_original)):", abs(det(R_original)) ! ---------------------------- ! Random scale and translation ! ---------------------------- call random_number(r1) - c_original = 0.5_dp + 2.0_dp * r1 + call random_number(r2) + c_original = cmplx(r1, r2) do i = 1, d call random_number(r1) - t_original(i) = r1 + call random_number(r2) + t_original(i) = cmplx(r1, r2, kind=dp) end do ! ---------------------------- ! Construct P = c*R*Q + t ! ---------------------------- do j = 1, N - P(:,j) = c_original * matmul(R_original, Q(:,j)) + t_original + P_original(:,j) = c_original * matmul(R_original, Q_original(:,j)) + t_original end do ! ---------------------------- - ! Call Kabsch–Umeyama + ! Call complex Kabsch–Umeyama ! ---------------------------- - call kabsch(P, Q, R, t, c, rmsd) + call kabsch(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) print *, "" print *, "Original rotation R_original:" @@ -86,7 +88,7 @@ program example_kabsch_real print *, "" print *, "Recovered rotation R:" do i = 1, d - print *, R(i,:) + print *, R_recovered(i,:) end do print *, "" @@ -95,28 +97,28 @@ program example_kabsch_real print *, "" print *, "Recovered translation t:" - print *, t + print *, t_recovered ! ---------------------------- ! Apply recovered transform ! ---------------------------- do j = 1, N - Q_1(:,j) = c * matmul(R, Q(:,j)) + t + P_recovered(:,j) = c_recovered * matmul(R_recovered, Q_original(:,j)) + t_recovered end do print *, "" - print *, "Check P ≈ c*R*Q + t: ", & - all_close(P, Q_1, rel_tol=1.0e-10_dp, abs_tol=1.0e-12_dp) + print *, "Check P_original ≈ P_recovered: ", & + all_close(P_original, P_recovered) print *, "" print *, "Original scale c_original:", c_original - print *, "Recovered scale c:", c + print *, "Recovered scale c:", c_recovered print *, "" print *, "RMSD:", rmsd print *, "" - print *, "Check R ≈ R_original: ", & - all_close(R, R_original, rel_tol=1.0e-10_dp, abs_tol=1.0e-12_dp) + print *, "Check c_recovered * R_recovered ≈ c_original * R_original: ", & + all_close(c_recovered*R_recovered, c_original*R_original) -end program example_kabsch_real +end program example_kabsch \ No newline at end of file diff --git a/src/spatial/stdlib_spatial.fypp b/src/spatial/stdlib_spatial.fypp index 93051a122..f6451a340 100644 --- a/src/spatial/stdlib_spatial.fypp +++ b/src/spatial/stdlib_spatial.fypp @@ -1,6 +1,7 @@ #:include "common.fypp" #:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) -#:set KINDS_TYPES = R_KINDS_TYPES +#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) +#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES module stdlib_spatial use stdlib_linalg_constants use stdlib_constants @@ -45,4 +46,4 @@ module stdlib_spatial end subroutine #:endfor end interface -end module stdlib_spatial +end module stdlib_spatial \ No newline at end of file diff --git a/src/spatial/stdlib_spatial_kabsch.fypp b/src/spatial/stdlib_spatial_kabsch.fypp index 5525a913c..64c6020d4 100644 --- a/src/spatial/stdlib_spatial_kabsch.fypp +++ b/src/spatial/stdlib_spatial_kabsch.fypp @@ -1,6 +1,7 @@ #:include "common.fypp" #:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) -#:set KINDS_TYPES = R_KINDS_TYPES +#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) +#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES submodule(stdlib_spatial) stdlib_spatial_kabsch contains @@ -30,6 +31,9 @@ contains real(${k}$) :: sum_w, variance_p real(${k}$), allocatable :: S(:) logical :: scale_ + #:if t.startswith('complex') + ${t}$ :: detR + #:endif ! Dimension checks if(size(P,dim=1)/=size(Q,dim=1) .or. size(P,dim=1)/=size(R,dim=1) .or. size(P,dim=1)/=size(R,dim=2) & @@ -52,8 +56,13 @@ contains sum_w = one_${k}$ / N if(present(W)) sum_w = one_${k}$ / stdlib_sum(W) + #:if t.startswith('complex') + allocate(c_P(d), source=zero_c${k}$) + allocate(c_Q(d), source=zero_c${k}$) + #:else allocate(c_P(d), source=zero_${k}$) allocate(c_Q(d), source=zero_${k}$) + #:endif ! Compute centroids of P and Q if(present(W)) then @@ -71,20 +80,32 @@ contains c_Q = c_Q * sum_w ! Compute covariance matrix H = (P - c_P) * (Q - c_Q)^T and variance of P + #:if t.startswith('complex') + allocate(covariance(d,d), source=zero_c${k}$) + #:else allocate(covariance(d,d), source=zero_${k}$) + #:endif variance_p = zero_${k}$ if (present(W)) then do point = 1, N do i = 1, d vp = P(i,point) - c_P(i) + #:if t.startswith('complex') + variance_p = variance_p + conjg(vp)*vp * W(point) + #:else variance_p = variance_p + vp*vp * W(point) + #:endif end do do j = 1, d vq = Q(j,point) - c_Q(j) do i = 1, d vp = P(i,point) - c_P(i) + #:if t.startswith('complex') + covariance(i,j) = covariance(i,j) + vp*conjg(vq) * W(point) + #:else covariance(i,j) = covariance(i,j) + vp*vq * W(point) + #:endif end do end do end do @@ -92,13 +113,21 @@ contains do point = 1, N do i = 1, d vp = P(i,point) - c_P(i) + #:if t.startswith('complex') + variance_p = variance_p + conjg(vp)*vp + #:else variance_p = variance_p + vp*vp + #:endif end do do j = 1, d vq = Q(j,point) - c_Q(j) do i = 1, d vp = P(i,point) - c_P(i) + #:if t.startswith('complex') + covariance(i,j) = covariance(i,j) + vp*conjg(vq) + #:else covariance(i,j) = covariance(i,j) + vp*vq + #:endif end do end do end do @@ -107,40 +136,52 @@ contains covariance = covariance * sum_w variance_p = variance_p * sum_w + #:if t.startswith('complex') + allocate(U(d,d), source=zero_c${k}$) + allocate(Vt(d,d), source=zero_c${k}$) + #:else allocate(U(d,d), source=zero_${k}$) allocate(Vt(d,d), source=zero_${k}$) + #:endif allocate(S(d), source=zero_${k}$) ! SVD of covariance matrix H -> H = U * S * Vt call svd(covariance, S, U, Vt) + #:if t.startswith('complex') + allocate(B(d,d), source=zero_c${k}$) + #:else allocate(B(d,d), source=zero_${k}$) - - ! Validate right-handed coordinate system - R = matmul(U, Vt) - B = eye(d,d) - B(d,d) = sign(one_${k}$, det(R)) + #:endif ! Optimal rotation matrix. - R = matmul(U,matmul(B,Vt)) + R = matmul(U, Vt) ! Scaling factor - c = variance_p / (sum(S(1:d-1))+B(d,d)*S(d)) + c = variance_p / (sum(S(1:d))) if (.not. scale_) c = one_${k}$ ! Translation vector - t = c_P - c*matmul(R , c_Q) + t = c_P - c*matmul(R , c_Q ) ! Compute RMSD + #:if t.startswith('complex') + allocate(vec(d), source=zero_c${k}$) + #:else allocate(vec(d), source=zero_${k}$) + #:endif rmsd = zero_${k}$ do i = 1, N - vec(1:d) = c * matmul(R, Q(1:d,i)) + vec(1:d) = c * matmul(R , Q(1:d,i)) vec(1:d) = vec(1:d) + t(1:d) vec(1:d) = vec(1:d) - P(1:d,i) + #:if t.startswith('complex') + rmsd = rmsd + real(dot_product(vec, vec), kind=${k}$) + #:else rmsd = rmsd + dot_product(vec, vec) + #:endif end do rmsd = sqrt(rmsd * sum_w) end subroutine #:endfor -end submodule stdlib_spatial_kabsch +end submodule stdlib_spatial_kabsch \ No newline at end of file From 25df8d279243c94bb5f42da46c979a1a5055aa5a Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Sat, 7 Feb 2026 17:17:58 +0530 Subject: [PATCH 06/29] modify example, add tests, use suffix in constants and move use stdlib det,svd,eye inside submodule file --- example/spatial/example_kabsch.f90 | 116 ++------------- src/spatial/stdlib_spatial.fypp | 3 +- src/spatial/stdlib_spatial_kabsch.fypp | 46 +++--- test/CMakeLists.txt | 1 + test/spatial/CMakeLists.txt | 7 + test/spatial/test_spatial_kabsch.fypp | 187 +++++++++++++++++++++++++ 6 files changed, 234 insertions(+), 126 deletions(-) create mode 100644 test/spatial/CMakeLists.txt create mode 100644 test/spatial/test_spatial_kabsch.fypp diff --git a/example/spatial/example_kabsch.f90 b/example/spatial/example_kabsch.f90 index b8deea26f..4ae37d794 100644 --- a/example/spatial/example_kabsch.f90 +++ b/example/spatial/example_kabsch.f90 @@ -1,124 +1,36 @@ program example_kabsch use stdlib_linalg_constants, only: dp - use stdlib_spatial - use stdlib_math, only: all_close - use stdlib_linalg, only: svd, det + use stdlib_spatial, only: kabsch implicit none - integer, parameter :: d = 3, N = 4 - complex(dp) :: P_original(d, N), Q_original(d, N), P_recovered(d, N) - complex(dp) :: R_recovered(d, d), R_original(d, d) - complex(dp) :: t_recovered(d), t_original(d) - complex(dp) :: c_recovered, c_original - real(dp) :: rmsd + integer, parameter :: d = 2, N = 3 + real(dp) :: P(d, N), Q(d, N), R(d, d), t(d), c, rmsd - integer :: i, j - real(dp) :: r1, r2 + integer :: i - call random_seed() + P(:,1) = [3.0_dp, -2.0_dp] + P(:,2) = [7.0_dp, 4.0_dp] + P(:,3) = [5.0_dp, 0.0_dp] - ! ---------------------------- - ! Random complex reference points Q - ! ---------------------------- - do j = 1, N - do i = 1, d - call random_number(r1) - call random_number(r2) - Q_original(i,j) = cmplx(r1, r2, kind=dp) - end do - end do - - ! ------------------------------------------------ - ! Random complex unitary matrix R_original - ! Constructed via SVD: R = U * V^H - ! ------------------------------------------------ - do i = 1, d - do j = 1, d - call random_number(r1) - call random_number(r2) - R_original(i,j) = cmplx(r1, r2, kind=dp) - end do - end do - - block - complex(dp) :: U(d,d), Vt(d,d) - real(dp) :: S(d) - - call svd(R_original, S, U, Vt) - R_original = matmul(U, Vt) ! unitary - - end block - - ! Complex reflection of pi phase flip - R_original(:, 1) = -R_original(:,1) - - print *, "abs(det(R_original)):", abs(det(R_original)) - - ! ---------------------------- - ! Random scale and translation - ! ---------------------------- - call random_number(r1) - call random_number(r2) - c_original = cmplx(r1, r2) - - do i = 1, d - call random_number(r1) - call random_number(r2) - t_original(i) = cmplx(r1, r2, kind=dp) - end do - - ! ---------------------------- - ! Construct P = c*R*Q + t - ! ---------------------------- - do j = 1, N - P_original(:,j) = c_original * matmul(R_original, Q_original(:,j)) + t_original - end do + Q(:,1) = [2.0_dp, 3.0_dp] + Q(:,2) = [-1.0_dp, 5.0_dp] + Q(:,3) = [1.0_dp, 4.0_dp] - ! ---------------------------- - ! Call complex Kabsch–Umeyama - ! ---------------------------- - call kabsch(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) - - print *, "" - print *, "Original rotation R_original:" - do i = 1, d - print *, R_original(i,:) - end do + call kabsch(P, Q, R, t, c, rmsd) print *, "" print *, "Recovered rotation R:" do i = 1, d - print *, R_recovered(i,:) + print *, R(i,:) end do - print *, "" - print *, "Original translation t_original:" - print *, t_original + print *, "Recovered scale c:", c print *, "" print *, "Recovered translation t:" - print *, t_recovered - - ! ---------------------------- - ! Apply recovered transform - ! ---------------------------- - do j = 1, N - P_recovered(:,j) = c_recovered * matmul(R_recovered, Q_original(:,j)) + t_recovered - end do - - print *, "" - print *, "Check P_original ≈ P_recovered: ", & - all_close(P_original, P_recovered) - - print *, "" - print *, "Original scale c_original:", c_original - print *, "Recovered scale c:", c_recovered + print *, t print *, "" print *, "RMSD:", rmsd - print *, "" - print *, "Check c_recovered * R_recovered ≈ c_original * R_original: ", & - all_close(c_recovered*R_recovered, c_original*R_original) - end program example_kabsch \ No newline at end of file diff --git a/src/spatial/stdlib_spatial.fypp b/src/spatial/stdlib_spatial.fypp index f6451a340..f0ee0f197 100644 --- a/src/spatial/stdlib_spatial.fypp +++ b/src/spatial/stdlib_spatial.fypp @@ -5,8 +5,7 @@ module stdlib_spatial use stdlib_linalg_constants use stdlib_constants - use stdlib_linalg, only: svd, det, eye, trace - use stdlib_intrinsics, only: stdlib_sum + use stdlib_error, only: error_stop implicit none private public :: kabsch diff --git a/src/spatial/stdlib_spatial_kabsch.fypp b/src/spatial/stdlib_spatial_kabsch.fypp index 64c6020d4..4e40788fb 100644 --- a/src/spatial/stdlib_spatial_kabsch.fypp +++ b/src/spatial/stdlib_spatial_kabsch.fypp @@ -3,6 +3,8 @@ #:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) #:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES submodule(stdlib_spatial) stdlib_spatial_kabsch + use stdlib_linalg, only: svd, det, eye + use stdlib_intrinsics, only: stdlib_sum contains #:for k, t, s in (KINDS_TYPES) @@ -38,14 +40,14 @@ contains ! Dimension checks if(size(P,dim=1)/=size(Q,dim=1) .or. size(P,dim=1)/=size(R,dim=1) .or. size(P,dim=1)/=size(R,dim=2) & .or. size(P,dim=1)/=size(t)) then - error stop "array sizes do not match" + call error_stop("array sizes do not match") end if if(size(P,dim=2)/=size(Q,dim=2)) then - error stop "array sizes do not match" + call error_stop("array sizes do not match") end if if (present(W)) then if (size(W) /= size(P,dim=2)) then - error stop "array sizes do not match" + call error_stop("array sizes do not match") end if end if d = size(P,dim=1) @@ -53,15 +55,15 @@ contains scale_ = .true. if(present(scale)) scale_ = scale - sum_w = one_${k}$ / N - if(present(W)) sum_w = one_${k}$ / stdlib_sum(W) + sum_w = one_${s}$ / N + if(present(W)) sum_w = one_${s}$ / stdlib_sum(W) #:if t.startswith('complex') - allocate(c_P(d), source=zero_c${k}$) - allocate(c_Q(d), source=zero_c${k}$) + allocate(c_P(d), source=zero_${s}$) + allocate(c_Q(d), source=zero_${s}$) #:else - allocate(c_P(d), source=zero_${k}$) - allocate(c_Q(d), source=zero_${k}$) + allocate(c_P(d), source=zero_${s}$) + allocate(c_Q(d), source=zero_${s}$) #:endif ! Compute centroids of P and Q @@ -81,11 +83,11 @@ contains ! Compute covariance matrix H = (P - c_P) * (Q - c_Q)^T and variance of P #:if t.startswith('complex') - allocate(covariance(d,d), source=zero_c${k}$) + allocate(covariance(d,d), source=zero_${s}$) #:else - allocate(covariance(d,d), source=zero_${k}$) + allocate(covariance(d,d), source=zero_${s}$) #:endif - variance_p = zero_${k}$ + variance_p = zero_${s}$ if (present(W)) then do point = 1, N @@ -137,11 +139,11 @@ contains variance_p = variance_p * sum_w #:if t.startswith('complex') - allocate(U(d,d), source=zero_c${k}$) - allocate(Vt(d,d), source=zero_c${k}$) + allocate(U(d,d), source=zero_${s}$) + allocate(Vt(d,d), source=zero_${s}$) #:else - allocate(U(d,d), source=zero_${k}$) - allocate(Vt(d,d), source=zero_${k}$) + allocate(U(d,d), source=zero_${s}$) + allocate(Vt(d,d), source=zero_${s}$) #:endif allocate(S(d), source=zero_${k}$) @@ -149,9 +151,9 @@ contains call svd(covariance, S, U, Vt) #:if t.startswith('complex') - allocate(B(d,d), source=zero_c${k}$) + allocate(B(d,d), source=zero_${s}$) #:else - allocate(B(d,d), source=zero_${k}$) + allocate(B(d,d), source=zero_${s}$) #:endif ! Optimal rotation matrix. @@ -159,18 +161,18 @@ contains ! Scaling factor c = variance_p / (sum(S(1:d))) - if (.not. scale_) c = one_${k}$ + if (.not. scale_) c = one_${s}$ ! Translation vector t = c_P - c*matmul(R , c_Q ) ! Compute RMSD #:if t.startswith('complex') - allocate(vec(d), source=zero_c${k}$) + allocate(vec(d), source=zero_${s}$) #:else - allocate(vec(d), source=zero_${k}$) + allocate(vec(d), source=zero_${s}$) #:endif - rmsd = zero_${k}$ + rmsd = zero_${s}$ do i = 1, N vec(1:d) = c * matmul(R , Q(1:d,i)) vec(1:d) = vec(1:d) + t(1:d) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index ea43bc8eb..caa46c947 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -49,6 +49,7 @@ endif() add_subdirectory(optval) add_subdirectory(selection) add_subdirectory(sorting) +add_subdirectory(spatial) add_subdirectory(specialfunctions) if (STDLIB_STATS) add_subdirectory(stats) diff --git a/test/spatial/CMakeLists.txt b/test/spatial/CMakeLists.txt new file mode 100644 index 000000000..9dc7aed59 --- /dev/null +++ b/test/spatial/CMakeLists.txt @@ -0,0 +1,7 @@ +set( + fppFiles + "test_spatial_kabsch.fypp" +) +fypp_f90pp("${fyppFlags}" "${fppFiles}" outFiles) + +ADDTESTPP(spatial_kabsch) \ No newline at end of file diff --git a/test/spatial/test_spatial_kabsch.fypp b/test/spatial/test_spatial_kabsch.fypp new file mode 100644 index 000000000..2f78bbbc4 --- /dev/null +++ b/test/spatial/test_spatial_kabsch.fypp @@ -0,0 +1,187 @@ +#:include "common.fypp" +#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) +#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) +module test_kabsch + use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test + use stdlib_kinds + use stdlib_math, only: all_close, is_close + use stdlib_linalg, only: svd + use stdlib_spatial + implicit none + +contains + + + !> Collect all exported unit tests + subroutine collect_suite(testsuite) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: testsuite(:) + + testsuite = [ & + new_unittest('kabsch_real', test_kabsch_real), & + new_unittest('kabsch_complex', test_kabsch_complex) & + ] + end subroutine + + subroutine test_kabsch_real(error) + type(error_type), allocatable, intent(out) :: error + #:for k, t, s in (R_KINDS_TYPES) + block + integer, parameter :: d = 3, N = 4 + ${t}$ :: P_original(d, N), Q_original(d, N), P_recovered(d, N) + ${t}$ :: R_recovered(d, d), R_original(d, d) + ${t}$ :: t_recovered(d), t_original(d) + ${t}$ :: U(d,d), Vt(d,d), S(d) + ${t}$ :: c_recovered, c_original + ${t}$ :: rmsd + + integer :: i, j + ${t}$ :: r1 + + call random_seed() + + ! Random reference points Q + call random_number(Q_original) + + ! Random proper rotation matrix R_original + ! Constructed via SVD: R = U * V^T + call random_number(R_original) + call svd(R_original, S, U, Vt) + R_original = matmul(U, Vt) + + ! Random scale and translation + call random_number(c_original) + call random_number(t_original) + + ! Construct P = c*R*Q + t + do j = 1, N + P_original(:,j) = c_original * matmul(R_original, Q_original(:,j)) + t_original + end do + + ! Call Kabsch–Umeyama + call kabsch(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) + + ! Apply recovered transform + do j = 1, N + P_recovered(:,j) = c_recovered * matmul(R_recovered, Q_original(:,j)) + t_recovered + end do + + call check(error, all_close(P_recovered, P_original), .true.) + if (allocated(error)) return + call check(error, all_close(R_recovered, R_original), .true.) + if (allocated(error)) return + call check(error, is_close(c_recovered, c_original), .true.) + if (allocated(error)) return + call check(error, all_close(t_recovered, t_original), .true.) + if (allocated(error)) return + end block + #:endfor + end subroutine + + subroutine test_kabsch_complex(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + #:for k, t, s in (C_KINDS_TYPES) + block + integer, parameter :: d = 3, N = 4 + ${t}$ :: P_original(d, N), Q_original(d, N), P_recovered(d, N) + ${t}$ :: R_recovered(d, d), R_original(d, d) + ${t}$ :: t_recovered(d), t_original(d) + ${t}$ :: c_recovered, c_original + ${t}$ :: U(d,d), Vt(d,d) + real(${k}$) :: S(d) + real(${k}$) :: rmsd + + integer :: i, j + real(${k}$) :: r1, r2 + + call random_seed() + + ! Random complex reference points Q + do j = 1, N + do i = 1, d + call random_number(r1) + call random_number(r2) + Q_original(i,j) = cmplx(r1, r2, kind=${k}$) + end do + end do + + ! Random complex unitary matrix R_original + ! Constructed via SVD: R = U * V^H + do i = 1, d + do j = 1, d + call random_number(r1) + call random_number(r2) + R_original(i,j) = cmplx(r1, r2, kind=${k}$) + end do + end do + call svd(R_original, S, U, Vt) + R_original = matmul(U, Vt) + + ! Random scale and translation + call random_number(r1) + call random_number(r2) + c_original = cmplx(r1,r2, kind=${k}$) + + do i = 1, d + call random_number(r1) + call random_number(r2) + t_original(i) = cmplx(r1, r2, kind=${k}$) + end do + + ! Construct P = c*R*Q + t + do j = 1, N + P_original(:,j) = c_original * matmul(R_original, Q_original(:,j)) + t_original + end do + + ! Call Kabsch–Umeyama + call kabsch(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) + + ! Apply recovered transform + do j = 1, N + P_recovered(:,j) = c_recovered * matmul(R_recovered, Q_original(:,j)) + t_recovered + end do + call check(error, all_close(P_recovered, P_original), .true.) + if (allocated(error)) print*, "return after P: " + if (allocated(error)) print*, "rmsd: ", rmsd + if (allocated(error)) return + call check(error, all_close(c_recovered * R_recovered, c_original * R_original), .true.) + if (allocated(error)) print*, "return after cR: " + if (allocated(error)) print*, "rmsd: ", rmsd + if (allocated(error)) return + call check(error, all_close(t_recovered, t_original), .true.) + if (allocated(error)) print*, "return after t: " + if (allocated(error)) print*, "rmsd: ", rmsd + if (allocated(error)) return + end block + #:endfor + end subroutine + +end module + + +program tester + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_kabsch, only : collect_suite + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("kabsch", collect_suite) & + ] + + do is = 1, size(testsuites) + write(error_unit, fmt) "Testing:", testsuites(is)%name + call run_testsuite(testsuites(is)%collect, error_unit, stat) + end do + + if (stat > 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program From 33859b604aa67961ecc3f282b304bdbee744b7c7 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Sat, 7 Feb 2026 17:25:44 +0530 Subject: [PATCH 07/29] remove: debug comments --- test/spatial/test_spatial_kabsch.fypp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/test/spatial/test_spatial_kabsch.fypp b/test/spatial/test_spatial_kabsch.fypp index 2f78bbbc4..33e07b2f5 100644 --- a/test/spatial/test_spatial_kabsch.fypp +++ b/test/spatial/test_spatial_kabsch.fypp @@ -142,16 +142,10 @@ contains P_recovered(:,j) = c_recovered * matmul(R_recovered, Q_original(:,j)) + t_recovered end do call check(error, all_close(P_recovered, P_original), .true.) - if (allocated(error)) print*, "return after P: " - if (allocated(error)) print*, "rmsd: ", rmsd if (allocated(error)) return call check(error, all_close(c_recovered * R_recovered, c_original * R_original), .true.) - if (allocated(error)) print*, "return after cR: " - if (allocated(error)) print*, "rmsd: ", rmsd if (allocated(error)) return call check(error, all_close(t_recovered, t_original), .true.) - if (allocated(error)) print*, "return after t: " - if (allocated(error)) print*, "rmsd: ", rmsd if (allocated(error)) return end block #:endfor From c7b3227887aaab71b57c6115b1dbede597154fc3 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Sat, 7 Feb 2026 18:04:21 +0530 Subject: [PATCH 08/29] remove: unnecessary fypp conditions --- src/spatial/stdlib_spatial_kabsch.fypp | 25 ------------------------- 1 file changed, 25 deletions(-) diff --git a/src/spatial/stdlib_spatial_kabsch.fypp b/src/spatial/stdlib_spatial_kabsch.fypp index 4e40788fb..1d014a453 100644 --- a/src/spatial/stdlib_spatial_kabsch.fypp +++ b/src/spatial/stdlib_spatial_kabsch.fypp @@ -33,9 +33,6 @@ contains real(${k}$) :: sum_w, variance_p real(${k}$), allocatable :: S(:) logical :: scale_ - #:if t.startswith('complex') - ${t}$ :: detR - #:endif ! Dimension checks if(size(P,dim=1)/=size(Q,dim=1) .or. size(P,dim=1)/=size(R,dim=1) .or. size(P,dim=1)/=size(R,dim=2) & @@ -58,13 +55,8 @@ contains sum_w = one_${s}$ / N if(present(W)) sum_w = one_${s}$ / stdlib_sum(W) - #:if t.startswith('complex') allocate(c_P(d), source=zero_${s}$) allocate(c_Q(d), source=zero_${s}$) - #:else - allocate(c_P(d), source=zero_${s}$) - allocate(c_Q(d), source=zero_${s}$) - #:endif ! Compute centroids of P and Q if(present(W)) then @@ -82,11 +74,7 @@ contains c_Q = c_Q * sum_w ! Compute covariance matrix H = (P - c_P) * (Q - c_Q)^T and variance of P - #:if t.startswith('complex') allocate(covariance(d,d), source=zero_${s}$) - #:else - allocate(covariance(d,d), source=zero_${s}$) - #:endif variance_p = zero_${s}$ if (present(W)) then @@ -138,23 +126,14 @@ contains covariance = covariance * sum_w variance_p = variance_p * sum_w - #:if t.startswith('complex') allocate(U(d,d), source=zero_${s}$) allocate(Vt(d,d), source=zero_${s}$) - #:else - allocate(U(d,d), source=zero_${s}$) - allocate(Vt(d,d), source=zero_${s}$) - #:endif allocate(S(d), source=zero_${k}$) ! SVD of covariance matrix H -> H = U * S * Vt call svd(covariance, S, U, Vt) - #:if t.startswith('complex') - allocate(B(d,d), source=zero_${s}$) - #:else allocate(B(d,d), source=zero_${s}$) - #:endif ! Optimal rotation matrix. R = matmul(U, Vt) @@ -167,11 +146,7 @@ contains t = c_P - c*matmul(R , c_Q ) ! Compute RMSD - #:if t.startswith('complex') - allocate(vec(d), source=zero_${s}$) - #:else allocate(vec(d), source=zero_${s}$) - #:endif rmsd = zero_${s}$ do i = 1, N vec(1:d) = c * matmul(R , Q(1:d,i)) From 78f35a4fe8fab1ab6dd5ffadf4468c1516bb26e9 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Mon, 9 Feb 2026 00:18:15 +0530 Subject: [PATCH 09/29] added instrinsics functions to spatial --- example/intrinsics/example_sum.f90 | 27 ++++++ src/spatial/stdlib_spatial.fypp | 2 +- src/spatial/stdlib_spatial_kabsch.fypp | 111 ++++++++++++++----------- test/spatial/test_spatial_kabsch.fypp | 58 ++++++++----- 4 files changed, 127 insertions(+), 71 deletions(-) diff --git a/example/intrinsics/example_sum.f90 b/example/intrinsics/example_sum.f90 index 0aec4835c..9ce7581ba 100644 --- a/example/intrinsics/example_sum.f90 +++ b/example/intrinsics/example_sum.f90 @@ -5,6 +5,8 @@ program example_sum real(sp), allocatable :: x(:) real(sp) :: total_sum(3) + real :: P(2,3), Q(2,3), res(2,2), tmp(3) + integer :: i, j allocate( x(1000) ) call random_number(x) @@ -14,4 +16,29 @@ program example_sum total_sum(3) = stdlib_sum_kahan(x)!> chunked kahan summation print *, total_sum(1:3) + print*, "my programme: " + P(:,1) = [1,2] + P(:,2) = [-1,3] + P(:,3) = [4,2] + Q(:,1) = [3,0] + Q(:,2) = [5,4] + Q(:,3) = [0,-1] + res = matmul(P, transpose(Q)) + print*, "" + do i = 1, 2 + print*, res(i, :) + end do + + + do j = 1, 2 + do i = 1, 2 + tmp(:) = P(i,:) * Q(j,:) + res(i,j) = stdlib_sum_kahan(tmp) + end do + end do + print*, "" + do i = 1, 2 + print*, res(i, :) + end do + end program example_sum \ No newline at end of file diff --git a/src/spatial/stdlib_spatial.fypp b/src/spatial/stdlib_spatial.fypp index f0ee0f197..ad3281687 100644 --- a/src/spatial/stdlib_spatial.fypp +++ b/src/spatial/stdlib_spatial.fypp @@ -39,7 +39,7 @@ module stdlib_spatial !> Root-mean-square deviation real(${k}$), intent(out) :: rmsd !> Optional weights - real(${k}$), intent(in), optional :: W(:) + ${t}$, intent(in), optional :: W(:) !> Enable scaling logical, intent(in), optional :: scale end subroutine diff --git a/src/spatial/stdlib_spatial_kabsch.fypp b/src/spatial/stdlib_spatial_kabsch.fypp index 1d014a453..57b6f1dde 100644 --- a/src/spatial/stdlib_spatial_kabsch.fypp +++ b/src/spatial/stdlib_spatial_kabsch.fypp @@ -4,7 +4,7 @@ #:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES submodule(stdlib_spatial) stdlib_spatial_kabsch use stdlib_linalg, only: svd, det, eye - use stdlib_intrinsics, only: stdlib_sum + use stdlib_intrinsics, only: stdlib_sum_kahan, stdlib_dot_product_kahan, stdlib_matmul, stdlib_sum, kahan_kernel contains #:for k, t, s in (KINDS_TYPES) @@ -22,17 +22,19 @@ contains !> Root-mean-square deviation real(${k}$), intent(out) :: rmsd !> Optional weights - real(${k}$), intent(in), optional :: W(:) + ${t}$, intent(in), optional :: W(:) !> Enable scaling logical, intent(in), optional :: scale ! Internal variables. integer(ilp) :: i, j, point, d, N - ${t}$, allocatable :: c_P(:), c_Q(:), covariance(:,:), U(:,:), Vt(:,:), B(:,:), vec(:) + ${t}$, allocatable :: covariance(:,:), U(:,:), Vt(:,:), B(:,:), vec(:), tmp_N(:), tmp_d(:), c_P(:), c_Q(:) ${t}$ :: vp, vq real(${k}$) :: sum_w, variance_p real(${k}$), allocatable :: S(:) logical :: scale_ + real(${k}$) :: rmsd_err + ! Dimension checks if(size(P,dim=1)/=size(Q,dim=1) .or. size(P,dim=1)/=size(R,dim=1) .or. size(P,dim=1)/=size(R,dim=2) & @@ -53,21 +55,21 @@ contains if(present(scale)) scale_ = scale sum_w = one_${s}$ / N - if(present(W)) sum_w = one_${s}$ / stdlib_sum(W) + if(present(W)) sum_w = one_${s}$ / stdlib_sum_kahan(W) allocate(c_P(d), source=zero_${s}$) allocate(c_Q(d), source=zero_${s}$) ! Compute centroids of P and Q if(present(W)) then - do point = 1, N - c_P(1:d) = c_P(1:d) + P(1:d,point)*W(point) - c_Q(1:d) = c_Q(1:d) + Q(1:d,point)*W(point) + do i = 1, d + c_P(i) = stdlib_dot_product_kahan(w,P(i, :)) + c_Q(i) = stdlib_dot_product_kahan(w,Q(i, :)) end do else - do point = 1, N - c_P(1:d) = c_P(1:d) + P(1:d,point) - c_Q(1:d) = c_Q(1:d) + Q(1:d,point) + do i = 1, d + c_P(i) = stdlib_sum_kahan(P(i, :)) + c_Q(i) = stdlib_sum_kahan(Q(i, :)) end do end if c_P = c_P * sum_w @@ -75,51 +77,43 @@ contains ! Compute covariance matrix H = (P - c_P) * (Q - c_Q)^T and variance of P allocate(covariance(d,d), source=zero_${s}$) + allocate(tmp_N(N), source=zero_${s}$) + allocate(tmp_d(d), source=zero_${s}$) variance_p = zero_${s}$ if (present(W)) then do point = 1, N + tmp_d = P(:, point) - c_P(:) + tmp_N(point) = stdlib_dot_product_kahan(tmp_d, tmp_d) + end do + variance_p = stdlib_dot_product_kahan(w, tmp_N) + do j = 1, d do i = 1, d - vp = P(i,point) - c_P(i) #:if t.startswith('complex') - variance_p = variance_p + conjg(vp)*vp * W(point) + tmp_N(:) = (P(i,:) - c_P(i)) * conjg(Q(j,:) - c_Q(j)) + covariance(i,j) = stdlib_dot_product_kahan(w, tmp_N) #:else - variance_p = variance_p + vp*vp * W(point) + tmp_N(:) = (P(i,:) - c_P(i)) * (Q(j,:) - c_Q(j)) + covariance(i,j) = stdlib_dot_product_kahan(w, tmp_N) #:endif end do - do j = 1, d - vq = Q(j,point) - c_Q(j) - do i = 1, d - vp = P(i,point) - c_P(i) - #:if t.startswith('complex') - covariance(i,j) = covariance(i,j) + vp*conjg(vq) * W(point) - #:else - covariance(i,j) = covariance(i,j) + vp*vq * W(point) - #:endif - end do - end do end do else do point = 1, N + tmp_d = P(:, point) - c_P(:) + tmp_N(point) = stdlib_dot_product_kahan(tmp_d, tmp_d) + end do + variance_p = stdlib_sum_kahan(tmp_N) + do j = 1, d do i = 1, d - vp = P(i,point) - c_P(i) #:if t.startswith('complex') - variance_p = variance_p + conjg(vp)*vp + tmp_N(:) = (P(i,:) - c_P(i)) * conjg(Q(j,:) - c_Q(j)) + covariance(i,j) = stdlib_sum_kahan(tmp_N) #:else - variance_p = variance_p + vp*vp + tmp_N(:) = (P(i,:) - c_P(i)) * (Q(j,:) - c_Q(j)) + covariance(i,j) = stdlib_sum_kahan(tmp_N) #:endif end do - do j = 1, d - vq = Q(j,point) - c_Q(j) - do i = 1, d - vp = P(i,point) - c_P(i) - #:if t.startswith('complex') - covariance(i,j) = covariance(i,j) + vp*conjg(vq) - #:else - covariance(i,j) = covariance(i,j) + vp*vq - #:endif - end do - end do end do end if @@ -136,27 +130,44 @@ contains allocate(B(d,d), source=zero_${s}$) ! Optimal rotation matrix. - R = matmul(U, Vt) + do i = 1,d + do j = 1,d + #:if t.startswith('complex') + R(i,j) = stdlib_dot_product_kahan(conjg(U(i,:)), Vt(:, j)) + #:else + R(i,j) = stdlib_dot_product_kahan(U(i,:), Vt(:, j)) + #:endif + end do + end do ! Scaling factor c = variance_p / (sum(S(1:d))) if (.not. scale_) c = one_${s}$ - ! Translation vector - t = c_P - c*matmul(R , c_Q ) + ! Translation vector t = c_P - c*R*c_Q + do i = 1, d + #:if t.startswith('complex') + t(i) = c_P(i) - c * stdlib_dot_product_kahan(conjg(R(i,1:d)), c_Q(1:d)) + #:else + t(i) = c_P(i) - c * stdlib_dot_product_kahan(R(i,1:d), c_Q(1:d)) + #:endif + end do ! Compute RMSD allocate(vec(d), source=zero_${s}$) rmsd = zero_${s}$ - do i = 1, N - vec(1:d) = c * matmul(R , Q(1:d,i)) - vec(1:d) = vec(1:d) + t(1:d) - vec(1:d) = vec(1:d) - P(1:d,i) - #:if t.startswith('complex') - rmsd = rmsd + real(dot_product(vec, vec), kind=${k}$) - #:else - rmsd = rmsd + dot_product(vec, vec) - #:endif + rmsd_err = zero_${s}$ + do point = 1, N + ! Calculate the k^th difference vector by the formula vec_k = c*R*Q_k + t - P_k + do i = 1, d + #:if t.startswith('complex') + vec(i) = c * stdlib_dot_product_kahan(conjg(R(i,1:d)), Q(1:d,point)) + #:else + vec(i) = c * stdlib_dot_product_kahan(R(i,1:d), Q(1:d,point)) + #:endif + end do + vec(1:d) = vec(1:d) + t(1:d) - P(1:d,point) + call kahan_kernel(real(stdlib_dot_product_kahan(vec,vec), kind=${k}$), rmsd, rmsd_err) end do rmsd = sqrt(rmsd * sum_w) end subroutine diff --git a/test/spatial/test_spatial_kabsch.fypp b/test/spatial/test_spatial_kabsch.fypp index 33e07b2f5..2d5ae6dbe 100644 --- a/test/spatial/test_spatial_kabsch.fypp +++ b/test/spatial/test_spatial_kabsch.fypp @@ -7,6 +7,7 @@ module test_kabsch use stdlib_math, only: all_close, is_close use stdlib_linalg, only: svd use stdlib_spatial + use stdlib_intrinsics, only: stdlib_sum_kahan, stdlib_dot_product_kahan, stdlib_matmul, stdlib_sum, kahan_kernel implicit none contains @@ -28,7 +29,7 @@ contains #:for k, t, s in (R_KINDS_TYPES) block integer, parameter :: d = 3, N = 4 - ${t}$ :: P_original(d, N), Q_original(d, N), P_recovered(d, N) + ${t}$ :: P_original(d, N), Q_original(d, N) ${t}$ :: R_recovered(d, d), R_original(d, d) ${t}$ :: t_recovered(d), t_original(d) ${t}$ :: U(d,d), Vt(d,d), S(d) @@ -47,7 +48,11 @@ contains ! Constructed via SVD: R = U * V^T call random_number(R_original) call svd(R_original, S, U, Vt) - R_original = matmul(U, Vt) + do i = 1,d + do j = 1,d + R_original(i,j) = stdlib_dot_product_kahan(U(i,:), Vt(:, j)) + end do + end do ! Random scale and translation call random_number(c_original) @@ -55,24 +60,28 @@ contains ! Construct P = c*R*Q + t do j = 1, N - P_original(:,j) = c_original * matmul(R_original, Q_original(:,j)) + t_original + do i = 1, d + P_original(i,j) = c_original * & + stdlib_dot_product_kahan(R_original(i,:), Q_original(:,j)) + & + t_original(i) + end do end do ! Call Kabsch–Umeyama call kabsch(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) - ! Apply recovered transform - do j = 1, N - P_recovered(:,j) = c_recovered * matmul(R_recovered, Q_original(:,j)) + t_recovered - end do - - call check(error, all_close(P_recovered, P_original), .true.) - if (allocated(error)) return call check(error, all_close(R_recovered, R_original), .true.) + if (allocated(error)) then + print*, "R did not match in real" + print*, "R_recovered: ", R_recovered + print*, "R_original: ", R_original + end if if (allocated(error)) return call check(error, is_close(c_recovered, c_original), .true.) + if (allocated(error)) print*, "c did not match in real", ", recovered: ", c_recovered, "original: ", c_original if (allocated(error)) return call check(error, all_close(t_recovered, t_original), .true.) + if (allocated(error)) print*, "t did not match in real", ", recovered: ", t_recovered, "original: ", t_original if (allocated(error)) return end block #:endfor @@ -84,8 +93,8 @@ contains #:for k, t, s in (C_KINDS_TYPES) block integer, parameter :: d = 3, N = 4 - ${t}$ :: P_original(d, N), Q_original(d, N), P_recovered(d, N) - ${t}$ :: R_recovered(d, d), R_original(d, d) + ${t}$ :: P_original(d, N), Q_original(d, N) + ${t}$ :: R_recovered(d, d), R_original(d, d), cR_original(d, d), cR_recovered(d, d) ${t}$ :: t_recovered(d), t_original(d) ${t}$ :: c_recovered, c_original ${t}$ :: U(d,d), Vt(d,d) @@ -116,7 +125,11 @@ contains end do end do call svd(R_original, S, U, Vt) - R_original = matmul(U, Vt) + do i = 1,d + do j = 1,d + R_original(i,j) = stdlib_dot_product_kahan(conjg(U(i,:)), Vt(:, j)) + end do + end do ! Random scale and translation call random_number(r1) @@ -131,21 +144,26 @@ contains ! Construct P = c*R*Q + t do j = 1, N - P_original(:,j) = c_original * matmul(R_original, Q_original(:,j)) + t_original + do i = 1, d + P_original(i,j) = c_original * & + stdlib_dot_product_kahan(conjg(R_original(i,:)), Q_original(:,j)) + & + t_original(i) + end do end do ! Call Kabsch–Umeyama call kabsch(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) - ! Apply recovered transform - do j = 1, N - P_recovered(:,j) = c_recovered * matmul(R_recovered, Q_original(:,j)) + t_recovered - end do - call check(error, all_close(P_recovered, P_original), .true.) - if (allocated(error)) return + print*, "RMSD: ",rmsd call check(error, all_close(c_recovered * R_recovered, c_original * R_original), .true.) + if (allocated(error)) then + print*, "cR did not match in real" + print*, "cR_recovered: ", c_recovered * R_recovered + print*, "cR_original: ", c_original * R_original + end if if (allocated(error)) return call check(error, all_close(t_recovered, t_original), .true.) + if (allocated(error)) print*, "t did not match in complex", ", recovered: ", t_recovered, "original: ", t_original if (allocated(error)) return end block #:endfor From a9c489cf3040ad44feeb9c868c7793690416a983 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Wed, 11 Feb 2026 11:22:26 +0530 Subject: [PATCH 10/29] modify: test only for P --- src/spatial/stdlib_spatial_kabsch.fypp | 2 - test/spatial/test_spatial_kabsch.fypp | 73 ++++++++++++++------------ 2 files changed, 39 insertions(+), 36 deletions(-) diff --git a/src/spatial/stdlib_spatial_kabsch.fypp b/src/spatial/stdlib_spatial_kabsch.fypp index 57b6f1dde..deab82df7 100644 --- a/src/spatial/stdlib_spatial_kabsch.fypp +++ b/src/spatial/stdlib_spatial_kabsch.fypp @@ -127,8 +127,6 @@ contains ! SVD of covariance matrix H -> H = U * S * Vt call svd(covariance, S, U, Vt) - allocate(B(d,d), source=zero_${s}$) - ! Optimal rotation matrix. do i = 1,d do j = 1,d diff --git a/test/spatial/test_spatial_kabsch.fypp b/test/spatial/test_spatial_kabsch.fypp index 2d5ae6dbe..6ea0dcd39 100644 --- a/test/spatial/test_spatial_kabsch.fypp +++ b/test/spatial/test_spatial_kabsch.fypp @@ -10,6 +10,10 @@ module test_kabsch use stdlib_intrinsics, only: stdlib_sum_kahan, stdlib_dot_product_kahan, stdlib_matmul, stdlib_sum, kahan_kernel implicit none + #:for k, t, s in (R_KINDS_TYPES) + real(${k}$), parameter :: ${s}$tol = 100 * epsilon(1._${k}$) + #:endfor + contains @@ -28,24 +32,21 @@ contains type(error_type), allocatable, intent(out) :: error #:for k, t, s in (R_KINDS_TYPES) block - integer, parameter :: d = 3, N = 4 - ${t}$ :: P_original(d, N), Q_original(d, N) + integer, parameter :: d = 3, N = 8 + ${t}$ :: P_original(d, N), P_recovered(d, N), Q_original(d, N) ${t}$ :: R_recovered(d, d), R_original(d, d) ${t}$ :: t_recovered(d), t_original(d) - ${t}$ :: U(d,d), Vt(d,d), S(d) ${t}$ :: c_recovered, c_original + ${t}$ :: U(d,d), Vt(d,d), S(d) ${t}$ :: rmsd integer :: i, j ${t}$ :: r1 - call random_seed() - ! Random reference points Q call random_number(Q_original) - ! Random proper rotation matrix R_original - ! Constructed via SVD: R = U * V^T + ! Random proper rotation matrix R_original constructed via SVD: R = U * V^T call random_number(R_original) call svd(R_original, S, U, Vt) do i = 1,d @@ -70,19 +71,21 @@ contains ! Call Kabsch–Umeyama call kabsch(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) - call check(error, all_close(R_recovered, R_original), .true.) - if (allocated(error)) then - print*, "R did not match in real" - print*, "R_recovered: ", R_recovered - print*, "R_original: ", R_original - end if - if (allocated(error)) return - call check(error, is_close(c_recovered, c_original), .true.) - if (allocated(error)) print*, "c did not match in real", ", recovered: ", c_recovered, "original: ", c_original - if (allocated(error)) return - call check(error, all_close(t_recovered, t_original), .true.) - if (allocated(error)) print*, "t did not match in real", ", recovered: ", t_recovered, "original: ", t_original + ! Apply recovered transform + do j = 1, N + do i = 1, d + P_recovered(i,j) = c_recovered * & + stdlib_dot_product_kahan(R_recovered(i,:), Q_original(:,j)) + & + t_recovered(i) + end do + end do + + call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) + if (allocated(error)) print*, "P did not match in real", ", recovered: ", P_recovered, "original: ", P_original if (allocated(error)) return + call check(error, is_close(rmsd, 0.0_${k}$, abs_tol = ${k}$tol), .true.) + if(allocated(error)) print*, "RMSD: ",rmsd + if(allocated(error)) return end block #:endfor end subroutine @@ -92,9 +95,9 @@ contains type(error_type), allocatable, intent(out) :: error #:for k, t, s in (C_KINDS_TYPES) block - integer, parameter :: d = 3, N = 4 - ${t}$ :: P_original(d, N), Q_original(d, N) - ${t}$ :: R_recovered(d, d), R_original(d, d), cR_original(d, d), cR_recovered(d, d) + integer, parameter :: d = 3, N = 8 + ${t}$ :: P_original(d, N), Q_original(d, N), P_recovered(d, N) + ${t}$ :: R_recovered(d, d), R_original(d, d) ${t}$ :: t_recovered(d), t_original(d) ${t}$ :: c_recovered, c_original ${t}$ :: U(d,d), Vt(d,d) @@ -104,8 +107,6 @@ contains integer :: i, j real(${k}$) :: r1, r2 - call random_seed() - ! Random complex reference points Q do j = 1, N do i = 1, d @@ -154,17 +155,21 @@ contains ! Call Kabsch–Umeyama call kabsch(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) - print*, "RMSD: ",rmsd - call check(error, all_close(c_recovered * R_recovered, c_original * R_original), .true.) - if (allocated(error)) then - print*, "cR did not match in real" - print*, "cR_recovered: ", c_recovered * R_recovered - print*, "cR_original: ", c_original * R_original - end if - if (allocated(error)) return - call check(error, all_close(t_recovered, t_original), .true.) - if (allocated(error)) print*, "t did not match in complex", ", recovered: ", t_recovered, "original: ", t_original + ! Apply recovered transform + do j = 1, N + do i = 1, d + P_recovered(i,j) = c_recovered * & + stdlib_dot_product_kahan(conjg(R_recovered(i,:)), Q_original(:,j)) + & + t_recovered(i) + end do + end do + + call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) + if (allocated(error)) print*, "P did not match in complex", ", recovered: ", P_recovered, "original: ", P_original if (allocated(error)) return + call check(error, is_close(rmsd, 0.0_${k}$, abs_tol = ${k}$tol), .true.) + if(allocated(error)) print*, "RMSD: ",rmsd + if(allocated(error)) return end block #:endfor end subroutine From fb08933043a3f47e93e14bd0bf5719709acd8e77 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Wed, 11 Feb 2026 12:26:44 +0530 Subject: [PATCH 11/29] modify: set arguments to pass to kabsh in test to zero --- test/spatial/test_spatial_kabsch.fypp | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/test/spatial/test_spatial_kabsch.fypp b/test/spatial/test_spatial_kabsch.fypp index 6ea0dcd39..2376f79b3 100644 --- a/test/spatial/test_spatial_kabsch.fypp +++ b/test/spatial/test_spatial_kabsch.fypp @@ -43,6 +43,14 @@ contains integer :: i, j ${t}$ :: r1 + R_recovered = 0.0_${k}$ + t_recovered = 0.0_${k}$ + c_recovered = 0.0_${k}$ + U = 0.0_${k}$ + S = 0.0_${k}$ + Vt = 0.0_${k}$ + rmsd = 0.0_${k}$ + ! Random reference points Q call random_number(Q_original) @@ -107,6 +115,14 @@ contains integer :: i, j real(${k}$) :: r1, r2 + R_recovered = (0.0_${k}$, 0.0_${k}$) + t_recovered = (0.0_${k}$, 0.0_${k}$) + c_recovered = (0.0_${k}$, 0.0_${k}$) + U = (0.0_${k}$, 0.0_${k}$) + S = (0.0_${k}$, 0.0_${k}$) + Vt = (0.0_${k}$, 0.0_${k}$) + rmsd = (0.0_${k}$, 0.0_${k}$) + ! Random complex reference points Q do j = 1, N do i = 1, d From 88e76fba5063902a6d1535011aedfa7387bb26ee Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Wed, 11 Feb 2026 12:35:01 +0530 Subject: [PATCH 12/29] remove: debug statements in test --- test/spatial/test_spatial_kabsch.fypp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/test/spatial/test_spatial_kabsch.fypp b/test/spatial/test_spatial_kabsch.fypp index 2376f79b3..7c1d1bddd 100644 --- a/test/spatial/test_spatial_kabsch.fypp +++ b/test/spatial/test_spatial_kabsch.fypp @@ -89,10 +89,8 @@ contains end do call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) - if (allocated(error)) print*, "P did not match in real", ", recovered: ", P_recovered, "original: ", P_original if (allocated(error)) return call check(error, is_close(rmsd, 0.0_${k}$, abs_tol = ${k}$tol), .true.) - if(allocated(error)) print*, "RMSD: ",rmsd if(allocated(error)) return end block #:endfor @@ -181,10 +179,8 @@ contains end do call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) - if (allocated(error)) print*, "P did not match in complex", ", recovered: ", P_recovered, "original: ", P_original if (allocated(error)) return call check(error, is_close(rmsd, 0.0_${k}$, abs_tol = ${k}$tol), .true.) - if(allocated(error)) print*, "RMSD: ",rmsd if(allocated(error)) return end block #:endfor From e73a93f6d6230b88bb84b9086d5a0021a3946bfb Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Wed, 11 Feb 2026 23:15:28 +0530 Subject: [PATCH 13/29] remove: unused imports and modify: S and rmsd init inside complex kabsch test --- src/spatial/stdlib_spatial_kabsch.fypp | 4 ++-- test/spatial/test_spatial_kabsch.fypp | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/spatial/stdlib_spatial_kabsch.fypp b/src/spatial/stdlib_spatial_kabsch.fypp index deab82df7..64afbf2e0 100644 --- a/src/spatial/stdlib_spatial_kabsch.fypp +++ b/src/spatial/stdlib_spatial_kabsch.fypp @@ -3,8 +3,8 @@ #:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) #:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES submodule(stdlib_spatial) stdlib_spatial_kabsch - use stdlib_linalg, only: svd, det, eye - use stdlib_intrinsics, only: stdlib_sum_kahan, stdlib_dot_product_kahan, stdlib_matmul, stdlib_sum, kahan_kernel + use stdlib_linalg, only: svd + use stdlib_intrinsics, only: stdlib_sum_kahan, stdlib_dot_product_kahan, kahan_kernel contains #:for k, t, s in (KINDS_TYPES) diff --git a/test/spatial/test_spatial_kabsch.fypp b/test/spatial/test_spatial_kabsch.fypp index 7c1d1bddd..5a77d6675 100644 --- a/test/spatial/test_spatial_kabsch.fypp +++ b/test/spatial/test_spatial_kabsch.fypp @@ -7,7 +7,7 @@ module test_kabsch use stdlib_math, only: all_close, is_close use stdlib_linalg, only: svd use stdlib_spatial - use stdlib_intrinsics, only: stdlib_sum_kahan, stdlib_dot_product_kahan, stdlib_matmul, stdlib_sum, kahan_kernel + use stdlib_intrinsics, only: stdlib_dot_product_kahan implicit none #:for k, t, s in (R_KINDS_TYPES) @@ -117,9 +117,9 @@ contains t_recovered = (0.0_${k}$, 0.0_${k}$) c_recovered = (0.0_${k}$, 0.0_${k}$) U = (0.0_${k}$, 0.0_${k}$) - S = (0.0_${k}$, 0.0_${k}$) + S = 0.0_${k}$ Vt = (0.0_${k}$, 0.0_${k}$) - rmsd = (0.0_${k}$, 0.0_${k}$) + rmsd = 0.0_${k}$ ! Random complex reference points Q do j = 1, N From e3a1cc2186b9a3020085a334a8737726bba38b5e Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Fri, 13 Feb 2026 23:06:35 +0530 Subject: [PATCH 14/29] change name: kabsch to kabsch_umeyama --- example/spatial/CMakeLists.txt | 2 +- ...e_kabsch.f90 => example_kabsch_umeyama.f90} | 8 ++++---- src/spatial/CMakeLists.txt | 2 +- src/spatial/stdlib_spatial.fypp | 6 +++--- ...fypp => stdlib_spatial_kabsch_umeyama.fypp} | 6 +++--- test/spatial/CMakeLists.txt | 4 ++-- ...h.fypp => test_spatial_kabsch_umeyama.fypp} | 18 +++++++++--------- 7 files changed, 23 insertions(+), 23 deletions(-) rename example/spatial/{example_kabsch.f90 => example_kabsch_umeyama.f90} (80%) rename src/spatial/{stdlib_spatial_kabsch.fypp => stdlib_spatial_kabsch_umeyama.fypp} (97%) rename test/spatial/{test_spatial_kabsch.fypp => test_spatial_kabsch_umeyama.fypp} (92%) diff --git a/example/spatial/CMakeLists.txt b/example/spatial/CMakeLists.txt index 54540eb25..867f35f11 100644 --- a/example/spatial/CMakeLists.txt +++ b/example/spatial/CMakeLists.txt @@ -1 +1 @@ -ADD_EXAMPLE(kabsch) \ No newline at end of file +ADD_EXAMPLE(kabsch_umeyama) \ No newline at end of file diff --git a/example/spatial/example_kabsch.f90 b/example/spatial/example_kabsch_umeyama.f90 similarity index 80% rename from example/spatial/example_kabsch.f90 rename to example/spatial/example_kabsch_umeyama.f90 index 4ae37d794..90c59b112 100644 --- a/example/spatial/example_kabsch.f90 +++ b/example/spatial/example_kabsch_umeyama.f90 @@ -1,6 +1,6 @@ -program example_kabsch +program example_kabsch_umeyama use stdlib_linalg_constants, only: dp - use stdlib_spatial, only: kabsch + use stdlib_spatial, only: kabsch_umeyama implicit none integer, parameter :: d = 2, N = 3 @@ -16,7 +16,7 @@ program example_kabsch Q(:,2) = [-1.0_dp, 5.0_dp] Q(:,3) = [1.0_dp, 4.0_dp] - call kabsch(P, Q, R, t, c, rmsd) + call kabsch_umeyama(P, Q, R, t, c, rmsd) print *, "" print *, "Recovered rotation R:" @@ -33,4 +33,4 @@ program example_kabsch print *, "" print *, "RMSD:", rmsd -end program example_kabsch \ No newline at end of file +end program example_kabsch_umeyama \ No newline at end of file diff --git a/src/spatial/CMakeLists.txt b/src/spatial/CMakeLists.txt index c78edd60e..58c671e7f 100644 --- a/src/spatial/CMakeLists.txt +++ b/src/spatial/CMakeLists.txt @@ -1,6 +1,6 @@ set(spatial_fppFiles stdlib_spatial.fypp - stdlib_spatial_kabsch.fypp + stdlib_spatial_kabsch_umeyama.fypp ) set(spatial_cppFiles diff --git a/src/spatial/stdlib_spatial.fypp b/src/spatial/stdlib_spatial.fypp index ad3281687..d12f6f383 100644 --- a/src/spatial/stdlib_spatial.fypp +++ b/src/spatial/stdlib_spatial.fypp @@ -8,9 +8,9 @@ module stdlib_spatial use stdlib_error, only: error_stop implicit none private - public :: kabsch + public :: kabsch_umeyama - interface kabsch + interface kabsch_umeyama !----------------------------------------------------------------------- !> Compute the optimal similarity transform (Kabsch–Umeyama): !> @@ -25,7 +25,7 @@ module stdlib_spatial !> of P and Q, optionally using weights. !----------------------------------------------------------------------- #:for k, t, s in (KINDS_TYPES) - module subroutine kabsch_${s}$(P, Q, R, t, c, rmsd, W, scale) + module subroutine kabsch_umeyama_${s}$(P, Q, R, t, c, rmsd, W, scale) !> Reference point set (d × N) ${t}$, intent(in) :: P(:, :) !> Target point set (d × N) diff --git a/src/spatial/stdlib_spatial_kabsch.fypp b/src/spatial/stdlib_spatial_kabsch_umeyama.fypp similarity index 97% rename from src/spatial/stdlib_spatial_kabsch.fypp rename to src/spatial/stdlib_spatial_kabsch_umeyama.fypp index 64afbf2e0..2c9f6dbba 100644 --- a/src/spatial/stdlib_spatial_kabsch.fypp +++ b/src/spatial/stdlib_spatial_kabsch_umeyama.fypp @@ -2,13 +2,13 @@ #:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) #:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) #:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES -submodule(stdlib_spatial) stdlib_spatial_kabsch +submodule(stdlib_spatial) stdlib_spatial_kabsch_umeyama use stdlib_linalg, only: svd use stdlib_intrinsics, only: stdlib_sum_kahan, stdlib_dot_product_kahan, kahan_kernel contains #:for k, t, s in (KINDS_TYPES) - module subroutine kabsch_${s}$(P, Q, R, t, c, rmsd, W, scale) + module subroutine kabsch_umeyama_${s}$(P, Q, R, t, c, rmsd, W, scale) !> Reference point set (d × N) ${t}$, intent(in) :: P(:, :) !> Target point set (d × N) @@ -170,4 +170,4 @@ contains rmsd = sqrt(rmsd * sum_w) end subroutine #:endfor -end submodule stdlib_spatial_kabsch \ No newline at end of file +end submodule stdlib_spatial_kabsch_umeyama \ No newline at end of file diff --git a/test/spatial/CMakeLists.txt b/test/spatial/CMakeLists.txt index 9dc7aed59..a64478b2d 100644 --- a/test/spatial/CMakeLists.txt +++ b/test/spatial/CMakeLists.txt @@ -1,7 +1,7 @@ set( fppFiles - "test_spatial_kabsch.fypp" + "test_spatial_kabsch_umeyama.fypp" ) fypp_f90pp("${fyppFlags}" "${fppFiles}" outFiles) -ADDTESTPP(spatial_kabsch) \ No newline at end of file +ADDTESTPP(spatial_kabsch_umeyama) \ No newline at end of file diff --git a/test/spatial/test_spatial_kabsch.fypp b/test/spatial/test_spatial_kabsch_umeyama.fypp similarity index 92% rename from test/spatial/test_spatial_kabsch.fypp rename to test/spatial/test_spatial_kabsch_umeyama.fypp index 5a77d6675..0a00dcc81 100644 --- a/test/spatial/test_spatial_kabsch.fypp +++ b/test/spatial/test_spatial_kabsch_umeyama.fypp @@ -1,7 +1,7 @@ #:include "common.fypp" #:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) #:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) -module test_kabsch +module test_kabsch_umeyama use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test use stdlib_kinds use stdlib_math, only: all_close, is_close @@ -23,12 +23,12 @@ contains type(unittest_type), allocatable, intent(out) :: testsuite(:) testsuite = [ & - new_unittest('kabsch_real', test_kabsch_real), & - new_unittest('kabsch_complex', test_kabsch_complex) & + new_unittest('kabsch_umeyama_real', test_kabsch_umeyama_real), & + new_unittest('kabsch_umeyama_complex', test_kabsch_umeyama_complex) & ] end subroutine - subroutine test_kabsch_real(error) + subroutine test_kabsch_umeyama_real(error) type(error_type), allocatable, intent(out) :: error #:for k, t, s in (R_KINDS_TYPES) block @@ -77,7 +77,7 @@ contains end do ! Call Kabsch–Umeyama - call kabsch(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) + call kabsch_umeyama(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) ! Apply recovered transform do j = 1, N @@ -96,7 +96,7 @@ contains #:endfor end subroutine - subroutine test_kabsch_complex(error) + subroutine test_kabsch_umeyama_complex(error) !> Error handling type(error_type), allocatable, intent(out) :: error #:for k, t, s in (C_KINDS_TYPES) @@ -167,7 +167,7 @@ contains end do ! Call Kabsch–Umeyama - call kabsch(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) + call kabsch_umeyama(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) ! Apply recovered transform do j = 1, N @@ -192,7 +192,7 @@ end module program tester use, intrinsic :: iso_fortran_env, only : error_unit use testdrive, only : run_testsuite, new_testsuite, testsuite_type - use test_kabsch, only : collect_suite + use test_kabsch_umeyama, only : collect_suite implicit none integer :: stat, is type(testsuite_type), allocatable :: testsuites(:) @@ -201,7 +201,7 @@ program tester stat = 0 testsuites = [ & - new_testsuite("kabsch", collect_suite) & + new_testsuite("kabsch_umeyama", collect_suite) & ] do is = 1, size(testsuites) From 27b68ca9fdb9f1813e1665f2a19758dca621b1de Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Fri, 13 Feb 2026 23:13:14 +0530 Subject: [PATCH 15/29] add stdlib_constants to test_kabsch_umeyama.fypp --- test/spatial/test_spatial_kabsch_umeyama.fypp | 35 ++++++++++--------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/test/spatial/test_spatial_kabsch_umeyama.fypp b/test/spatial/test_spatial_kabsch_umeyama.fypp index 0a00dcc81..9c532fc6b 100644 --- a/test/spatial/test_spatial_kabsch_umeyama.fypp +++ b/test/spatial/test_spatial_kabsch_umeyama.fypp @@ -8,10 +8,11 @@ module test_kabsch_umeyama use stdlib_linalg, only: svd use stdlib_spatial use stdlib_intrinsics, only: stdlib_dot_product_kahan + use stdlib_constants implicit none #:for k, t, s in (R_KINDS_TYPES) - real(${k}$), parameter :: ${s}$tol = 100 * epsilon(1._${k}$) + real(${k}$), parameter :: ${s}$tol = 100 * epsilon(one_${k}$) #:endfor contains @@ -43,13 +44,13 @@ contains integer :: i, j ${t}$ :: r1 - R_recovered = 0.0_${k}$ - t_recovered = 0.0_${k}$ - c_recovered = 0.0_${k}$ - U = 0.0_${k}$ - S = 0.0_${k}$ - Vt = 0.0_${k}$ - rmsd = 0.0_${k}$ + R_recovered = zero_${s}$ + t_recovered = zero_${s}$ + c_recovered = zero_${s}$ + U = zero_${s}$ + S = zero_${k}$ + Vt = zero_${s}$ + rmsd = zero_${k}$ ! Random reference points Q call random_number(Q_original) @@ -90,7 +91,7 @@ contains call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) if (allocated(error)) return - call check(error, is_close(rmsd, 0.0_${k}$, abs_tol = ${k}$tol), .true.) + call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) if(allocated(error)) return end block #:endfor @@ -113,13 +114,13 @@ contains integer :: i, j real(${k}$) :: r1, r2 - R_recovered = (0.0_${k}$, 0.0_${k}$) - t_recovered = (0.0_${k}$, 0.0_${k}$) - c_recovered = (0.0_${k}$, 0.0_${k}$) - U = (0.0_${k}$, 0.0_${k}$) - S = 0.0_${k}$ - Vt = (0.0_${k}$, 0.0_${k}$) - rmsd = 0.0_${k}$ + R_recovered = zero_${s}$ + t_recovered = zero_${s}$ + c_recovered = zero_${s}$ + U = zero_${s}$ + S = zero_${k}$ + Vt = zero_${s}$ + rmsd = zero_${k}$ ! Random complex reference points Q do j = 1, N @@ -180,7 +181,7 @@ contains call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) if (allocated(error)) return - call check(error, is_close(rmsd, 0.0_${k}$, abs_tol = ${k}$tol), .true.) + call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) if(allocated(error)) return end block #:endfor From c895ed56bc6c8f89cb28e130f91091b5ba6aaccf Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Fri, 13 Feb 2026 23:20:40 +0530 Subject: [PATCH 16/29] move tol declaration inside blocks --- test/spatial/test_spatial_kabsch_umeyama.fypp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/test/spatial/test_spatial_kabsch_umeyama.fypp b/test/spatial/test_spatial_kabsch_umeyama.fypp index 9c532fc6b..97b891a2b 100644 --- a/test/spatial/test_spatial_kabsch_umeyama.fypp +++ b/test/spatial/test_spatial_kabsch_umeyama.fypp @@ -11,10 +11,6 @@ module test_kabsch_umeyama use stdlib_constants implicit none - #:for k, t, s in (R_KINDS_TYPES) - real(${k}$), parameter :: ${s}$tol = 100 * epsilon(one_${k}$) - #:endfor - contains @@ -43,6 +39,7 @@ contains integer :: i, j ${t}$ :: r1 + real(${k}$), parameter :: ${k}$tol = 100 * epsilon(one_${k}$) R_recovered = zero_${s}$ t_recovered = zero_${s}$ @@ -113,6 +110,7 @@ contains integer :: i, j real(${k}$) :: r1, r2 + real(${k}$), parameter :: ${k}$tol = 100 * epsilon(one_${k}$) R_recovered = zero_${s}$ t_recovered = zero_${s}$ From d722366c8a4f5af75712a85e8117aaf0a0cb0a03 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Sat, 14 Feb 2026 00:20:32 +0530 Subject: [PATCH 17/29] revert accidental changes to example_sum.f90 --- example/intrinsics/example_sum.f90 | 27 --------------------------- 1 file changed, 27 deletions(-) diff --git a/example/intrinsics/example_sum.f90 b/example/intrinsics/example_sum.f90 index 9ce7581ba..0aec4835c 100644 --- a/example/intrinsics/example_sum.f90 +++ b/example/intrinsics/example_sum.f90 @@ -5,8 +5,6 @@ program example_sum real(sp), allocatable :: x(:) real(sp) :: total_sum(3) - real :: P(2,3), Q(2,3), res(2,2), tmp(3) - integer :: i, j allocate( x(1000) ) call random_number(x) @@ -16,29 +14,4 @@ program example_sum total_sum(3) = stdlib_sum_kahan(x)!> chunked kahan summation print *, total_sum(1:3) - print*, "my programme: " - P(:,1) = [1,2] - P(:,2) = [-1,3] - P(:,3) = [4,2] - Q(:,1) = [3,0] - Q(:,2) = [5,4] - Q(:,3) = [0,-1] - res = matmul(P, transpose(Q)) - print*, "" - do i = 1, 2 - print*, res(i, :) - end do - - - do j = 1, 2 - do i = 1, 2 - tmp(:) = P(i,:) * Q(j,:) - res(i,j) = stdlib_sum_kahan(tmp) - end do - end do - print*, "" - do i = 1, 2 - print*, res(i, :) - end do - end program example_sum \ No newline at end of file From b9b9832a50ef372b3eb0f834ac9dac7eec91326a Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Thu, 26 Feb 2026 21:09:03 +0530 Subject: [PATCH 18/29] change complex weights to real weights --- example/spatial/example_kabsch_umeyama.f90 | 1 + src/spatial/stdlib_spatial.fypp | 2 +- .../stdlib_spatial_kabsch_umeyama.fypp | 28 +++++++++++-------- 3 files changed, 19 insertions(+), 12 deletions(-) diff --git a/example/spatial/example_kabsch_umeyama.f90 b/example/spatial/example_kabsch_umeyama.f90 index 90c59b112..c0c3b71d4 100644 --- a/example/spatial/example_kabsch_umeyama.f90 +++ b/example/spatial/example_kabsch_umeyama.f90 @@ -24,6 +24,7 @@ program example_kabsch_umeyama print *, R(i,:) end do + print *, "" print *, "Recovered scale c:", c print *, "" diff --git a/src/spatial/stdlib_spatial.fypp b/src/spatial/stdlib_spatial.fypp index d12f6f383..1b56db2de 100644 --- a/src/spatial/stdlib_spatial.fypp +++ b/src/spatial/stdlib_spatial.fypp @@ -39,7 +39,7 @@ module stdlib_spatial !> Root-mean-square deviation real(${k}$), intent(out) :: rmsd !> Optional weights - ${t}$, intent(in), optional :: W(:) + real(${k}$), intent(in), optional :: W(:) !> Enable scaling logical, intent(in), optional :: scale end subroutine diff --git a/src/spatial/stdlib_spatial_kabsch_umeyama.fypp b/src/spatial/stdlib_spatial_kabsch_umeyama.fypp index 2c9f6dbba..ec4e3a7ae 100644 --- a/src/spatial/stdlib_spatial_kabsch_umeyama.fypp +++ b/src/spatial/stdlib_spatial_kabsch_umeyama.fypp @@ -22,7 +22,7 @@ contains !> Root-mean-square deviation real(${k}$), intent(out) :: rmsd !> Optional weights - ${t}$, intent(in), optional :: W(:) + real(${k}$), intent(in), optional :: W(:) !> Enable scaling logical, intent(in), optional :: scale @@ -59,12 +59,15 @@ contains allocate(c_P(d), source=zero_${s}$) allocate(c_Q(d), source=zero_${s}$) + allocate(tmp_N(N), source=zero_${s}$) ! Compute centroids of P and Q if(present(W)) then do i = 1, d - c_P(i) = stdlib_dot_product_kahan(w,P(i, :)) - c_Q(i) = stdlib_dot_product_kahan(w,Q(i, :)) + tmp_N(:) = W(:) * P(i,:) + c_P(i) = stdlib_sum_kahan(tmp_N) + tmp_N(:) = W(:) * Q(i,:) + c_Q(i) = stdlib_sum_kahan(tmp_N) end do else do i = 1, d @@ -77,7 +80,6 @@ contains ! Compute covariance matrix H = (P - c_P) * (Q - c_Q)^T and variance of P allocate(covariance(d,d), source=zero_${s}$) - allocate(tmp_N(N), source=zero_${s}$) allocate(tmp_d(d), source=zero_${s}$) variance_p = zero_${s}$ @@ -86,15 +88,16 @@ contains tmp_d = P(:, point) - c_P(:) tmp_N(point) = stdlib_dot_product_kahan(tmp_d, tmp_d) end do - variance_p = stdlib_dot_product_kahan(w, tmp_N) + tmp_N(:) = W(:) * tmp_N(:) + variance_p = stdlib_sum_kahan(tmp_N) do j = 1, d do i = 1, d #:if t.startswith('complex') - tmp_N(:) = (P(i,:) - c_P(i)) * conjg(Q(j,:) - c_Q(j)) - covariance(i,j) = stdlib_dot_product_kahan(w, tmp_N) + tmp_N(:) = W(:) * (P(i,:) - c_P(i)) * conjg(Q(j,:) - c_Q(j)) + covariance(i,j) = stdlib_sum_kahan(tmp_N) #:else - tmp_N(:) = (P(i,:) - c_P(i)) * (Q(j,:) - c_Q(j)) - covariance(i,j) = stdlib_dot_product_kahan(w, tmp_N) + tmp_N(:) = W(:) * (P(i,:) - c_P(i)) * (Q(j,:) - c_Q(j)) + covariance(i,j) = stdlib_sum_kahan(tmp_N) #:endif end do end do @@ -139,8 +142,11 @@ contains end do ! Scaling factor - c = variance_p / (sum(S(1:d))) - if (.not. scale_) c = one_${s}$ + if(scale_) then + c = variance_p / (sum(S(1:d))) + else + c = one_${s}$ + end if ! Translation vector t = c_P - c*R*c_Q do i = 1, d From 37b3cba8e5ad112f27c4123dc9eab5e6f1db9ff1 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Thu, 26 Feb 2026 21:22:11 +0530 Subject: [PATCH 19/29] minor changes --- src/spatial/stdlib_spatial_kabsch_umeyama.fypp | 13 ++++++------- test/spatial/CMakeLists.txt | 2 +- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/spatial/stdlib_spatial_kabsch_umeyama.fypp b/src/spatial/stdlib_spatial_kabsch_umeyama.fypp index ec4e3a7ae..3e081482b 100644 --- a/src/spatial/stdlib_spatial_kabsch_umeyama.fypp +++ b/src/spatial/stdlib_spatial_kabsch_umeyama.fypp @@ -28,8 +28,7 @@ contains ! Internal variables. integer(ilp) :: i, j, point, d, N - ${t}$, allocatable :: covariance(:,:), U(:,:), Vt(:,:), B(:,:), vec(:), tmp_N(:), tmp_d(:), c_P(:), c_Q(:) - ${t}$ :: vp, vq + ${t}$, allocatable :: covariance(:,:), U(:,:), Vt(:,:), vec(:), tmp_N(:), tmp_d(:), c_P(:), c_Q(:) real(${k}$) :: sum_w, variance_p real(${k}$), allocatable :: S(:) logical :: scale_ @@ -54,8 +53,8 @@ contains scale_ = .true. if(present(scale)) scale_ = scale - sum_w = one_${s}$ / N - if(present(W)) sum_w = one_${s}$ / stdlib_sum_kahan(W) + sum_w = one_${k}$ / N + if(present(W)) sum_w = one_${k}$ / stdlib_sum_kahan(W) allocate(c_P(d), source=zero_${s}$) allocate(c_Q(d), source=zero_${s}$) @@ -81,7 +80,7 @@ contains ! Compute covariance matrix H = (P - c_P) * (Q - c_Q)^T and variance of P allocate(covariance(d,d), source=zero_${s}$) allocate(tmp_d(d), source=zero_${s}$) - variance_p = zero_${s}$ + variance_p = zero_${k}$ if (present(W)) then do point = 1, N @@ -159,8 +158,8 @@ contains ! Compute RMSD allocate(vec(d), source=zero_${s}$) - rmsd = zero_${s}$ - rmsd_err = zero_${s}$ + rmsd = zero_${k}$ + rmsd_err = zero_${k}$ do point = 1, N ! Calculate the k^th difference vector by the formula vec_k = c*R*Q_k + t - P_k do i = 1, d diff --git a/test/spatial/CMakeLists.txt b/test/spatial/CMakeLists.txt index a64478b2d..6cd88dfd1 100644 --- a/test/spatial/CMakeLists.txt +++ b/test/spatial/CMakeLists.txt @@ -2,6 +2,6 @@ set( fppFiles "test_spatial_kabsch_umeyama.fypp" ) -fypp_f90pp("${fyppFlags}" "${fppFiles}" outFiles) +fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) ADDTESTPP(spatial_kabsch_umeyama) \ No newline at end of file From 723661b9de88cd7883c9040941bc5191ee646de0 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Thu, 26 Feb 2026 21:27:11 +0530 Subject: [PATCH 20/29] revert cmake change --- test/spatial/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/spatial/CMakeLists.txt b/test/spatial/CMakeLists.txt index 6cd88dfd1..a64478b2d 100644 --- a/test/spatial/CMakeLists.txt +++ b/test/spatial/CMakeLists.txt @@ -2,6 +2,6 @@ set( fppFiles "test_spatial_kabsch_umeyama.fypp" ) -fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) +fypp_f90pp("${fyppFlags}" "${fppFiles}" outFiles) ADDTESTPP(spatial_kabsch_umeyama) \ No newline at end of file From 3a7221fa5404b97317d797e5bee846d6b7aa8c7d Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Sat, 7 Mar 2026 18:31:29 +0530 Subject: [PATCH 21/29] add sum_w error check --- src/spatial/stdlib_spatial_kabsch_umeyama.fypp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/spatial/stdlib_spatial_kabsch_umeyama.fypp b/src/spatial/stdlib_spatial_kabsch_umeyama.fypp index 3e081482b..907732c69 100644 --- a/src/spatial/stdlib_spatial_kabsch_umeyama.fypp +++ b/src/spatial/stdlib_spatial_kabsch_umeyama.fypp @@ -55,6 +55,7 @@ contains sum_w = one_${k}$ / N if(present(W)) sum_w = one_${k}$ / stdlib_sum_kahan(W) + if(sum_w Date: Sun, 8 Mar 2026 16:03:06 +0530 Subject: [PATCH 22/29] reflection handling in real cases and modify tests --- src/spatial/stdlib_spatial.fypp | 4 +- .../stdlib_spatial_kabsch_umeyama.fypp | 26 +- test/spatial/test_spatial_kabsch_umeyama.fypp | 274 +++++++++++++----- 3 files changed, 225 insertions(+), 79 deletions(-) diff --git a/src/spatial/stdlib_spatial.fypp b/src/spatial/stdlib_spatial.fypp index 1b56db2de..4792d9a63 100644 --- a/src/spatial/stdlib_spatial.fypp +++ b/src/spatial/stdlib_spatial.fypp @@ -26,9 +26,9 @@ module stdlib_spatial !----------------------------------------------------------------------- #:for k, t, s in (KINDS_TYPES) module subroutine kabsch_umeyama_${s}$(P, Q, R, t, c, rmsd, W, scale) - !> Reference point set (d × N) - ${t}$, intent(in) :: P(:, :) !> Target point set (d × N) + ${t}$, intent(in) :: P(:, :) + !> Reference point set (d × N) ${t}$, intent(in) :: Q(:, :) !> Optimal rotation matrix (d × d) ${t}$, intent(out) :: R(:,:) diff --git a/src/spatial/stdlib_spatial_kabsch_umeyama.fypp b/src/spatial/stdlib_spatial_kabsch_umeyama.fypp index 907732c69..1ab1e0e58 100644 --- a/src/spatial/stdlib_spatial_kabsch_umeyama.fypp +++ b/src/spatial/stdlib_spatial_kabsch_umeyama.fypp @@ -3,15 +3,15 @@ #:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) #:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES submodule(stdlib_spatial) stdlib_spatial_kabsch_umeyama - use stdlib_linalg, only: svd + use stdlib_linalg, only: svd, det use stdlib_intrinsics, only: stdlib_sum_kahan, stdlib_dot_product_kahan, kahan_kernel contains #:for k, t, s in (KINDS_TYPES) module subroutine kabsch_umeyama_${s}$(P, Q, R, t, c, rmsd, W, scale) - !> Reference point set (d × N) - ${t}$, intent(in) :: P(:, :) !> Target point set (d × N) + ${t}$, intent(in) :: P(:, :) + !> Reference point set (d × N) ${t}$, intent(in) :: Q(:, :) !> Optimal rotation matrix (d × d) ${t}$, intent(out) :: R(:,:) @@ -32,6 +32,9 @@ contains real(${k}$) :: sum_w, variance_p real(${k}$), allocatable :: S(:) logical :: scale_ + #:if t.startswith('real') + logical :: reflect_ + #:endif real(${k}$) :: rmsd_err @@ -130,6 +133,15 @@ contains ! SVD of covariance matrix H -> H = U * S * Vt call svd(covariance, S, U, Vt) + ! Check for reflections in case of real entries. + #:if t.startswith('real') + reflect_ = det(matmul(U,Vt)) < zero_${s}$ + #:endif + + #:if t.startswith('real') + if(reflect_) Vt(d,:) = -Vt(d,:) + #:endif + ! Optimal rotation matrix. do i = 1,d do j = 1,d @@ -143,7 +155,15 @@ contains ! Scaling factor if(scale_) then + #:if t.startswith('real') + if(reflect_) then + c = variance_p / (sum(S(1:d-1)) - S(d)) + else + c = variance_p / (sum(S(1:d))) + end if + #:else c = variance_p / (sum(S(1:d))) + #:endif else c = one_${s}$ end if diff --git a/test/spatial/test_spatial_kabsch_umeyama.fypp b/test/spatial/test_spatial_kabsch_umeyama.fypp index 97b891a2b..2096ee3b0 100644 --- a/test/spatial/test_spatial_kabsch_umeyama.fypp +++ b/test/spatial/test_spatial_kabsch_umeyama.fypp @@ -5,7 +5,7 @@ module test_kabsch_umeyama use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test use stdlib_kinds use stdlib_math, only: all_close, is_close - use stdlib_linalg, only: svd + use stdlib_linalg, only: svd, det, eye use stdlib_spatial use stdlib_intrinsics, only: stdlib_dot_product_kahan use stdlib_constants @@ -25,58 +25,175 @@ contains ] end subroutine + #:for k, t, s in (R_KINDS_TYPES) + subroutine generate_random_set_real_${s}$(P, Q, d, N, scale) + integer, intent(in) :: d, N + ${t}$, intent(out) :: P(d, N), Q(d, N) + logical, intent(in) :: scale + + ! Internal variables. + ${t}$ :: R(d, d) + ${t}$ :: U(d,d), Vt(d,d) + real(${k}$) :: S(d) + ${t}$ :: t(d) + ${t}$ :: c + integer :: i, j + + c = one_${s}$ + ! Random reference points Q + call random_number(Q) + + ! Random proper rotation matrix R constructed via SVD: R = U * V^T + call random_number(R) + call svd(R, S, U, Vt) + do i = 1,d + do j = 1,d + R(i,j) = stdlib_dot_product_kahan(U(i,:), Vt(:, j)) + end do + end do + if(det(R) < zero_${s}$) R(:, d) = -R(:, d) + + ! Random scale and translation + if(scale) call random_number(c) + call random_number(t) + + ! Construct P = c*R*Q + t + do j = 1, N + do i = 1, d + P(i,j) = c * & + stdlib_dot_product_kahan(R(i,:), Q(:,j)) + & + t(i) + end do + end do + end subroutine + #:endfor + #:for k, t, s in (C_KINDS_TYPES) + subroutine generate_random_set_complex_${s}$(P, Q, d, N, scale) + integer, intent(in) :: d, N + ${t}$, intent(out) :: P(d, N), Q(d, N) + logical, intent(in) :: scale + + ! Internal variables. + ${t}$ :: R(d,d), t(d) + ${t}$ :: U(d,d), Vt(d,d) + real(${k}$) :: S(d) + ${t}$ :: c + real(${k}$) :: tempdn(d,N,2) + real(${k}$) :: tempdd(d,d,2) + real(${k}$) :: r1, r2 + integer :: i, j + + c = one_${s}$ + ! Random reference points Q + call random_number(tempdn) + Q = cmplx(tempdn(:,:,1), tempdn(:,:,2), kind=${k}$) + + ! Random proper rotation matrix R constructed via SVD: R = U * V^H + call random_number(tempdd) + R = cmplx(tempdd(:,:,1), tempdd(:,:,2), kind=${k}$) + call svd(R, S, U, Vt) + do i = 1,d + do j = 1,d + R(i,j) = stdlib_dot_product_kahan(U(i,:), Vt(:, j)) + end do + end do + + ! Random scale and translation + if(scale) then + call random_number(r1) + call random_number(r2) + c = cmplx(r1, r2, kind=${k}$) + end if + call random_number(tempdd) + t = cmplx(tempdd(:,1,1), tempdd(:,1,2), kind=${k}$) + + ! Construct P = c*R*Q + t + do j = 1, N + do i = 1, d + P(i,j) = c * & + stdlib_dot_product_kahan(R(i,:), Q(:,j)) + & + t(i) + end do + end do + end subroutine + #:endfor + subroutine test_kabsch_umeyama_real(error) type(error_type), allocatable, intent(out) :: error #:for k, t, s in (R_KINDS_TYPES) block integer, parameter :: d = 3, N = 8 ${t}$ :: P_original(d, N), P_recovered(d, N), Q_original(d, N) - ${t}$ :: R_recovered(d, d), R_original(d, d) - ${t}$ :: t_recovered(d), t_original(d) - ${t}$ :: c_recovered, c_original - ${t}$ :: U(d,d), Vt(d,d), S(d) + ${t}$ :: R_recovered(d, d) + ${t}$ :: t_recovered(d) + ${t}$ :: c_recovered ${t}$ :: rmsd + ${t}$ :: Identity(d,d) + ${t}$ :: zero_t(d) + ${t}$ :: W(N) integer :: i, j - ${t}$ :: r1 real(${k}$), parameter :: ${k}$tol = 100 * epsilon(one_${k}$) R_recovered = zero_${s}$ t_recovered = zero_${s}$ c_recovered = zero_${s}$ - U = zero_${s}$ - S = zero_${k}$ - Vt = zero_${s}$ rmsd = zero_${k}$ + Identity = eye(d,d) + zero_t = zero_${s}$ - ! Random reference points Q - call random_number(Q_original) - - ! Random proper rotation matrix R_original constructed via SVD: R = U * V^T - call random_number(R_original) - call svd(R_original, S, U, Vt) - do i = 1,d - do j = 1,d - R_original(i,j) = stdlib_dot_product_kahan(U(i,:), Vt(:, j)) + ! Random set of transformation + call generate_random_set_real_${s}$(P_original, Q_original, d, N, .true.) + ! Call Kabsch–Umeyama + call kabsch_umeyama(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) + ! Apply recovered transform + do j = 1, N + do i = 1, d + P_recovered(i,j) = c_recovered * & + stdlib_dot_product_kahan(R_recovered(i,:), Q_original(:,j)) + & + t_recovered(i) end do end do + call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) + if (allocated(error)) return + call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) + if(allocated(error)) return - ! Random scale and translation - call random_number(c_original) - call random_number(t_original) + ! Identity transformation (Q -> Q) + call random_number(Q_original) + ! Call Kabsch–Umeyama + call kabsch_umeyama(Q_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) + call check(error, all_close(R_recovered, Identity, abs_tol=${k}$tol), .true.) + if(allocated(error)) return + call check(error, all_close(t_recovered, zero_t, abs_tol = ${k}$tol), .true.) + if (allocated(error)) return + call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) + if(allocated(error)) return - ! Construct P = c*R*Q + t + ! Random set of transformation but with scale disabled. + call generate_random_set_real_${s}$(P_original, Q_original, d, N, .false.) + ! Call Kabsch–Umeyama + call kabsch_umeyama(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd, scale=.false.) + ! Apply recovered transform do j = 1, N do i = 1, d - P_original(i,j) = c_original * & - stdlib_dot_product_kahan(R_original(i,:), Q_original(:,j)) + & - t_original(i) + P_recovered(i,j) = c_recovered * & + stdlib_dot_product_kahan(R_recovered(i,:), Q_original(:,j)) + & + t_recovered(i) end do end do + call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) + if (allocated(error)) return + call check(error, is_close(c_recovered, one_${k}$, abs_tol = ${k}$tol), .true.) + if(allocated(error)) return + call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) + if(allocated(error)) return + ! Random set of transformation but with weights. Since P = c*R*Q+t, weights should not matter + call random_number(W) + call generate_random_set_real_${s}$(P_original, Q_original, d, N, .true.) ! Call Kabsch–Umeyama - call kabsch_umeyama(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) - + call kabsch_umeyama(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd, W=W) ! Apply recovered transform do j = 1, N do i = 1, d @@ -85,7 +202,8 @@ contains t_recovered(i) end do end do - + call check(error, is_close(det(R_recovered), one_${s}$, abs_tol=${k}$tol), .true.) + if(allocated(error)) return call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) if (allocated(error)) return call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) @@ -101,73 +219,80 @@ contains block integer, parameter :: d = 3, N = 8 ${t}$ :: P_original(d, N), Q_original(d, N), P_recovered(d, N) - ${t}$ :: R_recovered(d, d), R_original(d, d) - ${t}$ :: t_recovered(d), t_original(d) - ${t}$ :: c_recovered, c_original - ${t}$ :: U(d,d), Vt(d,d) - real(${k}$) :: S(d) + ${t}$ :: R_recovered(d, d) + ${t}$ :: t_recovered(d) + ${t}$ :: c_recovered real(${k}$) :: rmsd + ${t}$ :: Identity(d,d) + ${t}$ :: zero_t(d) + real(${k}$) :: W(N) integer :: i, j - real(${k}$) :: r1, r2 real(${k}$), parameter :: ${k}$tol = 100 * epsilon(one_${k}$) R_recovered = zero_${s}$ t_recovered = zero_${s}$ c_recovered = zero_${s}$ - U = zero_${s}$ - S = zero_${k}$ - Vt = zero_${s}$ rmsd = zero_${k}$ + Identity = eye(d,d) + zero_t = zero_${s}$ - ! Random complex reference points Q + ! Random set of transformation + call generate_random_set_complex_${s}$(P_original, Q_original, d, N, .true.) + ! Call Kabsch–Umeyama + call kabsch_umeyama(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) + ! Apply recovered transform do j = 1, N do i = 1, d - call random_number(r1) - call random_number(r2) - Q_original(i,j) = cmplx(r1, r2, kind=${k}$) - end do - end do - - ! Random complex unitary matrix R_original - ! Constructed via SVD: R = U * V^H - do i = 1, d - do j = 1, d - call random_number(r1) - call random_number(r2) - R_original(i,j) = cmplx(r1, r2, kind=${k}$) - end do - end do - call svd(R_original, S, U, Vt) - do i = 1,d - do j = 1,d - R_original(i,j) = stdlib_dot_product_kahan(conjg(U(i,:)), Vt(:, j)) + P_recovered(i,j) = c_recovered * & + stdlib_dot_product_kahan(conjg(R_recovered(i,:)), Q_original(:,j)) + & + t_recovered(i) end do end do + call check(error, is_close(abs(det(R_recovered)), one_${k}$, abs_tol=${k}$tol), .true.) + if(allocated(error)) return + call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) + if (allocated(error)) return + call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) + if(allocated(error)) return - ! Random scale and translation - call random_number(r1) - call random_number(r2) - c_original = cmplx(r1,r2, kind=${k}$) - - do i = 1, d - call random_number(r1) - call random_number(r2) - t_original(i) = cmplx(r1, r2, kind=${k}$) - end do + ! Identity transformation (Q -> Q) + call generate_random_set_complex_${s}$(P_original, Q_original, d, N, .true.) + ! Call Kabsch–Umeyama + call kabsch_umeyama(Q_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) + call check(error, all_close(R_recovered, Identity, abs_tol=${k}$tol), .true.) + if(allocated(error)) return + call check(error, all_close(t_recovered, zero_t, abs_tol = ${k}$tol), .true.) + if (allocated(error)) return + call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) + if(allocated(error)) return - ! Construct P = c*R*Q + t + ! Random set of transformation but with scale disabled. + call generate_random_set_complex_${s}$(P_original, Q_original, d, N, .false.) + ! Call Kabsch–Umeyama + call kabsch_umeyama(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd,scale=.false.) + ! Apply recovered transform do j = 1, N do i = 1, d - P_original(i,j) = c_original * & - stdlib_dot_product_kahan(conjg(R_original(i,:)), Q_original(:,j)) + & - t_original(i) + P_recovered(i,j) = c_recovered * & + stdlib_dot_product_kahan(conjg(R_recovered(i,:)), Q_original(:,j)) + & + t_recovered(i) end do end do + call check(error, is_close(abs(det(R_recovered)), one_${k}$, abs_tol=${k}$tol), .true.) + if(allocated(error)) return + call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) + if (allocated(error)) return + call check(error, is_close(c_recovered, one_${s}$, abs_tol = ${k}$tol), .true.) + if(allocated(error)) return + call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) + if(allocated(error)) return + ! Random set of transformation but with weights. Since P = c*R*Q+t, weights should not matter + call random_number(W) + call generate_random_set_complex_${s}$(P_original, Q_original, d, N, .true.) ! Call Kabsch–Umeyama - call kabsch_umeyama(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) - + call kabsch_umeyama(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd, W=W) ! Apply recovered transform do j = 1, N do i = 1, d @@ -176,7 +301,8 @@ contains t_recovered(i) end do end do - + call check(error, is_close(abs(det(R_recovered)), one_${k}$, abs_tol=${k}$tol), .true.) + if(allocated(error)) return call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) if (allocated(error)) return call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) From 7554914e2e74c57447b161d221725f32df42d042 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Sun, 8 Mar 2026 16:11:05 +0530 Subject: [PATCH 23/29] add target stdlib_core to spatial/CMakeLists.txt --- src/spatial/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spatial/CMakeLists.txt b/src/spatial/CMakeLists.txt index 58c671e7f..cebe6e111 100644 --- a/src/spatial/CMakeLists.txt +++ b/src/spatial/CMakeLists.txt @@ -11,4 +11,4 @@ set(spatial_f90Files configure_stdlib_target(${PROJECT_NAME}_spatial spatial_f90Files spatial_fppFiles spatial_cppFiles) -target_link_libraries(${PROJECT_NAME}_spatial PUBLIC ${PROJECT_NAME}_constants ${PROJECT_NAME}_linalg_core ${PROJECT_NAME}_linalg ${PROJECT_NAME}_intrinsics) \ No newline at end of file +target_link_libraries(${PROJECT_NAME}_spatial PUBLIC ${PROJECT_NAME}_core ${PROJECT_NAME}_constants ${PROJECT_NAME}_linalg_core ${PROJECT_NAME}_linalg ${PROJECT_NAME}_intrinsics) \ No newline at end of file From 12f1513fe5f16fa4bad72252612cd37d6394f859 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Sun, 8 Mar 2026 16:21:43 +0530 Subject: [PATCH 24/29] modify test/spatial/CMakeLists.txt --- test/spatial/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/spatial/CMakeLists.txt b/test/spatial/CMakeLists.txt index a64478b2d..081d83afb 100644 --- a/test/spatial/CMakeLists.txt +++ b/test/spatial/CMakeLists.txt @@ -2,6 +2,6 @@ set( fppFiles "test_spatial_kabsch_umeyama.fypp" ) -fypp_f90pp("${fyppFlags}" "${fppFiles}" outFiles) +fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) -ADDTESTPP(spatial_kabsch_umeyama) \ No newline at end of file +ADDTEST(spatial_kabsch_umeyama) \ No newline at end of file From fe7d96e0f4fe43dada2c4c896a8db7ca0fcde833 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Wed, 1 Apr 2026 21:16:31 +0530 Subject: [PATCH 25/29] modify: example, add: comment in test, minor changes --- example/spatial/example_kabsch_umeyama.f90 | 29 +++++++++++-------- src/spatial/stdlib_spatial.fypp | 2 +- .../stdlib_spatial_kabsch_umeyama.fypp | 3 +- test/spatial/test_spatial_kabsch_umeyama.fypp | 2 ++ 4 files changed, 22 insertions(+), 14 deletions(-) diff --git a/example/spatial/example_kabsch_umeyama.f90 b/example/spatial/example_kabsch_umeyama.f90 index c0c3b71d4..8986b9a01 100644 --- a/example/spatial/example_kabsch_umeyama.f90 +++ b/example/spatial/example_kabsch_umeyama.f90 @@ -3,35 +3,40 @@ program example_kabsch_umeyama use stdlib_spatial, only: kabsch_umeyama implicit none - integer, parameter :: d = 2, N = 3 - real(dp) :: P(d, N), Q(d, N), R(d, d), t(d), c, rmsd + integer, parameter :: d = 2, N = 7 + real(dp) :: P(d,N), Q(d,N), R(d, d), t(d), c, rmsd integer :: i - P(:,1) = [3.0_dp, -2.0_dp] - P(:,2) = [7.0_dp, 4.0_dp] - P(:,3) = [5.0_dp, 0.0_dp] + ! 2x7 matrices. + P(1,:) = [23.0_dp, 66.0_dp, 88.0_dp, 119.0_dp, 122.0_dp, 170.0_dp, 179.0_dp] + P(2,:) = [178.0_dp, 173.0_dp, 187.0_dp, 202.0_dp, 229.0_dp, 232.0_dp, 199.0_dp] - Q(:,1) = [2.0_dp, 3.0_dp] - Q(:,2) = [-1.0_dp, 5.0_dp] - Q(:,3) = [1.0_dp, 4.0_dp] + Q(1,:) = [232.0_dp, 208.0_dp, 181.0_dp, 155.0_dp, 142.0_dp, 121.0_dp, 139.0_dp] + Q(2,:) = [ 38.0_dp, 32.0_dp, 31.0_dp, 45.0_dp, 33.0_dp, 59.0_dp, 69.0_dp] call kabsch_umeyama(P, Q, R, t, c, rmsd) - print *, "" + print * print *, "Recovered rotation R:" do i = 1, d print *, R(i,:) end do - print *, "" + print * print *, "Recovered scale c:", c - print *, "" + print * print *, "Recovered translation t:" print *, t - print *, "" + print * print *, "RMSD:", rmsd + print * + print *, "Recovered P:" + do i = 1, N + print*, c*matmul(R, Q(:,i)) + t + end do + end program example_kabsch_umeyama \ No newline at end of file diff --git a/src/spatial/stdlib_spatial.fypp b/src/spatial/stdlib_spatial.fypp index 4792d9a63..b8638958f 100644 --- a/src/spatial/stdlib_spatial.fypp +++ b/src/spatial/stdlib_spatial.fypp @@ -31,7 +31,7 @@ module stdlib_spatial !> Reference point set (d × N) ${t}$, intent(in) :: Q(:, :) !> Optimal rotation matrix (d × d) - ${t}$, intent(out) :: R(:,:) + ${t}$, intent(out) :: R(:, :) !> Translation vector (d) ${t}$, intent(out) :: t(:) !> Scale factor diff --git a/src/spatial/stdlib_spatial_kabsch_umeyama.fypp b/src/spatial/stdlib_spatial_kabsch_umeyama.fypp index 1ab1e0e58..851d9288f 100644 --- a/src/spatial/stdlib_spatial_kabsch_umeyama.fypp +++ b/src/spatial/stdlib_spatial_kabsch_umeyama.fypp @@ -14,7 +14,7 @@ contains !> Reference point set (d × N) ${t}$, intent(in) :: Q(:, :) !> Optimal rotation matrix (d × d) - ${t}$, intent(out) :: R(:,:) + ${t}$, intent(out) :: R(:, :) !> Translation vector (d) ${t}$, intent(out) :: t(:) !> Scale factor @@ -105,6 +105,7 @@ contains end do end do else + ! Calculate variance by the formula (1/n)*sigma(P - c_P)^2 do point = 1, N tmp_d = P(:, point) - c_P(:) tmp_N(point) = stdlib_dot_product_kahan(tmp_d, tmp_d) diff --git a/test/spatial/test_spatial_kabsch_umeyama.fypp b/test/spatial/test_spatial_kabsch_umeyama.fypp index 2096ee3b0..c2b83ad2b 100644 --- a/test/spatial/test_spatial_kabsch_umeyama.fypp +++ b/test/spatial/test_spatial_kabsch_umeyama.fypp @@ -51,6 +51,8 @@ contains R(i,j) = stdlib_dot_product_kahan(U(i,:), Vt(:, j)) end do end do + ! For real cases, only det(R) = +1 can be checked. Because if det(R) = -1, then algorithm will + ! find a rotation matrix without reflection in which rmsd is minimized. if(det(R) < zero_${s}$) R(:, d) = -R(:, d) ! Random scale and translation From 81a228c7478cad1613e95a358b47061a6b187c34 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Wed, 1 Apr 2026 22:33:55 +0530 Subject: [PATCH 26/29] modify: dependencies --- src/spatial/CMakeLists.txt | 2 +- src/spatial/stdlib_spatial.fypp | 3 +-- src/spatial/stdlib_spatial_kabsch_umeyama.fypp | 3 ++- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/spatial/CMakeLists.txt b/src/spatial/CMakeLists.txt index cebe6e111..1cfd0e111 100644 --- a/src/spatial/CMakeLists.txt +++ b/src/spatial/CMakeLists.txt @@ -11,4 +11,4 @@ set(spatial_f90Files configure_stdlib_target(${PROJECT_NAME}_spatial spatial_f90Files spatial_fppFiles spatial_cppFiles) -target_link_libraries(${PROJECT_NAME}_spatial PUBLIC ${PROJECT_NAME}_core ${PROJECT_NAME}_constants ${PROJECT_NAME}_linalg_core ${PROJECT_NAME}_linalg ${PROJECT_NAME}_intrinsics) \ No newline at end of file +target_link_libraries(${PROJECT_NAME}_spatial PUBLIC ${PROJECT_NAME}_core ${PROJECT_NAME}_constants ${PROJECT_NAME}_linalg ${PROJECT_NAME}_intrinsics) \ No newline at end of file diff --git a/src/spatial/stdlib_spatial.fypp b/src/spatial/stdlib_spatial.fypp index b8638958f..bdb0014d0 100644 --- a/src/spatial/stdlib_spatial.fypp +++ b/src/spatial/stdlib_spatial.fypp @@ -3,9 +3,8 @@ #:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) #:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES module stdlib_spatial - use stdlib_linalg_constants + use stdlib_kinds, only: sp, dp, xdp, qp use stdlib_constants - use stdlib_error, only: error_stop implicit none private public :: kabsch_umeyama diff --git a/src/spatial/stdlib_spatial_kabsch_umeyama.fypp b/src/spatial/stdlib_spatial_kabsch_umeyama.fypp index 851d9288f..59fa659a1 100644 --- a/src/spatial/stdlib_spatial_kabsch_umeyama.fypp +++ b/src/spatial/stdlib_spatial_kabsch_umeyama.fypp @@ -5,6 +5,7 @@ submodule(stdlib_spatial) stdlib_spatial_kabsch_umeyama use stdlib_linalg, only: svd, det use stdlib_intrinsics, only: stdlib_sum_kahan, stdlib_dot_product_kahan, kahan_kernel + use stdlib_error, only: error_stop contains #:for k, t, s in (KINDS_TYPES) @@ -27,7 +28,7 @@ contains logical, intent(in), optional :: scale ! Internal variables. - integer(ilp) :: i, j, point, d, N + integer :: i, j, point, d, N ${t}$, allocatable :: covariance(:,:), U(:,:), Vt(:,:), vec(:), tmp_N(:), tmp_d(:), c_P(:), c_Q(:) real(${k}$) :: sum_w, variance_p real(${k}$), allocatable :: S(:) From a9f0fc19d759e54d7fa2c6424a9ed6a5491129d3 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Thu, 2 Apr 2026 22:50:03 +0530 Subject: [PATCH 27/29] modify: docs, ford comments and example --- doc/specs/stdlib_spatial.md | 59 +++++++++++++++++++ example/spatial/example_kabsch_umeyama.f90 | 10 +--- src/spatial/stdlib_spatial.fypp | 36 +++++------ .../stdlib_spatial_kabsch_umeyama.fypp | 16 ++--- 4 files changed, 84 insertions(+), 37 deletions(-) create mode 100644 doc/specs/stdlib_spatial.md diff --git a/doc/specs/stdlib_spatial.md b/doc/specs/stdlib_spatial.md new file mode 100644 index 000000000..61fe8c219 --- /dev/null +++ b/doc/specs/stdlib_spatial.md @@ -0,0 +1,59 @@ +--- +title: spatial +--- + +# The `stdlib_spatial` module + +This module provides implementations of algorithms for spatial data processing. + +[TOC] + +## `kabsch_umeyama` - Finding optimal rotation matrix + +### Status + +Experimental + +### Description + +Compute the optimal similarity transform (Kabsch–Umeyama): +\[ + P \approx c \, R \, Q + t +\] +where: + +- R is an orthogonal rotation matrix, +- c is an optional scaling factor, +- t is a translation vector. + +The transformation minimizes the root-mean-square deviation (RMSD) between corresponding columns +of P and Q, optionally using weights and with optional scaling. +The implementation is based on the algorithm described here: [Aligning point patterns with Kabsch–Umeyama algorithm](https://zpl.fi/aligning-point-patterns-with-kabsch-umeyama-algorithm) + +### Syntax + +`call ` [[stdlib_spatial(module):kabsch_umeyama(interface)]] `(P, Q, R, t, c, rmsd [, W, scale])` + +### Arguments + +`P`: Shall be a `real` or `complex` rank-2 array. It is an `intent(in)` argument. + +`Q`: Shall be a rank-2 array with same kind as `P`. It is an `intent(in)` argument. + +`R`: Shall be a rank-2 array with same kind as `P`. For `real` kinds, the algorithm returns a proper rotation matrix, meaning `det(R) = 1`. It is an `intent(out)` argument. + +`t`: Shall be a rank-1 array with same kind as `P`. It is an `intent(out)` argument. + +`c`: Scalar value of the same type as `P`. It is an `intent(out)` argument. If `scale` is disabled `c` will be returned with a value of `1`. + +`rmsd`: Scalar value of `real` kind. It is an `intent(out)` argument. + +`W` (optional): Shall be a rank-1 array of `real` kind. It is an `intent(in)` argument. By default, `W` is an array of `1`s. + +`scale` (optional): Shall be a logical type. It is an `intent(in)` argument. By default, `scale = .true.`. + +### Example + +```fortran +{!example/spatial/example_kabsch_umeyama.f90!} +``` \ No newline at end of file diff --git a/example/spatial/example_kabsch_umeyama.f90 b/example/spatial/example_kabsch_umeyama.f90 index 8986b9a01..daa9f5d1f 100644 --- a/example/spatial/example_kabsch_umeyama.f90 +++ b/example/spatial/example_kabsch_umeyama.f90 @@ -9,11 +9,11 @@ program example_kabsch_umeyama integer :: i ! 2x7 matrices. - P(1,:) = [23.0_dp, 66.0_dp, 88.0_dp, 119.0_dp, 122.0_dp, 170.0_dp, 179.0_dp] + P(1,:) = [23.0_dp, 66.0_dp, 88.0_dp, 119.0_dp, 122.0_dp, 170.0_dp, 179.0_dp] P(2,:) = [178.0_dp, 173.0_dp, 187.0_dp, 202.0_dp, 229.0_dp, 232.0_dp, 199.0_dp] Q(1,:) = [232.0_dp, 208.0_dp, 181.0_dp, 155.0_dp, 142.0_dp, 121.0_dp, 139.0_dp] - Q(2,:) = [ 38.0_dp, 32.0_dp, 31.0_dp, 45.0_dp, 33.0_dp, 59.0_dp, 69.0_dp] + Q(2,:) = [38.0_dp, 32.0_dp, 31.0_dp, 45.0_dp, 33.0_dp, 59.0_dp, 69.0_dp] call kabsch_umeyama(P, Q, R, t, c, rmsd) @@ -33,10 +33,4 @@ program example_kabsch_umeyama print * print *, "RMSD:", rmsd - print * - print *, "Recovered P:" - do i = 1, N - print*, c*matmul(R, Q(:,i)) + t - end do - end program example_kabsch_umeyama \ No newline at end of file diff --git a/src/spatial/stdlib_spatial.fypp b/src/spatial/stdlib_spatial.fypp index bdb0014d0..33a15192a 100644 --- a/src/spatial/stdlib_spatial.fypp +++ b/src/spatial/stdlib_spatial.fypp @@ -10,37 +10,31 @@ module stdlib_spatial public :: kabsch_umeyama interface kabsch_umeyama - !----------------------------------------------------------------------- - !> Compute the optimal similarity transform (Kabsch–Umeyama): - !> - !> P ≈ c * R * Q + t - !> - !> where: - !> - R is an orthogonal rotation matrix - !> - c is an optional scale factor - !> - t is a translation vector - !> - !> The transformation minimizes the RMSD between corresponding columns - !> of P and Q, optionally using weights. - !----------------------------------------------------------------------- + !! ([Specifications](../page/specs/stdlib_spatial.html#kabsch_umeyama)) + !! This interface computes the optimal similarity transform (Kabsch–Umeyama): + !! \[ + !! P \approx c \, R \, Q + t + !! \] + !! The transformation minimizes the RMSD between corresponding columns + !! of P and Q, optionally using weights and with optional scaling. #:for k, t, s in (KINDS_TYPES) module subroutine kabsch_umeyama_${s}$(P, Q, R, t, c, rmsd, W, scale) - !> Target point set (d × N) ${t}$, intent(in) :: P(:, :) - !> Reference point set (d × N) + !! Target point set (d × N) ${t}$, intent(in) :: Q(:, :) - !> Optimal rotation matrix (d × d) + !! Reference point set (d × N) ${t}$, intent(out) :: R(:, :) - !> Translation vector (d) + !! Optimal rotation matrix (d × d) ${t}$, intent(out) :: t(:) - !> Scale factor + !! Translation vector (d) ${t}$, intent(out) :: c - !> Root-mean-square deviation + !! Scale factor real(${k}$), intent(out) :: rmsd - !> Optional weights + !! Root-mean-square deviation real(${k}$), intent(in), optional :: W(:) - !> Enable scaling + !! Optional weights logical, intent(in), optional :: scale + !! Enable scaling end subroutine #:endfor end interface diff --git a/src/spatial/stdlib_spatial_kabsch_umeyama.fypp b/src/spatial/stdlib_spatial_kabsch_umeyama.fypp index 59fa659a1..5de05309a 100644 --- a/src/spatial/stdlib_spatial_kabsch_umeyama.fypp +++ b/src/spatial/stdlib_spatial_kabsch_umeyama.fypp @@ -10,22 +10,22 @@ submodule(stdlib_spatial) stdlib_spatial_kabsch_umeyama contains #:for k, t, s in (KINDS_TYPES) module subroutine kabsch_umeyama_${s}$(P, Q, R, t, c, rmsd, W, scale) - !> Target point set (d × N) ${t}$, intent(in) :: P(:, :) - !> Reference point set (d × N) + !! Target point set (d × N) ${t}$, intent(in) :: Q(:, :) - !> Optimal rotation matrix (d × d) + !! Reference point set (d × N) ${t}$, intent(out) :: R(:, :) - !> Translation vector (d) + !! Optimal rotation matrix (d × d) ${t}$, intent(out) :: t(:) - !> Scale factor + !! Translation vector (d) ${t}$, intent(out) :: c - !> Root-mean-square deviation + !! Scale factor real(${k}$), intent(out) :: rmsd - !> Optional weights + !! Root-mean-square deviation real(${k}$), intent(in), optional :: W(:) - !> Enable scaling + !! Optional weights logical, intent(in), optional :: scale + !! Enable scaling ! Internal variables. integer :: i, j, point, d, N From 42db95b84c407008247e1683628c1c31f0537ed1 Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Thu, 2 Apr 2026 22:57:30 +0530 Subject: [PATCH 28/29] minor changes --- test/spatial/test_spatial_kabsch_umeyama.fypp | 1 - 1 file changed, 1 deletion(-) diff --git a/test/spatial/test_spatial_kabsch_umeyama.fypp b/test/spatial/test_spatial_kabsch_umeyama.fypp index c2b83ad2b..c8812287b 100644 --- a/test/spatial/test_spatial_kabsch_umeyama.fypp +++ b/test/spatial/test_spatial_kabsch_umeyama.fypp @@ -13,7 +13,6 @@ module test_kabsch_umeyama contains - !> Collect all exported unit tests subroutine collect_suite(testsuite) !> Collection of tests From fa159bdddde992ab20eb1e2a34c31b1323606f9f Mon Sep 17 00:00:00 2001 From: Mahmood-Sinan Date: Fri, 3 Apr 2026 17:33:26 +0530 Subject: [PATCH 29/29] add: debug messages for test --- test/spatial/test_spatial_kabsch_umeyama.fypp | 94 ++++++++++++------- 1 file changed, 59 insertions(+), 35 deletions(-) diff --git a/test/spatial/test_spatial_kabsch_umeyama.fypp b/test/spatial/test_spatial_kabsch_umeyama.fypp index c8812287b..7c2e3ec00 100644 --- a/test/spatial/test_spatial_kabsch_umeyama.fypp +++ b/test/spatial/test_spatial_kabsch_umeyama.fypp @@ -155,20 +155,25 @@ contains t_recovered(i) end do end do - call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) - if (allocated(error)) return - call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) + call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation, P_recovered not equal to P_original") + if(allocated(error)) return + call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation, RMSD not equal to zero") if(allocated(error)) return ! Identity transformation (Q -> Q) call random_number(Q_original) ! Call Kabsch–Umeyama call kabsch_umeyama(Q_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) - call check(error, all_close(R_recovered, Identity, abs_tol=${k}$tol), .true.) + call check(error, all_close(R_recovered, Identity, abs_tol=${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Identity transformation, R_recovered not equal to Identity") + if(allocated(error)) return + call check(error, all_close(t_recovered, zero_t, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Identity transformation, t_recovered not equal to zero") if(allocated(error)) return - call check(error, all_close(t_recovered, zero_t, abs_tol = ${k}$tol), .true.) - if (allocated(error)) return - call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) + call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Identity transformation, RMSD not equal to zero") if(allocated(error)) return ! Random set of transformation but with scale disabled. @@ -183,14 +188,17 @@ contains t_recovered(i) end do end do - call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) - if (allocated(error)) return - call check(error, is_close(c_recovered, one_${k}$, abs_tol = ${k}$tol), .true.) + call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation (scale disabled), P_recovered not equal to P_original") + if(allocated(error)) return + call check(error, is_close(c_recovered, one_${k}$, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation (scale disabled), c_recovered not equal to one") if(allocated(error)) return - call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) + call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation (scale disabled), RMSD not equal to zero") if(allocated(error)) return - ! Random set of transformation but with weights. Since P = c*R*Q+t, weights should not matter + ! Random set of transformation but with weights. Since P = c*R*Q + t, weights should not matter call random_number(W) call generate_random_set_real_${s}$(P_original, Q_original, d, N, .true.) ! Call Kabsch–Umeyama @@ -203,11 +211,14 @@ contains t_recovered(i) end do end do - call check(error, is_close(det(R_recovered), one_${s}$, abs_tol=${k}$tol), .true.) + call check(error, is_close(det(R_recovered), one_${s}$, abs_tol=${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation (weighted), det(R_recovered) not equal to one") if(allocated(error)) return - call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) - if (allocated(error)) return - call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) + call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation (weighted), P_recovered not equal to P_original") + if(allocated(error)) return + call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation (weighted), RMSD not equal to zero") if(allocated(error)) return end block #:endfor @@ -250,28 +261,34 @@ contains t_recovered(i) end do end do - call check(error, is_close(abs(det(R_recovered)), one_${k}$, abs_tol=${k}$tol), .true.) + call check(error, is_close(abs(det(R_recovered)), one_${k}$, abs_tol=${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation, abs(det(R_recovered)) not equal to one") + if(allocated(error)) return + call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation, P_recovered not equal to P_original") if(allocated(error)) return - call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) - if (allocated(error)) return - call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) + call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation, RMSD not equal to zero") if(allocated(error)) return ! Identity transformation (Q -> Q) call generate_random_set_complex_${s}$(P_original, Q_original, d, N, .true.) ! Call Kabsch–Umeyama call kabsch_umeyama(Q_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd) - call check(error, all_close(R_recovered, Identity, abs_tol=${k}$tol), .true.) + call check(error, all_close(R_recovered, Identity, abs_tol=${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Identity transformation, R_recovered not equal to Identity") if(allocated(error)) return - call check(error, all_close(t_recovered, zero_t, abs_tol = ${k}$tol), .true.) - if (allocated(error)) return - call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) + call check(error, all_close(t_recovered, zero_t, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Identity transformation, t_recovered not equal to zero") + if(allocated(error)) return + call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Identity transformation, RMSD not equal to zero") if(allocated(error)) return ! Random set of transformation but with scale disabled. call generate_random_set_complex_${s}$(P_original, Q_original, d, N, .false.) ! Call Kabsch–Umeyama - call kabsch_umeyama(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd,scale=.false.) + call kabsch_umeyama(P_original, Q_original, R_recovered, t_recovered, c_recovered, rmsd, scale=.false.) ! Apply recovered transform do j = 1, N do i = 1, d @@ -280,16 +297,20 @@ contains t_recovered(i) end do end do - call check(error, is_close(abs(det(R_recovered)), one_${k}$, abs_tol=${k}$tol), .true.) + call check(error, is_close(abs(det(R_recovered)), one_${k}$, abs_tol=${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation (scale disabled), abs(det(R_recovered)) not equal to one") if(allocated(error)) return - call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) - if (allocated(error)) return - call check(error, is_close(c_recovered, one_${s}$, abs_tol = ${k}$tol), .true.) + call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation (scale disabled), P_recovered not equal to P_original") if(allocated(error)) return - call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) + call check(error, is_close(c_recovered, one_${s}$, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation (scale disabled), c_recovered not equal to one") + if(allocated(error)) return + call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation (scale disabled), RMSD not equal to zero") if(allocated(error)) return - ! Random set of transformation but with weights. Since P = c*R*Q+t, weights should not matter + ! Random set of transformation but with weights. Since P = c*R*Q + t, weights should not matter call random_number(W) call generate_random_set_complex_${s}$(P_original, Q_original, d, N, .true.) ! Call Kabsch–Umeyama @@ -302,11 +323,14 @@ contains t_recovered(i) end do end do - call check(error, is_close(abs(det(R_recovered)), one_${k}$, abs_tol=${k}$tol), .true.) + call check(error, is_close(abs(det(R_recovered)), one_${k}$, abs_tol=${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation (weighted), abs(det(R_recovered)) not equal to one") + if(allocated(error)) return + call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation (weighted), P_recovered not equal to P_original") if(allocated(error)) return - call check(error, all_close(P_recovered, P_original, abs_tol = ${k}$tol), .true.) - if (allocated(error)) return - call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true.) + call check(error, is_close(rmsd, zero_${k}$, abs_tol = ${k}$tol), .true., & + "Kabsch_Umeyama ${t}$ failed: Random transformation (weighted), RMSD not equal to zero") if(allocated(error)) return end block #:endfor