diff --git a/autd3-gain-holo/src/combinatorial/greedy.rs b/autd3-gain-holo/src/combinatorial/greedy.rs index f5797427..434aecf6 100644 --- a/autd3-gain-holo/src/combinatorial/greedy.rs +++ b/autd3-gain-holo/src/combinatorial/greedy.rs @@ -4,7 +4,7 @@ * Created Date: 03/06/2021 * Author: Shun Suzuki * ----- - * Last Modified: 02/12/2023 + * Last Modified: 09/01/2024 * Modified By: Shun Suzuki (suzuki@hapis.k.u-tokyo.ac.jp) * ----- * Copyright (c) 2021 Shun Suzuki. All rights reserved. @@ -19,7 +19,7 @@ use autd3_derive::Gain; use autd3_driver::{ acoustics::{directivity::Sphere, propagate}, common::{EmitIntensity, Phase}, - defined::{PI, T4010A1_AMPLITUDE}, + defined::PI, derive::prelude::*, geometry::{Geometry, Vector3}, }; @@ -63,7 +63,7 @@ impl Greedy { res: &mut [Complex], ) { res.iter_mut().zip(foci.iter()).for_each(|(r, f)| { - *r = propagate::(trans, attenuation, sound_speed, f) * T4010A1_AMPLITUDE; + *r = propagate::(trans, attenuation, sound_speed, f); }); } diff --git a/autd3-gain-holo/src/linear_synthesis/gs.rs b/autd3-gain-holo/src/linear_synthesis/gs.rs index 66ae3152..9cc428e1 100644 --- a/autd3-gain-holo/src/linear_synthesis/gs.rs +++ b/autd3-gain-holo/src/linear_synthesis/gs.rs @@ -4,7 +4,7 @@ * Created Date: 29/05/2021 * Author: Shun Suzuki * ----- - * Last Modified: 23/11/2023 + * Last Modified: 09/01/2024 * Modified By: Shun Suzuki (suzuki@hapis.k.u-tokyo.ac.jp) * ----- * Copyright (c) 2021 Shun Suzuki. All rights reserved. @@ -20,7 +20,6 @@ use crate::{ use autd3_derive::Gain; use autd3_driver::{ - defined::T4010A1_AMPLITUDE, derive::prelude::*, geometry::{Geometry, Vector3}, }; @@ -70,19 +69,19 @@ impl Gain for GS { .backend .generate_propagation_matrix(geometry, &self.foci, &filter)?; - let m = self.backend.cols_c(&g)?; - let n = self.foci.len(); - let ones = vec![T4010A1_AMPLITUDE; m]; + let m = self.foci.len(); + let n = self.backend.cols_c(&g)?; + let ones = vec![1.; n]; - let mut b = self.backend.alloc_cm(m, n)?; - self.backend.gen_back_prop(m, n, &g, &mut b)?; + let mut b = self.backend.alloc_cm(n, m)?; + self.backend.gen_back_prop(n, m, &g, &mut b)?; let mut q = self.backend.from_slice_cv(&ones)?; let q0 = self.backend.from_slice_cv(&ones)?; let amps = self.backend.from_slice_cv(self.amps_as_slice())?; - let mut p = self.backend.alloc_zeros_cv(n)?; + let mut p = self.backend.alloc_zeros_cv(m)?; for _ in 0..self.repeat { self.backend.gemv_c( Trans::NoTrans, @@ -106,9 +105,6 @@ impl Gain for GS { self.backend.scaled_to_assign_cv(&q0, &mut q)?; } - self.backend - .scale_assign_cv(Complex::new(1. / T4010A1_AMPLITUDE, 0.0), &mut q)?; - generate_result( geometry, self.backend.to_host_cv(q)?, diff --git a/autd3-gain-holo/src/linear_synthesis/gspat.rs b/autd3-gain-holo/src/linear_synthesis/gspat.rs index 2f778056..114cd420 100644 --- a/autd3-gain-holo/src/linear_synthesis/gspat.rs +++ b/autd3-gain-holo/src/linear_synthesis/gspat.rs @@ -4,7 +4,7 @@ * Created Date: 29/05/2021 * Author: Shun Suzuki * ----- - * Last Modified: 23/11/2023 + * Last Modified: 09/01/2024 * Modified By: Shun Suzuki (suzuki@hapis.k.u-tokyo.ac.jp) * ----- * Copyright (c) 2021 Shun Suzuki. All rights reserved. @@ -20,7 +20,6 @@ use crate::{ use autd3_derive::Gain; use autd3_driver::{ - defined::T4010A1_AMPLITUDE, derive::prelude::*, geometry::{Geometry, Vector3}, }; @@ -70,17 +69,17 @@ impl Gain for GSPAT { .backend .generate_propagation_matrix(geometry, &self.foci, &filter)?; - let m = self.backend.cols_c(&g)?; - let n = self.foci.len(); + let m = self.foci.len(); + let n = self.backend.cols_c(&g)?; - let mut q = self.backend.alloc_zeros_cv(m)?; + let mut q = self.backend.alloc_zeros_cv(n)?; let amps = self.backend.from_slice_cv(self.amps_as_slice())?; - let mut b = self.backend.alloc_cm(m, n)?; - self.backend.gen_back_prop(m, n, &g, &mut b)?; + let mut b = self.backend.alloc_cm(n, m)?; + self.backend.gen_back_prop(n, m, &g, &mut b)?; - let mut r = self.backend.alloc_zeros_cm(n, n)?; + let mut r = self.backend.alloc_zeros_cm(m, m)?; self.backend.gemm_c( Trans::NoTrans, Trans::NoTrans, @@ -121,8 +120,6 @@ impl Gain for GSPAT { Complex::new(0., 0.), &mut q, )?; - self.backend - .scale_assign_cv(Complex::new(1. / T4010A1_AMPLITUDE, 0.0), &mut q)?; generate_result( geometry, diff --git a/autd3-gain-holo/src/linear_synthesis/naive.rs b/autd3-gain-holo/src/linear_synthesis/naive.rs index 5a2d1aba..85601ba5 100644 --- a/autd3-gain-holo/src/linear_synthesis/naive.rs +++ b/autd3-gain-holo/src/linear_synthesis/naive.rs @@ -4,7 +4,7 @@ * Created Date: 28/05/2021 * Author: Shun Suzuki * ----- - * Last Modified: 23/11/2023 + * Last Modified: 09/01/2024 * Modified By: Shun Suzuki (suzuki@hapis.k.u-tokyo.ac.jp) * ----- * Copyright (c) 2021 Shun Suzuki. All rights reserved. @@ -19,7 +19,6 @@ use crate::{ }; use autd3_driver::{ - defined::T4010A1_AMPLITUDE, derive::prelude::*, geometry::{Geometry, Vector3}, }; @@ -58,14 +57,14 @@ impl Gain for Naive { .backend .generate_propagation_matrix(geometry, &self.foci, &filter)?; - let m = self.backend.cols_c(&g)?; - let n = self.foci.len(); + let m = self.foci.len(); + let n = self.backend.cols_c(&g)?; - let mut b = self.backend.alloc_cm(m, n)?; - self.backend.gen_back_prop(m, n, &g, &mut b)?; + let mut b = self.backend.alloc_cm(n, m)?; + self.backend.gen_back_prop(n, m, &g, &mut b)?; let p = self.backend.from_slice_cv(self.amps_as_slice())?; - let mut q = self.backend.alloc_zeros_cv(m)?; + let mut q = self.backend.alloc_zeros_cv(n)?; self.backend.gemv_c( Trans::NoTrans, Complex::new(1., 0.), @@ -74,8 +73,6 @@ impl Gain for Naive { Complex::new(0., 0.), &mut q, )?; - self.backend - .scale_assign_cv(Complex::new(1.0 / T4010A1_AMPLITUDE, 0.), &mut q)?; generate_result( geometry, diff --git a/autd3-gain-holo/src/matrix/sdp.rs b/autd3-gain-holo/src/matrix/sdp.rs index 5edbcb1a..7a4afc38 100644 --- a/autd3-gain-holo/src/matrix/sdp.rs +++ b/autd3-gain-holo/src/matrix/sdp.rs @@ -4,7 +4,7 @@ * Created Date: 28/05/2021 * Author: Shun Suzuki * ----- - * Last Modified: 23/11/2023 + * Last Modified: 09/01/2024 * Modified By: Shun Suzuki (suzuki@hapis.k.u-tokyo.ac.jp) * ----- * Copyright (c) 2021 Shun Suzuki. All rights reserved. @@ -22,7 +22,6 @@ use crate::{ use autd3_derive::Gain; use autd3_driver::{ - defined::T4010A1_AMPLITUDE, derive::prelude::*, geometry::{Geometry, Vector3}, }; @@ -93,131 +92,140 @@ impl Gain for SDP { .backend .generate_propagation_matrix(geometry, &self.foci, &filter)?; - let m = self.backend.cols_c(&G)?; - let n = self.foci.len(); - - let zeros = self.backend.alloc_zeros_cv(n)?; - let ones = self.backend.from_slice_cv(&vec![1.; n])?; - - let mut a = self.backend.alloc_zeros_cm(n, n)?; - let amps = self.backend.from_slice_cv(self.amps_as_slice())?; - self.backend.create_diagonal_c(&s, &mut a)?; - - let mut pseudo_inv_G = self.backend.alloc_zeros_cm(m, n)?; - let mut u_ = self.backend.alloc_cm(n, n)?; - let mut s = self.backend.alloc_cm(m, n)?; - let mut vt = self.backend.alloc_cm(m, m)?; - let mut buf = self.backend.alloc_zeros_cm(m, n)?; - let G_ = self.backend.clone_cm(&G)?; - self.backend.pseudo_inverse_svd( - G_, - self.alpha, - &mut u_, - &mut s, - &mut vt, - &mut buf, - &mut pseudo_inv_G, - )?; + let m = self.foci.len(); + let n = self.backend.cols_c(&G)?; + + let zeros = self.backend.alloc_zeros_cv(m)?; + let ones = self.backend.from_slice_cv(&vec![1.; m])?; + + let P = { + let mut P = self.backend.alloc_zeros_cm(m, m)?; + let amps = self.backend.from_slice_cv(self.amps_as_slice())?; + self.backend.create_diagonal_c(&s, &mut P)?; + P + }; + + let G_inv = { + let mut G_inv = self.backend.alloc_zeros_cm(n, m)?; + let mut u_ = self.backend.alloc_cm(m, m)?; + let mut s = self.backend.alloc_cm(n, m)?; + let mut vt = self.backend.alloc_cm(n, n)?; + let mut buf = self.backend.alloc_zeros_cm(n, m)?; + let G_ = self.backend.clone_cm(&G)?; + self.backend.pseudo_inverse_svd( + G_, self.alpha, &mut u_, &mut s, &mut vt, &mut buf, &mut G_inv, + )?; + G_inv + }; let M = { - let mut M = self.backend.alloc_cm(n, n)?; + let mut M = self.backend.alloc_cm(m, m)?; + + // M = I self.backend.create_diagonal_c(&ones, &mut M)?; + // M = I - GG^{-1} self.backend.gemm_c( Trans::NoTrans, Trans::NoTrans, Complex::new(-1., 0.), &G, - &pseudo_inv_G, + &G_inv, Complex::new(1., 0.), &mut M, )?; - let mut tmp = self.backend.alloc_zeros_cm(n, n)?; + // tmp = P (I - GG^{-1}) + let mut tmp = self.backend.alloc_zeros_cm(m, m)?; self.backend.gemm_c( Trans::NoTrans, Trans::NoTrans, Complex::new(1., 0.), - &a, + &P, &M, Complex::new(0., 0.), &mut tmp, )?; + + // M = P (I - GG^{-1}) P self.backend.gemm_c( Trans::NoTrans, Trans::NoTrans, Complex::new(1., 0.), &tmp, - &a, + &P, Complex::new(0., 0.), &mut M, )?; M }; - let mut U = self.backend.alloc_cm(n, n)?; - self.backend.create_diagonal_c(&ones, &mut U)?; - - let mut rng = rand::thread_rng(); - - let mut x = self.backend.alloc_zeros_cv(n)?; - let mut Mc = self.backend.alloc_cv(n)?; - for _ in 0..self.repeat { - let i = rng.gen_range(0..n); - - self.backend.get_col_c(&M, i, &mut Mc)?; - self.backend.set_cv(i, Complex::new(0., 0.), &mut Mc)?; - - self.backend.gemv_c( - Trans::NoTrans, - Complex::new(1., 0.), - &U, - &Mc, - Complex::new(0., 0.), - &mut x, - )?; - - let gamma = self.backend.dot_c(&x, &Mc)?; - if gamma.re > 0. { - self.backend - .scale_assign_cv(Complex::new(-(self.lambda / gamma.re).sqrt(), 0.), &mut x)?; - - self.backend.set_col_c(&x, i, 0, i, &mut U)?; - self.backend.set_col_c(&x, i, i + 1, n, &mut U)?; - self.backend.conj_assign_v(&mut x)?; - self.backend.set_row_c(&x, i, 0, i, &mut U)?; - self.backend.set_row_c(&x, i, i + 1, n, &mut U)?; - } else { - self.backend.set_col_c(&zeros, i, 0, i, &mut U)?; - self.backend.set_col_c(&zeros, i, i + 1, n, &mut U)?; - self.backend.set_row_c(&zeros, i, 0, i, &mut U)?; - self.backend.set_row_c(&zeros, i, i + 1, n, &mut U)?; + // Block coordinate descent + let u = { + let mut U = self.backend.alloc_cm(m, m)?; + self.backend.create_diagonal_c(&ones, &mut U)?; + + let mut rng = rand::thread_rng(); + + let mut x = self.backend.alloc_zeros_cv(m)?; + let mut Mc = self.backend.alloc_cv(m)?; + for _ in 0..self.repeat { + let i = rng.gen_range(0..m); + + self.backend.get_col_c(&M, i, &mut Mc)?; + self.backend.set_cv(i, Complex::new(0., 0.), &mut Mc)?; + + self.backend.gemv_c( + Trans::NoTrans, + Complex::new(1., 0.), + &U, + &Mc, + Complex::new(0., 0.), + &mut x, + )?; + + let gamma = self.backend.dot_c(&x, &Mc)?; + if gamma.re > 0. { + self.backend.scale_assign_cv( + Complex::new(-(self.lambda / gamma.re).sqrt(), 0.), + &mut x, + )?; + + self.backend.set_col_c(&x, i, 0, i, &mut U)?; + self.backend.set_col_c(&x, i, i + 1, m, &mut U)?; + self.backend.conj_assign_v(&mut x)?; + self.backend.set_row_c(&x, i, 0, i, &mut U)?; + self.backend.set_row_c(&x, i, i + 1, m, &mut U)?; + } else { + self.backend.set_col_c(&zeros, i, 0, i, &mut U)?; + self.backend.set_col_c(&zeros, i, i + 1, m, &mut U)?; + self.backend.set_row_c(&zeros, i, 0, i, &mut U)?; + self.backend.set_row_c(&zeros, i, i + 1, m, &mut U)?; + } } - } - let u = self.backend.max_eigen_vector_c(U)?; + self.backend.max_eigen_vector_c(U)? + }; - let mut ut = self.backend.alloc_zeros_cv(n)?; + let mut ut = self.backend.alloc_zeros_cv(m)?; self.backend.gemv_c( Trans::NoTrans, Complex::new(1., 0.), - &a, + &P, &u, Complex::new(0., 0.), &mut ut, )?; - let mut q = self.backend.alloc_zeros_cv(m)?; + let mut q = self.backend.alloc_zeros_cv(n)?; self.backend.gemv_c( Trans::NoTrans, Complex::new(1., 0.), - &pseudo_inv_G, + &G_inv, &ut, Complex::new(0., 0.), &mut q, )?; - self.backend - .scale_assign_cv(Complex::new(1. / T4010A1_AMPLITUDE, 0.0), &mut q)?; generate_result( geometry, diff --git a/autd3-gain-holo/src/nls/lm.rs b/autd3-gain-holo/src/nls/lm.rs index 3dcfec7a..d924ed3b 100644 --- a/autd3-gain-holo/src/nls/lm.rs +++ b/autd3-gain-holo/src/nls/lm.rs @@ -4,7 +4,7 @@ * Created Date: 29/05/2021 * Author: Shun Suzuki * ----- - * Last Modified: 02/12/2023 + * Last Modified: 09/01/2024 * Modified By: Shun Suzuki (suzuki@hapis.k.u-tokyo.ac.jp) * ----- * Copyright (c) 2021 Shun Suzuki. All rights reserved. @@ -20,7 +20,7 @@ use autd3_derive::Gain; use autd3_driver::{ common::{EmitIntensity, Phase}, - defined::{PI, T4010A1_AMPLITUDE}, + defined::PI, derive::prelude::*, geometry::{Geometry, Vector3}, }; @@ -172,28 +172,26 @@ impl Gain for LM { geometry: &Geometry, filter: GainFilter, ) -> Result>, AUTDInternalError> { - let mut g = self + let g = self .backend .generate_propagation_matrix(geometry, &self.foci, &filter)?; - self.backend - .scale_assign_cm(Complex::new(T4010A1_AMPLITUDE, 0.), &mut g)?; - let m = self.backend.cols_c(&g)?; - let n = self.foci.len(); + let n = self.backend.cols_c(&g)?; + let m = self.foci.len(); - let n_param = n + m; + let n_param = m + n; let bhb = { let mut bhb = self.backend.alloc_zeros_cm(n_param, n_param)?; let mut amps = self.backend.from_slice_cv(self.amps_as_slice())?; - let mut p = self.backend.alloc_cm(n, n)?; + let mut p = self.backend.alloc_cm(m, m)?; self.backend .scale_assign_cv(Complex::new(-1., 0.), &mut amps)?; self.backend.create_diagonal_c(&s, &mut p)?; - let mut b = self.backend.alloc_cm(n, n_param)?; + let mut b = self.backend.alloc_cm(m, n_param)?; self.backend.concat_col_cm(&g, &p, &mut b)?; self.backend.gemm_c( diff --git a/examples/src/tests/holo.rs b/examples/src/tests/holo.rs index e20389b2..6b6367cd 100644 --- a/examples/src/tests/holo.rs +++ b/examples/src/tests/holo.rs @@ -4,7 +4,7 @@ * Created Date: 29/05/2021 * Author: Shun Suzuki * ----- - * Last Modified: 26/12/2023 + * Last Modified: 09/01/2024 * Modified By: Shun Suzuki (suzuki@hapis.k.u-tokyo.ac.jp) * ----- * Copyright (c) 2021 Shun Suzuki. All rights reserved. @@ -40,7 +40,7 @@ pub async fn holo(autd: &mut Controller) -> anyhow::Result { let backend = NalgebraBackend::new()?; - let target_amp = 5e3 * autd.geometry.num_devices() as float * Pascal; + let target_amp = 2.5e3 * autd.geometry.num_devices() as float * Pascal; match s.trim().parse::() { Ok(0) => { let g = SDP::new(backend)