diff --git a/Cargo.lock b/Cargo.lock index a7bad1b..3afd84f 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -252,6 +252,12 @@ version = "0.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2da6da31387c7e4ef160ffab6d5e7f00c42626fe39aea70a7b0f1773f7dd6c1b" +[[package]] +name = "closure" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d6173fd61b610d15a7566dd7b7620775627441c4ab9dac8906e17cb93a24b782" + [[package]] name = "colorchoice" version = "1.0.1" @@ -875,6 +881,7 @@ dependencies = [ "assert", "blas", "blas-sys", + "closure", "colored", "criterion", "csv", diff --git a/Cargo.toml b/Cargo.toml index bfb75ff..ade86d5 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -28,6 +28,7 @@ intel-mkl-src = {version = "0.8.1", features = ["mkl-dynamic-lp64-iomp"]} #openblas-src = "0.10.8" # For python interface pyo3 = {version = "0.21.2", features = ["extension-module"], optional = true} +# Parallelism rayon = "1.10.0" # Mersenne Twister rand_mt = {version = "4.2.2", features = ["rand_core"]} @@ -35,6 +36,8 @@ rand_mt = {version = "4.2.2", features = ["rand_core"]} indicatif = "0.17.8" # Colored terminal text colored = "3.0.0" +# Selectively move into closure +closure = "0.3.0" [dependencies.rand] version = "0.8.5" diff --git a/src/bin/boboche2000.rs b/src/bin/boboche2000.rs index 14d1995..853e513 100644 --- a/src/bin/boboche2000.rs +++ b/src/bin/boboche2000.rs @@ -84,26 +84,28 @@ fn main() { hopping_bitmask: &bitmask, clean_update_frequency: CLEAN_UPDATE_FREQUENCY, nmcwarmup: NMCWARMUP, + nwarmupchains: 1, nmcsample: NMCSAMP, tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, tolerance_singularity: TOLERANCE_SINGULARITY, pair_wavefunction: false, + _opt_iter: 0, }; log_system_parameters(&sys); - let mut otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; - let mut expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut visited: Vec = vec![0; NMCSAMP]; + let otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited: Vec = vec![0; NMCSAMP]; let mut der = DerivativeOperator { - o_tilde: &mut otilde, - expval_o: &mut expvalo, - ho: &mut expval_ho, + o_tilde: otilde.into_boxed_slice(), + expval_o: expvalo.into_boxed_slice(), + ho: expval_ho.into_boxed_slice(), n: (NFIJ + NVIJ + NGI) as i32, nsamp: NMCSAMP as f64, nsamp_int: 1, mu: -1, - visited: &mut visited, + visited: visited.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: 0.0, diff --git a/src/bin/dvmc_fastupdate.rs b/src/bin/dvmc_fastupdate.rs index f979ab5..5c29566 100644 --- a/src/bin/dvmc_fastupdate.rs +++ b/src/bin/dvmc_fastupdate.rs @@ -19,6 +19,7 @@ const NELEC: usize = SIZE; const NMCSAMP: usize = 10_000; const NBOOTSTRAP: usize = 1; const NMCWARMUP: usize = 1000; +const NWARMUPCHAINS: usize = 1; const MCSAMPLE_INTERVAL: usize = 1; const _NTHREADS: usize = 6; const CLEAN_UPDATE_FREQUENCY: usize = 2; @@ -160,10 +161,12 @@ fn main() { nmcsample: NMCSAMP, nbootstrap: NBOOTSTRAP, nmcwarmup: NMCWARMUP, + nwarmupchains: NWARMUPCHAINS, mcsample_interval: MCSAMPLE_INTERVAL, tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, tolerance_singularity: TOLERENCE_SINGULARITY, pair_wavefunction: false, + _opt_iter: 0, }; let mut rng = Mt64::new(SEED); @@ -270,19 +273,19 @@ fn main() { tmp }; - let mut otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; - let mut expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut visited: Vec = vec![0; NMCSAMP]; + let otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited: Vec = vec![0; NMCSAMP]; let mut derivative = DerivativeOperator { - o_tilde: &mut otilde, - expval_o: &mut expvalo, - ho: &mut expval_ho, + o_tilde: otilde.into_boxed_slice(), + expval_o: expvalo.into_boxed_slice(), + ho: expval_ho.into_boxed_slice(), n: (NFIJ + NVIJ + NGI) as i32, nsamp: NMCSAMP as f64, nsamp_int: MCSAMPLE_INTERVAL, mu: -1, - visited: &mut visited, + visited: visited.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: EPSILON_SPREAD, @@ -320,8 +323,8 @@ fn main() { unsafe { let incx = 1; let incy = 1; - daxpy(derivative.n, -mean_energy, derivative.expval_o, incx, derivative.ho, incy); - dcopy(derivative.n, derivative.ho, incx, &mut b, incy); + daxpy(derivative.n, -mean_energy, &derivative.expval_o, incx, &mut derivative.ho, incy); + dcopy(derivative.n, &derivative.ho, incx, &mut b, incy); } spread_eigenvalues(&mut derivative); exact_overlap_inverse(&derivative, &mut b, EPSILON_CG, NPARAMS as i32, PARAMTHRESHOLD); diff --git a/src/bin/dvmc_pairwf.rs b/src/bin/dvmc_pairwf.rs index 55df29f..0f5bf6b 100644 --- a/src/bin/dvmc_pairwf.rs +++ b/src/bin/dvmc_pairwf.rs @@ -19,6 +19,7 @@ const NELEC: usize = SIZE; const NMCSAMP: usize = 500; const NBOOTSTRAP: usize = 1; const NMCWARMUP: usize = 1000; +const NWARMUPCHAINS: usize = 1; const MCSAMPLE_INTERVAL: usize = 1; const _NTHREADS: usize = 1; const CLEAN_UPDATE_FREQUENCY: usize = 32; @@ -87,7 +88,7 @@ fn _save_otilde(fp: &mut File, der: &DerivativeOperator) { let mut c = vec![0.0; (der.n * der.n) as usize]; println!("dim = {}", der.n * der.n); unsafe { - dgemm(b"N"[0], b"T"[0], der.n, der.n, der.mu, 1.0, der.o_tilde, der.n, der.o_tilde, der.n, 0.0, &mut c, der.n); + dgemm(b"N"[0], b"T"[0], der.n, der.n, der.mu, 1.0, &der.o_tilde, der.n, &der.o_tilde, der.n, 0.0, &mut c, der.n); } let mut outstr = "".to_owned(); outstr.push_str(&format!(" = ")); @@ -192,7 +193,7 @@ fn main() { // Initialize logger env_logger::init(); let bitmask = generate_bitmask(&HOPPINGS, SIZE); - let system_params = SysParams { + let mut system_params = SysParams { size: SIZE, nelec: NELEC, array_size: (SIZE + 7) / 8, @@ -207,10 +208,12 @@ fn main() { nmcsample: NMCSAMP, nbootstrap: NBOOTSTRAP, nmcwarmup: NMCWARMUP, + nwarmupchains: NWARMUPCHAINS, mcsample_interval: MCSAMPLE_INTERVAL, tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, tolerance_singularity: TOLERENCE_SINGULARITY, pair_wavefunction: true, + _opt_iter: 0, }; let mut rng = Mt64::new(SEED); @@ -253,18 +256,18 @@ fn main() { tmp }; - let mut otilde: Vec = vec![0.0; (4*NFIJ + NVIJ + NGI) * (NMCSAMP + 1)]; - let mut work_otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * (NMCSAMP + 1)]; - let mut expvalo: Vec = vec![0.0; 4*NFIJ + NVIJ + NGI]; - let mut work_expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut expval_ho: Vec = vec![0.0; 4*NFIJ + NVIJ + NGI]; - let mut work_expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut visited: Vec = vec![0; NMCSAMP + 1]; - let mut work_visited: Vec = vec![0; NMCSAMP + 1]; + let otilde: Vec = vec![0.0; (4*NFIJ + NVIJ + NGI) * (NMCSAMP + 1)]; + let work_otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * (NMCSAMP + 1)]; + let expvalo: Vec = vec![0.0; 4*NFIJ + NVIJ + NGI]; + let work_expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho: Vec = vec![0.0; 4*NFIJ + NVIJ + NGI]; + let work_expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited: Vec = vec![0; NMCSAMP + 1]; + let work_visited: Vec = vec![0; NMCSAMP + 1]; let mut derivative = DerivativeOperator { - o_tilde: &mut otilde, - expval_o: &mut expvalo, - ho: &mut expval_ho, + o_tilde: otilde.into_boxed_slice(), + expval_o: expvalo.into_boxed_slice(), + ho: expval_ho.into_boxed_slice(), n: (4*NFIJ + NVIJ + NGI) as i32, nsamp: match COMPUTE_ENERGY_METHOD { EnergyComputationMethod::ExactSum => 1.0, @@ -272,15 +275,15 @@ fn main() { }, nsamp_int: MCSAMPLE_INTERVAL, mu: -1, - visited: &mut visited, + visited: visited.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: EPSILON_SHIFT, }; let mut work_derivative = DerivativeOperator { - o_tilde: &mut work_otilde, - expval_o: &mut work_expvalo, - ho: &mut work_expval_ho, + o_tilde: work_otilde.into_boxed_slice(), + expval_o: work_expvalo.into_boxed_slice(), + ho: work_expval_ho.into_boxed_slice(), n: (NFIJ + NVIJ + NGI) as i32, nsamp: match COMPUTE_ENERGY_METHOD { EnergyComputationMethod::ExactSum => 1.0, @@ -288,7 +291,7 @@ fn main() { }, nsamp_int: MCSAMPLE_INTERVAL, mu: -1, - visited: &mut work_visited, + visited: work_visited.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: EPSILON_SHIFT, @@ -305,6 +308,7 @@ fn main() { .progress_chars("##-")); for opt_iter in 0..NOPTITER { + system_params._opt_iter = opt_iter; let (mean_energy, accumulated_states, deltae, correlation_time) = { match COMPUTE_ENERGY_METHOD { EnergyComputationMethod::MonteCarlo => compute_mean_energy(&mut rng, state, ¶meters, &system_params, &mut derivative), @@ -347,8 +351,8 @@ fn main() { unsafe { let incx = 1; let incy = 1; - daxpy(work_derivative.n, -mean_energy, work_derivative.expval_o, incx, work_derivative.ho, incy); - dcopy(work_derivative.n, work_derivative.ho, incx, &mut b, incy); + daxpy(work_derivative.n, -mean_energy, &work_derivative.expval_o, incx, &mut work_derivative.ho, incy); + dcopy(work_derivative.n, &work_derivative.ho, incx, &mut b, incy); } //save_otilde(&mut der_fp, &derivative); //save_otilde(&mut wder_fp, &work_derivative); @@ -463,5 +467,5 @@ fn main() { opt_progress_bar.inc(1); opt_progress_bar.set_message(format!("Changed parameters by norm: {:+>.05e} Current energy: {:+>.05e}", opt_delta, mean_energy)); } - opt_progress_bar.finish() + opt_progress_bar.finish(); } diff --git a/src/bin/dvmc_pairwf2.rs b/src/bin/dvmc_pairwf2.rs index 426c69e..955be9f 100644 --- a/src/bin/dvmc_pairwf2.rs +++ b/src/bin/dvmc_pairwf2.rs @@ -1,7 +1,11 @@ +use impurity::optimisation::GenParameterMap; //use blas::dcopy; //use log::{debug, info}; -use rand_mt::Mt64; +use rand_mt::{Mt19937GenRand64, Mt64}; use rand::Rng; +use std::mem; + +use rayon::ThreadPoolBuilder; //use indicatif::{ProgressBar, ProgressStyle}; use std::fs::File; use std::io::Write; @@ -9,33 +13,36 @@ use std::io::Write; use impurity::{VarParams, SysParams, generate_bitmask, FockState, RandomStateGeneration}; use impurity::dvmc::{variationnal_monte_carlo, EnergyOptimisationMethod, EnergyComputationMethod, VMCParams}; +type BitSize = u32; + const SEED: u64 = 14; -const SIZE: usize = 8; +const SIZE: usize = 2; const NFIJ: usize = 4*SIZE*SIZE; const NVIJ: usize = SIZE*(SIZE - 1) / 2; const NGI: usize = SIZE; const NPARAMS: usize = NFIJ + NGI + NVIJ; const NELEC: usize = SIZE; -const NMCSAMP: usize = 1_000; +const NMCSAMP: usize = 200; const NBOOTSTRAP: usize = 1; -const NMCWARMUP: usize = 1000; +const NMCWARMUP: usize = 100; +const NWARMUPCHAINS: usize = 1; const MCSAMPLE_INTERVAL: usize = 1; -const _NTHREADS: usize = 1; +const NTHREADS: usize = 12; const CLEAN_UPDATE_FREQUENCY: usize = 32; const TOLERENCE_SHERMAN_MORRISSON: f64 = 1e-12; const TOLERENCE_SINGULARITY: f64 = 1e-12; const CONS_U: f64 = 1.0; const CONS_T: f64 = 1.0; -const INITIAL_RATIO_UT: f64 = 1.0; +const INITIAL_RATIO_UT: f64 = 8.0; const FINAL_RATIO_UT: f64 = 32.0; const NRATIO_POINTS: usize = 1; const EPSILON_CG: f64 = 1e-16; -const EPSILON_SHIFT: f64 = 1e-2; +const EPSILON_SHIFT: f64 = 1e-3; const OPTIMISATION_TIME_STEP: f64 = 1e-2; const OPTIMISATION_DECAY: f64 = 0.0; const NOPTITER: usize = 1000; const KMAX: usize = NPARAMS; -const PARAM_THRESHOLD: f64 = ::MIN; +const PARAM_THRESHOLD: f64 = 1e-3; //const PARAM_THRESHOLD: f64 = 0.0; const OPTIMISE: bool = true; const OPTIMISE_GUTZ: bool = true; @@ -45,15 +52,20 @@ const SET_EXPVALO_ZERO: bool = false; const COMPUTE_ENERGY_METHOD: EnergyComputationMethod = EnergyComputationMethod::MonteCarlo; const OPTIMISE_ENERGY_METHOD: EnergyOptimisationMethod = EnergyOptimisationMethod::ConjugateGradiant; const ENERGY_CONV_AVERAGE_SAMPLE: usize = 20; +const N_INDEP_PARAMS: usize = SIZE*SIZE + NGI + NVIJ; +const N_GUTZ: usize = NGI; +const N_JAST: usize = NVIJ; +const PAIRWF: bool = false; +const CONV_PARAM_THRESHOLD: f64 = 1e-10; // 4 et 2 Sites -//pub const HOPPINGS: [f64; SIZE*SIZE] = [ -// //0.0, 1.0, 1.0, 0.0, -// 0.0, 1.0, 1.0, 0.0, -// 1.0, 0.0, 0.0, 1.0, -// 1.0, 0.0, 0.0, 1.0, -// 0.0, 1.0, 1.0, 0.0 -//]; +pub const HOPPINGS: [f64; SIZE*SIZE] = [ + 0.0, 1.0, 1.0, 0.0, + //0.0, 1.0, 1.0, 0.0, + //1.0, 0.0, 0.0, 1.0, + //1.0, 0.0, 0.0, 1.0, + //0.0, 1.0, 1.0, 0.0 +]; // 7 Sites //pub const HOPPINGS: [f64; SIZE*SIZE] = [ @@ -67,16 +79,43 @@ const ENERGY_CONV_AVERAGE_SAMPLE: usize = 20; //]; // 8 Sites -pub const HOPPINGS: [f64; SIZE*SIZE] = [ - 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, - 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, - 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, - 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, - 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, - 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, -]; +//pub const HOPPINGS: [f64; SIZE*SIZE] = [ +// 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, +// 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, +// 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, +// 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, +// 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, +// 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, +// 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, +// 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, +//]; + +// 9 Sites +//pub const HOPPINGS: [f64; SIZE*SIZE] = [ +// 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, +// 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, +// 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, +// 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, +// 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, +// 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, +// 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, +// 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +//]; + +// 10 Sites +//pub const HOPPINGS: [f64; SIZE*SIZE] = [ +// 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, +// 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, +// 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, +// 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +//]; // 16 Sites //pub const HOPPINGS: [f64; SIZE*SIZE] = [ @@ -97,6 +136,71 @@ pub const HOPPINGS: [f64; SIZE*SIZE] = [ // 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, // 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, //]; +// + +// Pairwf +//const PARAMS_PROJECTOR: [f64; (NFIJ + NVIJ + NGI) * (NFIJ + NVIJ + NGI - 1) / 2 + NFIJ + NVIJ + NGI] = [ +// /* g0 */ 1.0, // 0, +// /* g1 */ 1.0, 0.0, // +// /* v01 */ 0.0, 0.0, 0.0, +// /* f00uu */ 0.0, 0.0, 0.0, 0.0, +// /* f01uu */ 0.0, 0.0, 0.0, 0.0, 0.0, +// /* f10uu */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* f11uu */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* f00ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, +// /* f01ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, +// /* f10ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, +// /* f11ud */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, +// /* f00du */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* f01du */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* f10du */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* f11du */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* f00dd */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* f01dd */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* f10dd */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +// /* f11dd */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, +//]; + +// General rep +//const PARAMS_PROJECTOR: [f64; (NFIJ + NVIJ + NGI) * (NFIJ + NVIJ + NGI - 1) / 2 + NFIJ + NVIJ + NGI] = { +// let mut param = [0.0; (NFIJ + NVIJ + NGI) * (NFIJ + NVIJ + NGI - 1) / 2 + NFIJ + NVIJ + NGI]; +// let mut j = 0; +// let mut n = 0; +// while j < NFIJ + NVIJ + NGI { +// param[j + (j * (j+1) / 2)] = 1.0; +// j += 1; +// n += 1; +// } +// if n != N_INDEP_PARAMS { +// panic!("Number of set independant params is not correct."); +// } +// param +//}; + + +// General pairwf rep +const PARAMS_PROJECTOR: [f64; (NFIJ + NVIJ + NGI) * (NFIJ + NVIJ + NGI - 1) / 2 + NFIJ + NVIJ + NGI] = { + let mut param = [0.0; (NFIJ + NVIJ + NGI) * (NFIJ + NVIJ + NGI - 1) / 2 + NFIJ + NVIJ + NGI]; + let mut j = 0; + let mut n = 0; + while j < NVIJ + NGI { + param[j + (j * (j+1) / 2)] = 1.0; + j += 1; + n += 1; + } + let mut j = NVIJ + NGI; + while j < SIZE * SIZE + NVIJ + NGI{ + j += SIZE * SIZE; + param[j + (j * (j+1) / 2)] = 1.0; + j -= SIZE * SIZE; + j += 1; + n += 1; + } + if n != N_INDEP_PARAMS { + panic!("Number of set independant params is not correct."); + } + param +}; fn sq(a: f64) -> f64 { @@ -127,6 +231,7 @@ fn log_energy_convs(e: &[f64], fp: &mut File) { fn main() { let mut fp = File::create("u_t_sweep").unwrap(); + let pool = ThreadPoolBuilder::new().num_threads(NTHREADS).build_global().unwrap(); writeln!(fp, "{}", format!("# {} {} {}", SIZE, NMCSAMP, NOPTITER)).unwrap(); let mut paramsfp = File::create("params").unwrap(); writeln!(paramsfp, "{}", format!("# {} {} {}", SIZE, NMCSAMP, NOPTITER)).unwrap(); @@ -134,10 +239,14 @@ fn main() { // Initialize logger env_logger::init(); let bitmask = generate_bitmask(&HOPPINGS, SIZE); - let mut rng = Mt64::new(SEED); + let mut rng = Vec::new(); + for i in 0..NTHREADS { + rng.push(Mt64::new(SEED + i as u64)); + } + let mut rngs: Vec<&mut Mt64> = rng.iter_mut().collect(); for nsweep_iter in 0..NRATIO_POINTS { - let system_params = SysParams { + let mut system_params = SysParams { size: SIZE, nelec: NELEC, array_size: (SIZE + 7) / 8, @@ -152,16 +261,18 @@ fn main() { nmcsample: NMCSAMP, nbootstrap: NBOOTSTRAP, nmcwarmup: NMCWARMUP, + nwarmupchains: NWARMUPCHAINS, mcsample_interval: MCSAMPLE_INTERVAL, tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, tolerance_singularity: TOLERENCE_SINGULARITY, - pair_wavefunction: true, + pair_wavefunction: PAIRWF, + _opt_iter: 0, }; - println!("U = {}, T = {}", system_params.cons_u, system_params.cons_t); + //println!("U = {}, T = {}", system_params.cons_u, system_params.cons_t); let mut all_params: Vec = Vec::with_capacity(NGI + NVIJ + NFIJ); for _ in 0..(NGI + NVIJ + NFIJ) { - all_params.push(rng.gen()); + all_params.push(rngs[0].gen()); } let (gi, params) = all_params.split_at_mut(NGI); let (vij, fij) = params.split_at_mut(NVIJ); @@ -171,6 +282,13 @@ fn main() { vij, size: SIZE }; + //for i in 0..SIZE*SIZE { + // parameters.fij[i] = 0.0; + // parameters.fij[i +SIZE*SIZE] = 0.5; + // parameters.fij[i +2*SIZE*SIZE] = 0.0; + // parameters.fij[i +3*SIZE*SIZE] = 0.0; + //} + //println!("{:?}", parameters.fij); let vmcparams = VMCParams { dt: OPTIMISATION_TIME_STEP, @@ -180,27 +298,42 @@ fn main() { epsilon: EPSILON_SHIFT, epsilon_cg: EPSILON_CG, noptiter: NOPTITER, - nparams: NGI + NVIJ + NFIJ, + nparams: N_INDEP_PARAMS, optimise: OPTIMISE, optimise_gutzwiller: OPTIMISE_GUTZ, optimise_jastrow: OPTIMISE_JAST, optimise_orbital: OPTIMISE_ORB, compute_energy_method: COMPUTE_ENERGY_METHOD, optimise_energy_method: OPTIMISE_ENERGY_METHOD, + conv_param_threshold: CONV_PARAM_THRESHOLD, + nthreads: NTHREADS }; - let state: FockState = { - let mut tmp: FockState = FockState::generate_from_nelec(&mut rng, NELEC, SIZE); + let state: FockState = { + let mut tmp: FockState = FockState::generate_from_nelec(&mut rngs[0], NELEC, SIZE); while tmp.spin_up.count_ones() != tmp.spin_down.count_ones() { - tmp = FockState::generate_from_nelec(&mut rng, NELEC, SIZE); + tmp = FockState::generate_from_nelec(&mut rngs[0], NELEC, SIZE); } tmp }; + let mut states_vec = Vec::with_capacity(NTHREADS); + for i in 0..NTHREADS { + states_vec.push(state); + } + let param_map = GenParameterMap { + dim: N_INDEP_PARAMS as i32, + gendim: (NFIJ + NGI + NVIJ) as i32, + n_genparams: (NFIJ + NGI + NVIJ) as i32, + n_independant_gutzwiller: N_GUTZ, + n_independant_jastrow: N_JAST, + projector: Box::new(PARAMS_PROJECTOR), + }; - let e_array = variationnal_monte_carlo(&mut rng, state, &mut parameters, &system_params, &vmcparams); - write_energy(&mut fp, &e_array); + let e_array = variationnal_monte_carlo(&mut rngs, &mut states_vec, &mut parameters, &mut system_params, &vmcparams, ¶m_map); + //write_energy(&mut fp, &e_array); log_energy_convs(&e_array, &mut paramsfp); } + mem::drop(rng); } diff --git a/src/constants.rs b/src/constants.rs index 97fece5..6acf1d2 100644 --- a/src/constants.rs +++ b/src/constants.rs @@ -1,4 +1,4 @@ -const SIZE: usize = 17; +const SIZE: usize = 32; const ARRAY_SIZE: usize = (SIZE + 7) / 8; /// System parameters /// TODOC @@ -16,12 +16,61 @@ pub struct SysParams<'a> { pub hopping_bitmask: &'a [SpinState], pub nmcwarmup: usize, pub nmcsample: usize, + pub nwarmupchains: usize, pub nbootstrap: usize, pub mcsample_interval: usize, pub clean_update_frequency: usize, pub tolerance_singularity: f64, pub tolerance_sherman_morrison: f64, pub pair_wavefunction: bool, + pub _opt_iter: usize, +} + +impl<'a> SysParams<'a> { + pub fn new( + size: usize, + nelec: usize, + array_size: usize, + cons_u: f64, + cons_t: f64, + nfij: usize, + nvij: usize, + ngi: usize, + transfert_matrix: &'a [f64], + hopping_bitmask: &'a [SpinState], + nmcwarmup: usize, + nmcsample: usize, + nwarmupchains: usize, + nbootstrap: usize, + mcsample_interval: usize, + clean_update_frequency: usize, + tolerance_singularity: f64, + tolerance_sherman_morrison: f64, + pair_wavefunction: bool, + ) -> Self { + SysParams { + size, + nelec, + array_size, + cons_t, + cons_u, + nfij, + nvij, + ngi, + transfert_matrix, + hopping_bitmask, + nmcwarmup, + nmcsample, + nwarmupchains, + nbootstrap, + mcsample_interval, + clean_update_frequency, + tolerance_sherman_morrison, + tolerance_singularity, + pair_wavefunction, + _opt_iter: 0, + } + } } pub type BitStruct = u8; diff --git a/src/density.rs b/src/density.rs index 3f5587e..6b3e911 100644 --- a/src/density.rs +++ b/src/density.rs @@ -12,7 +12,7 @@ pub fn compute_internal_product( params: &VarParams, sys: &SysParams, ) -> f64 -where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl +where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl + Send { let (mut pfaffian_state, jastrow_exp, gutz_exp) = { ( @@ -34,7 +34,7 @@ pub fn compute_internal_product_parts( params: &VarParams, sys: &SysParams, ) -> (PfaffianState, f64) -where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl +where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl + Send { let (mut pfaffian_state, jastrow_exp, gutz_exp) = { ( @@ -55,7 +55,7 @@ pub fn fast_internal_product_no_otilde( previous_proj: &mut f64, params: &VarParams, ) -> (f64, Vec, usize) -where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl +where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl + Send { // Rename things. let previous_i = hopping.0; @@ -94,7 +94,7 @@ pub fn fast_internal_product( previous_proj: &mut f64, params: &VarParams, ) -> (f64, Vec, usize) -where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl +where T: BitOps + std::fmt::Debug + std::fmt::Display + From + std::ops::Shl + Send { // Rename things. let previous_i = hopping.0; diff --git a/src/dvmc.rs b/src/dvmc.rs index f28d43d..424c0e8 100644 --- a/src/dvmc.rs +++ b/src/dvmc.rs @@ -1,12 +1,15 @@ use log::{info, error}; use rand::distributions::{Distribution, Standard}; use rand::Rng; +use rayon::iter::PanicFuse; use std::fmt::{Display, Debug}; -use blas::{idamax, daxpy, dcopy, dscal}; +use blas::{idamax, daxpy, dcopy, dscal, dnrm2}; +use rayon::ThreadPool; +use std::thread; use crate::{VarParams, DerivativeOperator, SysParams, FockState, BitOps}; use crate::monte_carlo::{compute_mean_energy, compute_mean_energy_exact}; -use crate::optimisation::{conjugate_gradiant, exact_overlap_inverse}; +use crate::optimisation::{conjugate_gradiant, exact_overlap_inverse, GenParameterMap, ReducibleGeneralRepresentation}; #[derive(Debug)] pub enum EnergyOptimisationMethod { @@ -31,18 +34,56 @@ pub struct VMCParams { pub dt: f64, pub optimisation_decay: f64, pub nparams: usize, + pub nthreads: usize, pub threshold: f64, pub optimise: bool, pub optimise_gutzwiller: bool, pub optimise_jastrow: bool, pub optimise_orbital: bool, + pub conv_param_threshold: f64, +} + +fn merge_derivatives(der_vec: &[DerivativeOperator], nmc: usize, nthreads: usize) -> DerivativeOperator { + let mut out_der = DerivativeOperator::new( + der_vec[0].n, + 0, + nmc as f64, + der_vec[0].nsamp_int, + der_vec[0].pfaff_off, + der_vec[0].jas_off, + der_vec[0].epsilon + ); + let mut k = 0; + for i in 0..nthreads { + out_der.mu += der_vec[i].mu; + unsafe { + let incx = 1; + let incy = 1; + dcopy( + der_vec[i].n * der_vec[i].mu, + &der_vec[i].o_tilde, + incx, + &mut out_der.o_tilde[(der_vec[i].n * der_vec[i].mu) as usize * i..der_vec[i].n as usize * nmc], + incy + ); + daxpy(der_vec[i].n, 1.0 / nthreads as f64, &der_vec[i].expval_o, incx, &mut out_der.expval_o, incy); + daxpy(der_vec[i].n, 1.0 / nthreads as f64, &der_vec[i].ho, incx, &mut out_der.ho, incy); + } + for j in 0..der_vec[i].mu as usize { + out_der.visited[k] = der_vec[i].visited[j]; + k += 1; + + } + } + + out_der } fn zero_out_derivatives(der: &mut DerivativeOperator, sys: &SysParams) { - for i in 0.. (sys.nfij + sys.nvij + sys.ngi) * sys.nmcsample { + for i in 0.. (der.n as usize) * sys.nmcsample { der.o_tilde[i] = 0.0; } - for i in 0..sys.nfij + sys.nvij + sys.ngi { + for i in 0..der.n as usize { der.expval_o[i] = 0.0; der.ho[i] = 0.0; } @@ -52,40 +93,124 @@ fn zero_out_derivatives(der: &mut DerivativeOperator, sys: &SysParams) { der.mu = -1; } -pub fn variationnal_monte_carlo( - rng: &mut R, - initial_state: FockState, +fn update_initial_state + Display + Debug + Send + Sync>( + states: &mut [FockState], + acc_states: &[FockState], + nmcsample: usize, + nthreads: usize +) { + for i in 0..nthreads { + states[i] = acc_states[i * nmcsample]; + } +} + +fn parallel_monte_carlo( + rngs: &mut [&mut R], + initial_state: &[FockState], + work_der_vec: &mut [DerivativeOperator], + vmcparams: &VMCParams, + params: & VarParams, + sys: & SysParams +) -> (f64, Vec>, f64, f64) +where + T: BitOps + From + Display + Debug + Send + Sync, + R: Rng + ?Sized + Send + Sync, + Standard: Distribution + Send +{ + let mut res_vec = Vec::new(); + // Scoped thread because threads need to accept not static stack parameters + thread::scope(|scope| { + let threads: Vec<_> = (0..vmcparams.nthreads) + .map(|idx| { + // Wrestle the borrow checker + let wder_ptr = &mut work_der_vec[idx] as *mut _; + let rng_ptr = &mut *rngs[idx] as *mut R; + let state_ptr = &initial_state[idx] as *const _; + let rng = unsafe { &mut *rng_ptr}; + let wder = unsafe { &mut *wder_ptr}; + let state = unsafe { & *state_ptr}; + scope.spawn( + || { + compute_mean_energy( + rng, + *state, + params, + sys, + wder + ) + } + ) + }) + .collect(); + + for handle in threads { + res_vec.push(handle.join().unwrap()); + } + }); + let mut out_energy = 0.0; + let mut out_states_vec = Vec::new(); + let mut out_de = 0.0; + let mut out_corrtime = 0.0; + for thread in 0..vmcparams.nthreads { + //println!("Energy of thread {} = {}", thread, res_vec[thread].0); + out_energy += res_vec[thread].0; + out_de += res_vec[thread].2; + out_corrtime += res_vec[thread].3; + out_states_vec.append(&mut res_vec[thread].1); + } + + //panic!("Stop!"); + //(res_vec[0].0, res_vec[0].1.clone(), res_vec[0].2, res_vec[0].3) + (out_energy / vmcparams.nthreads as f64, out_states_vec, out_de / vmcparams.nthreads as f64, out_corrtime / vmcparams.nthreads as f64) +} + +pub fn variationnal_monte_carlo( + rng: &mut [&mut R], + initial_state: &mut [FockState], params: &mut VarParams, - sys: &SysParams, - vmcparams: &VMCParams + sys: &mut SysParams, + vmcparams: &VMCParams, + param_map: &GenParameterMap, ) -> Vec -where T: BitOps + From + Display + Debug, Standard: Distribution +where T: BitOps + From + Display + Debug + Send + Sync, Standard: Distribution + std::marker::Send { let mut output_energy_array = vec![0.0; vmcparams.noptiter * 3]; - let mut otilde: Vec = vec![0.0; (4*sys.nfij + sys.nvij + sys.ngi) * (sys.nmcsample + 1)]; - //let mut work_otilde: Vec = vec![0.0; (sys.nfij + sys.nvij + sys.ngi) * (sys.nmcsample + 1)]; - let mut expvalo: Vec = vec![0.0; 4*sys.nfij + sys.nvij + sys.ngi]; - //let mut work_expvalo: Vec = vec![0.0; sys.nfij + sys.nvij + sys.ngi]; - let mut expval_ho: Vec = vec![0.0; 4*sys.nfij + sys.nvij + sys.ngi]; - //let mut work_expval_ho: Vec = vec![0.0; sys.nfij + sys.nvij + sys.ngi]; - let mut visited: Vec = vec![0; sys.nmcsample + 1]; - //let mut work_visited: Vec = vec![0; sys.nmcsample + 1]; - let mut der = DerivativeOperator { - o_tilde: &mut otilde, - expval_o: &mut expvalo, - ho: &mut expval_ho, - n: (4*sys.nfij + sys.nvij + sys.ngi) as i32, - nsamp: match vmcparams.compute_energy_method { + + let mut x0 = vec![0.0; sys.ngi + sys.nvij + sys.nfij].into_boxed_slice(); + let mut b = vec![0.0; sys.ngi + sys.nvij + sys.nfij].into_boxed_slice(); + let mut der = DerivativeOperator::new( + vmcparams.nparams as i32, + -1, + match vmcparams.compute_energy_method { EnergyComputationMethod::ExactSum => 1.0, - EnergyComputationMethod::MonteCarlo => sys.nmcsample as f64, + EnergyComputationMethod::MonteCarlo => (sys.nmcsample * vmcparams.nthreads) as f64, }, - nsamp_int: sys.mcsample_interval, - mu: -1, - visited: &mut visited, - pfaff_off: sys.ngi + sys.nvij, - jas_off: sys.ngi, - epsilon: vmcparams.epsilon, - }; + sys.mcsample_interval, + param_map.n_independant_jastrow + param_map.n_independant_gutzwiller, + param_map.n_independant_gutzwiller, + vmcparams.epsilon + ); + let mut work_der_vec = Vec::new(); + for i in 0..vmcparams.nthreads { + let mut work_der = param_map.mapto_general_representation(&der, &mut x0); + work_der_vec.push(work_der); + } + //let mut der = DerivativeOperator { + // o_tilde: &mut otilde, + // expval_o: &mut expvalo, + // ho: &mut expval_ho, + // n: (4*sys.nfij + sys.nvij + sys.ngi) as i32, + // nsamp: match vmcparams.compute_energy_method { + // EnergyComputationMethod::ExactSum => 1.0, + // EnergyComputationMethod::MonteCarlo => sys.nmcsample as f64, + // }, + // nsamp_int: sys.mcsample_interval, + // mu: -1, + // visited: &mut visited, + // pfaff_off: sys.ngi + sys.nvij, + // jas_off: sys.ngi, + // epsilon: vmcparams.epsilon, + //}; //let mut work_derivative = DerivativeOperator { // o_tilde: &mut work_otilde, // expval_o: &mut work_expvalo, @@ -103,14 +228,18 @@ where T: BitOps + From + Display + Debug, Standard: Distribution // epsilon: vmcparams.epsilon, //}; for opt_iter in 0..vmcparams.noptiter { + sys._opt_iter = opt_iter; let (mean_energy, _accumulated_states, deltae, correlation_time) = { match vmcparams.compute_energy_method { - EnergyComputationMethod::MonteCarlo => compute_mean_energy(rng, initial_state, params, sys, &mut der), + EnergyComputationMethod::MonteCarlo => { + parallel_monte_carlo(rng, initial_state, &mut work_der_vec, vmcparams, params, sys) + }, EnergyComputationMethod::ExactSum => { - (compute_mean_energy_exact(params, sys, &mut der), Vec::with_capacity(0), 0.0, 0.0) + (compute_mean_energy_exact(params, sys, &mut work_der_vec[0]), Vec::with_capacity(0), 0.0, 0.0) }, } }; + update_initial_state(initial_state, &_accumulated_states, sys.nmcsample, vmcparams.nthreads); // Save energy, error and correlation_time. output_energy_array[opt_iter * 3] = mean_energy; output_energy_array[opt_iter * 3 + 1] = deltae; @@ -125,18 +254,19 @@ where T: BitOps + From + Display + Debug, Standard: Distribution //} //mapto_pairwf(&derivative, &mut work_derivative, sys); - let mut x0 = vec![0.0; sys.nfij + sys.nvij + sys.ngi]; x0[(sys.ngi + sys.nvij)..(sys.ngi + sys.nvij + sys.nfij)].copy_from_slice(params.fij); x0[sys.ngi..(sys.ngi + sys.nvij)].copy_from_slice(params.vij); x0[0..sys.ngi].copy_from_slice(params.gi); + //let work_der = merge_derivatives(&work_der_vec, sys.nmcsample, vmcparams.nthreads); + param_map.update_reduced_representation(&work_der_vec[0], &mut der, &mut x0); // 68 misawa - let mut b: Vec = vec![0.0; der.n as usize]; + //let mut b: Vec = vec![0.0; der.n as usize]; unsafe { let incx = 1; let incy = 1; - daxpy(der.n, -mean_energy, der.expval_o, incx, der.ho, incy); - dcopy(der.n, der.ho, incx, &mut b, incy); + daxpy(der.n, -mean_energy, &der.expval_o, incx, &mut der.ho, incy); + dcopy(der.n, &der.ho, incx, &mut b, incy); } let mut _flag: bool = true; @@ -148,7 +278,7 @@ where T: BitOps + From + Display + Debug, Standard: Distribution conjugate_gradiant(&der, &mut b, &mut x0, vmcparams.epsilon, vmcparams.kmax, vmcparams.nparams as i32, vmcparams.threshold, vmcparams.epsilon_cg) }, }; - let mut delta_alpha = vec![0.0; vmcparams.nparams]; + let mut delta_alpha = vec![0.0; param_map.gendim as usize]; let mut j: usize = 0; for i in 0..vmcparams.nparams { if ignored_columns[i] { @@ -160,6 +290,8 @@ where T: BitOps + From + Display + Debug, Standard: Distribution _flag = false; } } + //panic!("Stop!"); + param_map.update_delta_alpha_reduced_to_gen(&mut delta_alpha); if vmcparams.optimise { unsafe { let incx = 1; @@ -208,6 +340,14 @@ where T: BitOps + From + Display + Debug, Standard: Distribution } // HARD CODE vij = vji // Slater Rescaling + let opt_delta = unsafe { + let incx = 1; + dnrm2(der.n, &delta_alpha, incx) + }; + if opt_delta <= vmcparams.conv_param_threshold { + println!("Exit early, achieved convergence within {} iteration, update now under supplied threshold.", opt_iter+1); + return output_energy_array; + } unsafe { let incx = 1; let max = params.fij[idamax(sys.nfij as i32, ¶ms.fij, incx) - 1]; @@ -228,6 +368,9 @@ where T: BitOps + From + Display + Debug, Standard: Distribution // ); //} zero_out_derivatives(&mut der, sys); + for i in 0..vmcparams.nthreads { + zero_out_derivatives(&mut work_der_vec[i], sys); + } //print_delta_alpha(&delta_alpha, sys.ngi, sys.nvij, sys.nfij); //let opt_delta = unsafe { // let incx = 1; diff --git a/src/fock_state.rs b/src/fock_state.rs index 3602c1d..c2df100 100644 --- a/src/fock_state.rs +++ b/src/fock_state.rs @@ -5,6 +5,7 @@ use rand::Rng; use rand::distributions::{Distribution, Standard}; use std::fmt; use log::error; +use std::marker::Send; use strings::{UPARROW, DOWNARROW}; @@ -40,7 +41,8 @@ pub trait BitOps: std::ops::Shr + std::ops::Shl + std::ops::BitOr + - From + From + + Sync { /// Provides the number of leading zeros in the bitstring. This gives the /// position of the first set bit in the string. This method is consistent @@ -72,7 +74,7 @@ pub trait BitOps: /// built-in methods. [BitOps::set] and [BitOps::check] are implemented by /// shifting a bitmask. impl BitOps for I - where I: PrimInt + std::ops::BitXorAssign + std::ops::BitAndAssign + From + From + where I: PrimInt + std::ops::BitXorAssign + std::ops::BitAndAssign + From + From + Sync { #[inline(always)] fn leading_zeros(self) -> u32 { @@ -127,6 +129,18 @@ impl From for u8 { } } +impl From for u16 { + fn from(other: SpinState) -> Self { + const NBYTES: usize = 2; + if ARRAY_SIZE > NBYTES {error!("Casting SpinState of size {} to discards overflowing bits.", SIZE)} + let mut tmp: [u8; NBYTES] = [0x00; NBYTES]; + for i in 0..ARRAY_SIZE { + tmp[i] = other.state[i]; + } + ::from_be_bytes(tmp) + } +} + impl From for u32 { fn from(other: SpinState) -> Self { const NBYTES: usize = 4; @@ -135,7 +149,7 @@ impl From for u32 { for i in 0..ARRAY_SIZE { tmp[i] = other.state[i]; } - ::from_ne_bytes(tmp) + ::from_be_bytes(tmp) } } @@ -147,7 +161,7 @@ impl From for u64 { for i in 0..ARRAY_SIZE { tmp[i] = other.state[i]; } - ::from_ne_bytes(tmp) + ::from_be_bytes(tmp) } } @@ -159,7 +173,7 @@ impl From for u128 { for i in 0..ARRAY_SIZE { tmp[i] = other.state[i]; } - ::from_ne_bytes(tmp) + ::from_be_bytes(tmp) } } @@ -375,6 +389,15 @@ impl std::ops::Shl for SpinState { } } +impl fmt::Display for SpinState { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + for i in 0..ARRAY_SIZE { + write!(f, "{:0b} ", self.state[i])?; + } + write!(f, "") + } +} + /// The Fock state structure. Encodes the spins positions. /// # Definition /// This structure has two different fields, the spin up and spin down component. @@ -392,14 +415,14 @@ impl std::ops::Shl for SpinState { /// assert_eq!(state_5_both << 2, state_20_both); /// ``` #[derive(Debug, Eq, PartialEq, Clone, Copy)] -pub struct FockState where T: From +pub struct FockState where T: From + std::fmt::Display + Send { pub spin_up: T, pub spin_down: T, pub n_sites: usize, } -impl + From> std::ops::Shl for FockState { +impl + From + std::fmt::Display + Send> std::ops::Shl for FockState { type Output = Self; fn shl(self, u: usize) -> Self::Output { @@ -411,7 +434,7 @@ impl + From> std::ops::Shl } } -impl + From> std::ops::Shr for FockState { +impl + From + std::fmt::Display + Send> std::ops::Shr for FockState { type Output = Self; fn shr(self, u: usize) -> Self::Output { @@ -423,7 +446,7 @@ impl + From> std::ops::Shr } } -impl fmt::Display for FockState { +impl fmt::Display for FockState { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { write!(f, "|")?; for i in 0..self.n_sites { @@ -449,7 +472,7 @@ pub trait Hopper { fn generate_all_exchange(self: &Self) -> Vec<(usize, usize)>; } -impl> Hopper for FockState { +impl + std::fmt::Display + Send> Hopper for FockState { fn generate_all_hoppings(self: &FockState, bitmask: &[SpinState]) -> Vec<(usize, usize, Spin)> { let sup = self.spin_up; let sdo = self.spin_down; @@ -460,6 +483,7 @@ impl> Hopper for FockState { for j in 1..self.n_sites /2 + 1 { // See note 2024-06-08, rotate right gives all the horizontal links. // We need to go to j up to SIZE / 2 + // let bitmask = ::from(bitmask[j - 1]); let mut possible_hoppings_up = (rot_right(sup, j, self.n_sites) ^ sup) & bitmask; let mut possible_hoppings_do = (rot_right(sdo, j, self.n_sites) ^ sdo) & bitmask; @@ -543,7 +567,7 @@ impl> Hopper for FockState { } // Interface for random state generation -impl Distribution> for Standard where Standard: Distribution { +impl Distribution> for Standard where Standard: Distribution { fn sample(&self, rng: &mut R) -> FockState { let su = rng.gen::(); let sd = rng.gen::(); @@ -562,7 +586,7 @@ pub trait RandomStateGeneration { } impl RandomStateGeneration for FockState -where T: BitOps, +where T: BitOps + std::fmt::Display + Send, Standard: Distribution { fn generate_from_nelec(rng: &mut R, nelec: usize, max_size: usize) -> FockState { diff --git a/src/gutzwiller.rs b/src/gutzwiller.rs index 85693fb..a79664b 100644 --- a/src/gutzwiller.rs +++ b/src/gutzwiller.rs @@ -34,7 +34,7 @@ pub fn compute_gutzwiller_exp( n_sites: usize, ) -> f64 where - T: BitOps + From + std::ops::Shl + Debug, + T: BitOps + From + std::ops::Shl + Debug + std::fmt::Display + Send, { // sum_i g_i n_i up n_i down let mut gutzwiller_sites = fock_state.spin_up & fock_state.spin_down; @@ -56,7 +56,7 @@ pub fn compute_gutzwiller_der( der: &mut DerivativeOperator, ) where - T: BitOps + From + std::ops::Shl + Debug, + T: BitOps + From + std::ops::Shl + Debug + std::fmt::Display + Send, { // sum_i g_i n_i up n_i down let mut gutzwiller_sites = fock_state.spin_up & fock_state.spin_down; diff --git a/src/hamiltonian.rs b/src/hamiltonian.rs index dff737b..38ca019 100644 --- a/src/hamiltonian.rs +++ b/src/hamiltonian.rs @@ -19,7 +19,7 @@ use crate::{VarParams, Spin, SysParams}; /// $$ pub fn potential(state: FockState, proj: f64, pstate: &PfaffianState, sys: &SysParams) -> f64 where - T: BitOps, + T: BitOps + std::fmt::Display + Send, { let pot = ((state.spin_up & state.spin_down).count_ones() as f64) * sys.cons_u; trace!("Output potential = {:.2} for state |x> = {}", pot, state); @@ -40,7 +40,7 @@ where /// $$ pub fn kinetic(state: FockState, previous_pstate: &PfaffianState, previous_proj: f64, params: &VarParams, sys: &SysParams) -> f64 where - T: BitOps + From + std::fmt::Debug + std::fmt::Display + T: BitOps + From + std::fmt::Debug + std::fmt::Display + Send { let hops = state.generate_all_hoppings(&sys.hopping_bitmask); diff --git a/src/hoppings.rs b/src/hoppings.rs index 656b27a..eee2c0d 100644 --- a/src/hoppings.rs +++ b/src/hoppings.rs @@ -1,3 +1,5 @@ +// Problem when compiling with N = ARRAY_SIZE, one array is out of bound, +// Should look into this, seems like a waste to use u32 for 16 Sites for example pub fn generate_bitmask(transfer_matrix: &[f64], size: usize) -> Vec { const WORD_SIZE: usize = 8; let all_zeros = SpinState{state: [0x00; ARRAY_SIZE]}; diff --git a/src/jastrow.rs b/src/jastrow.rs index 87d966e..ba90bcf 100644 --- a/src/jastrow.rs +++ b/src/jastrow.rs @@ -34,7 +34,7 @@ pub fn compute_jastrow_exp( n_sites: usize, ) -> f64 where - T: BitOps + std::fmt::Display, + T: BitOps + std::fmt::Display + Send, { let mut jastrow_out = 0.0; let mut regular_nor = !(fock_state.spin_up ^ fock_state.spin_down); @@ -73,7 +73,7 @@ pub fn compute_jastrow_der( n_sites: usize, ) where - T: BitOps + std::fmt::Display, + T: BitOps + std::fmt::Display + Send, { let mut regular_nor = !(fock_state.spin_up ^ fock_state.spin_down); regular_nor.mask_bits(n_sites); @@ -113,7 +113,7 @@ fn jastrow_undo_update( index_skip: usize, n_sites: usize, ) where - T: BitOps, + T: BitOps + std::fmt::Display + Send, { if spin_mask.check(index_j) { //*spin_mask = *spin_mask & (::ones() >> (index_j + 1)); @@ -152,7 +152,7 @@ fn jastrow_do_update( index_skip: usize, n_sites: usize, ) where - T: BitOps, + T: BitOps + std::fmt::Display + Send, { if spin_mask.check(index_j) { //*spin_mask = *spin_mask & (::ones() >> (index_j + 1)); @@ -192,7 +192,7 @@ fn jastrow_single_update( index_i: usize, sign: bool, ) where - T: BitOps, + T: BitOps + std::fmt::Display + Send, { if spin_mask.check(index_j) & spin_mask.check(index_i) { let (n1, n2) = (fock_state.spin_up, fock_state.spin_down); @@ -272,7 +272,7 @@ pub fn fast_update_jastrow( previous_j: usize, new_j: usize, ) where - T: BitOps + std::fmt::Display, + T: BitOps + std::fmt::Display + Send, { // Undo previous fock state // 1/2 sum {i\neq j} v_{ij}(n_i - 1)(n_j - 1) diff --git a/src/lib.rs b/src/lib.rs index a3620cf..faa762d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -51,23 +51,41 @@ impl<'a> std::fmt::Display for VarParams<'a> { /// The operator $O_k$ /// TODOC -pub struct DerivativeOperator<'a> { - pub o_tilde: &'a mut [f64], - pub expval_o: &'a mut [f64], - pub ho: &'a mut [f64], +pub struct DerivativeOperator { /// Number of variationnal parameters pub n: i32, /// Number of monte-carlo sampling pub mu: i32, pub nsamp: f64, pub nsamp_int: usize, - pub visited: &'a mut [usize], pub pfaff_off: usize, pub jas_off: usize, pub epsilon: f64, + pub o_tilde: Box<[f64]>, + pub expval_o: Box<[f64]>, + pub ho: Box<[f64]>, + pub visited: Box<[usize]>, } -impl<'a> std::fmt::Display for DerivativeOperator<'a> { +impl DerivativeOperator { + fn new(n: i32, mu: i32, nsamp: f64, nsamp_int: usize, pfaff_off: usize, jas_off: usize, epsilon: f64) -> Self { + DerivativeOperator { + o_tilde: vec![0.0; n as usize * nsamp as usize].into_boxed_slice(), + expval_o: vec![0.0; n as usize].into_boxed_slice(), + ho: vec![0.0; n as usize].into_boxed_slice(), + visited: vec![0; nsamp as usize].into_boxed_slice(), + mu, + n, + nsamp, + nsamp_int, + pfaff_off, + jas_off, + epsilon + } + } +} + +impl<'a> std::fmt::Display for DerivativeOperator { fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::result::Result<(), std::fmt::Error> { let width = 10; let mut expval = " = ".to_owned(); @@ -99,7 +117,7 @@ fn _save_otilde(der: &DerivativeOperator) { let mut c = vec![0.0; (der.n * der.n) as usize]; println!("dim = {}", der.n * der.n); unsafe { - dgemm(b"N"[0], b"T"[0], der.n, der.n, der.mu, 1.0, der.o_tilde, der.n, der.o_tilde, der.n, 0.0, &mut c, der.n); + dgemm(b"N"[0], b"T"[0], der.n, der.n, der.mu, 1.0, &der.o_tilde, der.n, &der.o_tilde, der.n, 0.0, &mut c, der.n); } let mut outstr = "".to_owned(); outstr.push_str(&format!(" = ")); diff --git a/src/monte_carlo.rs b/src/monte_carlo.rs index 4ffcc83..198724c 100644 --- a/src/monte_carlo.rs +++ b/src/monte_carlo.rs @@ -2,6 +2,7 @@ use blas::daxpy; use log::{error, info, trace, warn}; use rand::distributions::{Distribution, Standard}; use rand::Rng; +use std::sync::Arc; use crate::gutzwiller::compute_gutzwiller_der; use crate::jastrow::compute_jastrow_der; @@ -34,7 +35,7 @@ use crate::hamiltonian::{kinetic, potential}; #[inline(always)] fn propose_hopping > +T: BitOps + std::fmt::Display + std::fmt::Debug + From + Send> ( state: &FockState, pfaff_state: &PfaffianState, @@ -54,7 +55,7 @@ T: BitOps + std::fmt::Display + std::fmt::Debug + From> } #[inline(always)] -fn compute_hamiltonian(state: FockState, pstate: &PfaffianState, proj: f64, params: &VarParams, sys: &SysParams) -> f64 { +fn compute_hamiltonian(state: FockState, pstate: &PfaffianState, proj: f64, params: &VarParams, sys: &SysParams) -> f64 { let kin = kinetic(state, pstate, proj, params, sys); let e = kin + potential(state, proj, pstate, sys); trace!("Hamiltonian application = {} for state: |x> = {}", e, state); @@ -62,7 +63,7 @@ fn compute_hamiltonian(state: F } #[inline(always)] -fn compute_derivative_operator> +fn compute_derivative_operator + Send> (state: FockState, pstate: &PfaffianState, der: &mut DerivativeOperator, sys: &SysParams) { compute_gutzwiller_der(state, sys.size, der); @@ -71,7 +72,7 @@ fn compute_derivative_operator>( +fn make_update + Send>( n: &mut usize, proj: &mut f64, proj_copy: &mut f64, @@ -126,7 +127,7 @@ fn warmup( params: &VarParams, sys: &SysParams ) -where T: BitOps + From + std::fmt::Debug + std::fmt::Display, +where T: BitOps + From + std::fmt::Debug + std::fmt::Display + Send, R: Rng + ?Sized, Standard: Distribution { @@ -284,7 +285,7 @@ pub fn compute_mean_energy_exact(params: &VarParams, sys: &SysParams, der: &mut pub fn compute_mean_energy > +T: BitOps + std::fmt::Debug + std::fmt::Display + From + Send> (rng: &mut R, initial_state: FockState, params: &VarParams, sys: &SysParams, derivatives: &mut DerivativeOperator) -> (f64, Vec>, f64, f64) where Standard: Distribution { @@ -313,9 +314,11 @@ where Standard: Distribution let mut n_values = vec![0; error_estimation_level]; let mut energy_bootstraped = 0.0; - info!("Starting the warmup phase."); - warmup(rng, &mut state, &mut hop, &mut proj, &mut proj_copy_persistent, &mut ratio_prod, &mut - pstate, &mut n_accepted_updates, params, sys); + if sys.nwarmupchains > sys._opt_iter { + info!("Starting the warmup phase."); + warmup(rng, &mut state, &mut hop, &mut proj, &mut proj_copy_persistent, &mut ratio_prod, &mut + pstate, &mut n_accepted_updates, params, sys); + } info!("Starting the sampling phase."); // MC Sampling diff --git a/src/optimisation.rs b/src/optimisation.rs index 0520a80..00b6d49 100644 --- a/src/optimisation.rs +++ b/src/optimisation.rs @@ -5,6 +5,349 @@ use colored::Colorize; use crate::DerivativeOperator; +/// Encoding of a map to a reduced representation +/// # Map encoding +/// This struct encodes a Map to and from the general representation with +/// all the parameters to and from a reduced representation where some parameters +/// are multiples of other. +/// This struct is usefull to reduce the dimension of the optimisation and also +/// reduces the amount of linearly dependant columns, which helps with conditionning. +pub struct GenParameterMap { + /// Number of independant parameters. + pub dim: i32, + /// Dimension of general representation + pub gendim: i32, + /// Number of general parameters + pub n_genparams: i32, + /// Number of independant gutzwiller parameters. + pub n_independant_gutzwiller: usize, + /// Number of independant jastrow parameters. + pub n_independant_jastrow: usize, + /// The projector from the general representation to the reduced one. + /// Should behave: + /// Ax => \tilde{x} where every exquivalent component of x are sumed in + /// the lowest index and zero elsewhere. + /// Example: + /// A = ((1, 1), (0, 0)) + /// x = (a, b) + /// \tilde{x} = (a + b, 0) + /// A will be triangular, so only the upper triangle needs to be stored. + /// If you transpose this matrix, it will reverse this operation somewhat + /// A^T = ((1, 0), (1, 0) + /// A\tilde{x} = (a + b, a + b) + /// This indeed copies the parameter x_0 to x_1. If x_1 is set to half x_0, + /// then + /// A = ((1, 2), (0, 0)) + /// x = (a, b), Ax = (a + 2b) + /// Then the coefficient wise-inverse of the transpose is usefull + /// A^T^{-1} = ((1, 0), (1/2, 0)) + /// x = (a, 0) + /// \tilde{x} = (a, 1/2 a) + /// It is then usefull to store only the lower triangle of the "inverse" of + /// the transpose, as coefficient don't matter when projecting, but do + /// matter when recovering. + /// + /// Indexing + /// n = dim (dim - 1) / 2 + i - (dim - j) * ( dim - j - 1) / 2 + /// Can only have one non zero value per row + /// + /// PERFORMANCE: + /// If performance is an issue, an idea would be to write this matrix inside + /// an integer. This is not done here as it removes to possibility to make + /// a parameter say half another one. + pub projector: Box<[f64]>, +} + +pub trait ReducibleGeneralRepresentation { + /// Overrides b and x0 + /// Safety: + /// b and x0 must have lenght self.gendim + fn mapto_general_representation(&self, der: &DerivativeOperator, x0: &mut [f64],) -> DerivativeOperator; + /// Overrides b and x0 + /// Safety: + /// b and x0 must have lenght self.gendim + /// On exit, the range [self.dim..self.gendim] contains garbage. + fn mapto_reduced_representation(&self, der: &DerivativeOperator, x0: &mut [f64],) -> DerivativeOperator; + /// Overrides b and x0 + /// Safety: + /// b and x0 must have lenght self.gendim + fn update_general_representation(&self, der: &DerivativeOperator, out_der: &mut DerivativeOperator, x0: &mut [f64],); + /// Overrides b and x0 + /// Safety: + /// b and x0 must have lenght self.gendim + /// On exit, the range [self.dim..self.gendim] contains garbage. + fn update_reduced_representation(&self, der: &DerivativeOperator, out_der: &mut DerivativeOperator, x0: &mut [f64],); + fn update_delta_alpha_reduced_to_gen(&self, x0: &mut [f64]); +} + +impl<'a> ReducibleGeneralRepresentation for GenParameterMap { + fn mapto_general_representation(&self, der: &DerivativeOperator, x0: &mut [f64]) -> DerivativeOperator + { + if self.dim == 0 { + error!("Tried to map from an empty parameter representation, this will + break conjugate gradient or diagonalisation. Dim = {}", self.dim); + panic!("Undefined behavior."); + } + + let mut out_x = vec![0.0; self.gendim as usize]; + let n = self.gendim as usize; + let mut out_der = DerivativeOperator::new(self.gendim, der.mu, der.nsamp, der.nsamp_int, + self.n_independant_gutzwiller + self.n_independant_jastrow, + self.n_independant_gutzwiller, der.epsilon); + if der.mu < 0 { + return out_der; + } + + // Copy the visited vector + for i in 0..der.mu as usize { + out_der.visited[i] = der.visited[i]; + } + + // Copy all other data + let mut k = 0; + let mut flag = false; + for j in 0..n { + for i in n..n { + let pij = self.projector[ + j + (i * (i+1) / 2) + ]; + if pij != 0.0 { + out_der.expval_o[i] = pij * der.expval_o[k]; + out_der.ho[i] = pij * der.ho[k]; + out_x[i] = pij * x0[k]; + unsafe { + let incx = self.dim; + let incy = self.gendim; + dcopy( + der.mu, + &der.o_tilde[i..(self.dim*der.mu) as usize], + incx, + &mut out_der.o_tilde[k..(self.gendim * der.mu) as usize], + incy + ); + } + flag = true; + //break; + } + } + if flag{ + k += 1; + flag = false; + } + } + unsafe { + let incx = 1; + let incy = 1; + dcopy(self.gendim, &out_x, incx, x0, incy); + } + + // Should return b and x0 as well + out_der + } + + fn mapto_reduced_representation(&self, der: &DerivativeOperator, x0: &mut [f64]) -> DerivativeOperator + { + if self.dim == 0 { + error!("Tried to map from an empty parameter representation, this will + break conjugate gradient or diagonalisation. Dim = {}", self.dim); + panic!("Undefined behavior."); + } + let n = self.gendim as usize; + let mut out_der = DerivativeOperator::new(self.dim, der.mu, der.nsamp, der.nsamp_int, + self.n_independant_gutzwiller + self.n_independant_jastrow, + self.n_independant_gutzwiller, der.epsilon); + if der.mu < 0 { + return out_der; + } + + // Copy the visited vector + for i in 0..der.mu as usize { + out_der.visited[i] = der.visited[i]; + } + + // Copy all other data + // iter over transpose's diagonal + let mut i: usize = 0; + for j in 0..n { + let pjj = self.projector[ + j + (j * (j+1) / 2) + ]; + if pjj != 0.0 { // If 0, parameter does not map to itself. + if pjj != 1.0 { + error!("Parameters must map to themselfs with multiple 1."); + panic!("Undefined behavior."); + } + out_der.expval_o[i] = der.expval_o[j]; + out_der.ho[i] = der.ho[j]; + // Modify the vectors inplace, we won't iterate on past values + x0[i] = x0[j]; + unsafe { + let incx = self.gendim; + let incy = self.dim; + dcopy( + der.mu, + &der.o_tilde[j..n*der.mu as usize], + incx, + &mut out_der.o_tilde[i..(self.dim * der.mu) as usize], + incy + ); + } + // Filling the independant parameter index + i += 1; + } + } + if i != self.dim as usize { + error!("Parameter map did not yield a dimension of {}, instead got {}", self.dim, i); + panic!("Undefined behavior."); + } + + out_der + } + + fn update_general_representation(&self, der: &DerivativeOperator, out_der: &mut DerivativeOperator, x0: &mut [f64]) { + if self.dim == 0 { + error!("Tried to map from an empty parameter representation, this will + break conjugate gradient or diagonalisation. Dim = {}", self.dim); + panic!("Undefined behavior."); + } + + let mut out_x = vec![0.0; self.gendim as usize]; + let n = self.gendim as usize; + out_der.mu = der.mu; + + // Copy the visited vector + for i in 0..der.mu as usize { + out_der.visited[i] = der.visited[i]; + } + + // Copy all other data + let mut k = 0; + let mut flag = false; + for j in 0..n { + for i in j..n { + let pij = self.projector[ + j + (i * (i+1) / 2) + ]; + if pij != 0.0 { + out_der.expval_o[i] = pij * der.expval_o[k]; + out_der.ho[i] = pij * der.ho[k]; + out_x[i] = pij * x0[k]; + unsafe { + let incx = self.dim; + let incy = self.gendim; + dcopy( + der.mu, + &der.o_tilde[i..(self.dim*der.mu) as usize], + incx, + &mut out_der.o_tilde[k..(self.gendim * der.mu) as usize], + incy + ); + } + flag = true; + //break; + } + } + if flag{ + k += 1; + flag = false; + } + } + unsafe { + let incx = 1; + let incy = 1; + dcopy(self.gendim, &out_x, incx, x0, incy); + } + } + + fn update_delta_alpha_reduced_to_gen(&self, x0: &mut [f64]) { + if self.dim == 0 { + error!("Tried to map from an empty parameter representation, this will + break conjugate gradient or diagonalisation. Dim = {}", self.dim); + panic!("Undefined behavior."); + } + let n = self.gendim as usize; + let mut out_x0 = vec![0.0; n]; + + // Copy all other data + // iter over transpose's diagonal + // Copy all other data + let mut k = 0; + let mut flag = false; + for j in 0..n { + for i in j..n { + let pij = self.projector[ + j + (i * (i+1) / 2) + ]; + if pij != 0.0 { + out_x0[i] = pij * x0[k]; + flag = true; + //break; + } + } + if flag{ + k += 1; + flag = false; + } + } + unsafe { + let incx = 1; + let incy = 1; + dcopy(self.gendim, &out_x0, incx, x0, incy); + } + } + + fn update_reduced_representation(&self, der: &DerivativeOperator, out_der: &mut DerivativeOperator, x0: &mut [f64]) { + if self.dim == 0 { + error!("Tried to map from an empty parameter representation, this will + break conjugate gradient or diagonalisation. Dim = {}", self.dim); + panic!("Undefined behavior."); + } + let n = self.gendim as usize; + out_der.mu = der.mu; + + // Copy the visited vector + for i in 0..der.mu as usize { + out_der.visited[i] = der.visited[i]; + } + + // Copy all other data + // iter over transpose's diagonal + let mut i = 0; + for j in 0..n { + let pjj = self.projector[ + j + (j * (j+1) / 2) + ]; + if pjj != 0.0 { // If 0, parameter does not map to itself. + if pjj != 1.0 { + error!("Parameters must map to themselfs with multiple 1."); + panic!("Undefined behavior."); + } + out_der.expval_o[i] = der.expval_o[j]; + out_der.ho[i] = der.ho[j]; + // Modify the vectors inplace, we won't iterate on past values + x0[i] = x0[j]; + unsafe { + let incx = self.gendim; + let incy = self.dim; + dcopy( + der.mu, + &der.o_tilde[j..n*der.mu as usize], + incx, + &mut out_der.o_tilde[i..(self.dim * der.mu) as usize], + incy + ); + } + // Filling the independant parameter index + i += 1; + } + } + if i != self.dim as usize { + error!("Parameter map did not yield a dimension of {}, instead got {}", self.dim, i); + panic!("Undefined behavior."); + } + } +} + fn gradient(x: &[f64], otilde: &[f64], visited: &[usize], expval_o: &[f64], b: &mut [f64], epsilon: f64, dim: i32, mu: i32, nsamp: f64) { let alpha = -1.0; let incx = 1; @@ -308,7 +651,7 @@ pub fn exact_overlap_inverse(a: &DerivativeOperator, b: &mut [f64], epsilon: f64 // PRE FILTER let mut ignore = vec![false; dim as usize]; //println!("{}", save_otilde(a.o_tilde, a.mu as usize, a.n as usize)); - let mut unfiltered_s = compute_s_explicit(a.o_tilde, a.expval_o, a.visited, dim, a.mu, a.nsamp, epsilon); + let mut unfiltered_s = compute_s_explicit(&a.o_tilde, &a.expval_o, &a.visited, dim, a.mu, a.nsamp, epsilon); //println!("dim = {}, Unfiltered S = ", dim); //println!("{}", save_otilde(&unfiltered_s, dim as usize, dim as usize)); let new_dim = prefilter_overlap_matrix(a, &mut ignore, dim, thresh); @@ -412,7 +755,7 @@ pub fn conjugate_gradiant(a: &DerivativeOperator, b: &mut [f64], x0: &mut [f64], e *= epsilon_convergence; trace!("Error threshold e = {}", e); //println!("Error threshold e = {}", e); - gradient(x0, õ, a.visited, &expvalo, b, epsilon, new_dim as i32, a.mu, a.nsamp); + gradient(x0, õ, &a.visited, &expvalo, b, epsilon, new_dim as i32, a.mu, a.nsamp); let mut p = vec![0.0; new_dim]; let mut j: usize = 0; for i in 0..dim as usize { @@ -432,7 +775,7 @@ pub fn conjugate_gradiant(a: &DerivativeOperator, b: &mut [f64], x0: &mut [f64], trace!("p_{} : {:?}", k, p); //println!("r_{} : {:?}", k, b); //println!("p_{} : {:?}", k, p); - compute_w(&mut w, õ, a.visited, &expvalo, &p, epsilon, new_dim as i32, a.mu, a.nsamp); + compute_w(&mut w, õ, &a.visited, &expvalo, &p, epsilon, new_dim as i32, a.mu, a.nsamp); //println!("w_{} : {:?}", k, w); let nrm2rk = alpha_k(b, &p, &w, &mut alpha, new_dim as i32); trace!("alpha_{} : {}", k, alpha); diff --git a/src/parse/mod.rs b/src/parse/mod.rs deleted file mode 100644 index 229ff7e..0000000 --- a/src/parse/mod.rs +++ /dev/null @@ -1,4 +0,0 @@ -pub use orbitale::parse_orbitale_def as orbitale_def; - -/// For orbitale.csv -pub mod orbitale; diff --git a/src/parse/orbitale.rs b/src/parse/orbitale.rs deleted file mode 100644 index 3185813..0000000 --- a/src/parse/orbitale.rs +++ /dev/null @@ -1,116 +0,0 @@ -use csv::{Reader, StringRecord}; -use derive_more::{Constructor, Error}; -use rand::Rng; -use std::fmt; -use std::fs::File; -use std::path::PathBuf; - -/// Parse the orbital variationnal definition -/// # Arguments -/// * __`fp`__ - File path to the definition file, in csv format. -/// * __`s`__ - Size of the system. -pub fn parse_orbitale_def(fp: &PathBuf, s: usize) -> Result> { - // Parameter reference vector. column packed. - let mut fij: Vec = vec![0.0; 4 * s * s]; - // Parameter vector. Ordered by "n". - let mut orbitale_params: Vec = Vec::new(); - let mut n_orbitale_params: usize = 0; - // RNG thread - let mut rng = rand::thread_rng(); - - // Parse input file. - let file = File::open(fp)?; - let mut reader = Reader::from_reader(file); - for (k, result) in reader.records().enumerate() { - let rec = result?; - // Should have 3 column. - if rec.len() != 3 { - let d = format!( - "Error at line {} in orbitale.csv, invalid number of elements.", - k - ); - let mut e = OrbitaleParseError::new("Invalid number of argument on a line.".to_owned()); - e.details = d; - return Err(e); - } - - // Parse element coordinates. - let i = parse_single_elem(&rec, 0, k)?; - let j = parse_single_elem(&rec, 1, k)?; - - let n = rec.get(2).unwrap().parse::(); - match n { - Ok(param) => { - if n_orbitale_params < param { - n_orbitale_params += 1; - orbitale_params.push(rng.gen::() * 2.0 - 1.0); - } - let var_param: f64 = orbitale_params[param - 1]; - fij[i + j * s] = var_param; - fij[j + i * s] = -var_param; - fij[i + j * s + s*s] = var_param; - fij[j + i * s + s*s] = -var_param; - fij[i + j * s + 2*s*s] = var_param; - fij[j + i * s + 2*s*s] = -var_param; - fij[i + j * s + 3*s*s] = var_param; - fij[j + i * s + 3*s*s] = -var_param; - } - Err(error) => { - // Add error message. - let d = format!("Error at line {}, invalid parameter identifier.", k); - let mut e = OrbitaleParseError::from(error); - e.details = d; - return Err(e); - } - }; - } - Ok(fij) -} - -fn parse_single_elem(line: &StringRecord, col: usize, l: usize) -> Result { - match line.get(col).unwrap().parse::() { - Ok(v) => Ok(v), - Err(error) => { - let details = format!( - "Expected valid coordinates in orbitale.csv at line {}, col {}", - l, col - ); - let mut e = OrbitaleParseError::from(error); - e.details = details; - Err(e) - } - } -} - -type Result = std::result::Result; - -/// Error in the orbital params definition. -#[derive(Debug, Clone, Error, Constructor)] -pub struct OrbitaleParseError { - pub details: String, -} - -impl fmt::Display for OrbitaleParseError { - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - write!(f, "Parsing error encountered in orbitale.csv.") - } -} - -impl From for OrbitaleParseError { - fn from(err: std::num::ParseIntError) -> Self { - println!("{}", err); - OrbitaleParseError::new("Expected to parse an integer.".to_owned()) - } -} - -impl From for OrbitaleParseError { - fn from(err: csv::Error) -> Self { - OrbitaleParseError::new(err.to_string()) - } -} - -impl From for OrbitaleParseError { - fn from(err: std::io::Error) -> Self { - OrbitaleParseError::new(err.to_string()) - } -} diff --git a/src/pfaffian.rs b/src/pfaffian.rs index ed769a6..684332b 100644 --- a/src/pfaffian.rs +++ b/src/pfaffian.rs @@ -127,7 +127,7 @@ fn transpose(a: &Vec, n: usize) -> Vec{ /// * __`state`__ - The state of the system. pub fn construct_matrix_a_from_state(fij: &[f64], state: FockState, sys: &SysParams) -> PfaffianState where - T: BitOps + std::fmt::Display, + T: BitOps + std::fmt::Display + Send, { // Fij upup, updown, downup, downdown let n = state.spin_up.count_ones() as usize + state.spin_down.count_ones() as usize; @@ -658,9 +658,11 @@ mod tests { clean_update_frequency: 0, nmcwarmup: 0, nmcsample: 0, + nwarmupchains: 0, tolerance_sherman_morrison: 0.0, tolerance_singularity: 0.0, pair_wavefunction: false, + _opt_iter: 0, }; let mut rng = SmallRng::seed_from_u64(42); // Size of the system @@ -797,9 +799,11 @@ mod tests { clean_update_frequency: 0, nmcwarmup: 0, nmcsample: 0, + nwarmupchains: 0, tolerance_sherman_morrison: 0.0, tolerance_singularity: 0.0, pair_wavefunction: false, + _opt_iter: 0, }; let mut params = vec![0.0; 4 * SIZE * SIZE]; // params[i+8*j] = f_ij @@ -840,9 +844,11 @@ mod tests { clean_update_frequency: 0, nmcwarmup: 0, nmcsample: 0, + nwarmupchains: 0, tolerance_sherman_morrison: 0.0, tolerance_singularity: 0.0, pair_wavefunction: false, + _opt_iter: 0, }; let mut params = vec![0.0; 4 * SIZE * SIZE]; // params[i+8*j] = f_ij @@ -902,9 +908,11 @@ mod tests { clean_update_frequency: 0, nmcwarmup: 0, nmcsample: 0, + nwarmupchains: 0, tolerance_sherman_morrison: 0.0, tolerance_singularity: 0.0, pair_wavefunction: false, + _opt_iter: 0, }; let mut params = vec![0.0; 4 * SIZE * SIZE]; // params[i+8*j] = f_ij @@ -982,9 +990,11 @@ mod tests { clean_update_frequency: 0, nmcwarmup: 0, nmcsample: 0, + nwarmupchains: 0, tolerance_sherman_morrison: 0.0, tolerance_singularity: 0.0, pair_wavefunction: false, + _opt_iter: 0, }; let mut params = vec![0.0; 4 * SIZE * SIZE]; // params[i+8*j] = f_ij diff --git a/tests/orbitales_def/test_csv_parse_error.csv b/tests/orbitales_def/test_csv_parse_error.csv deleted file mode 100644 index 96e3e80..0000000 --- a/tests/orbitales_def/test_csv_parse_error.csv +++ /dev/null @@ -1,2 +0,0 @@ -i,j,n -0,1,1,1 diff --git a/tests/orbitales_def/test_int_parse_error.csv b/tests/orbitales_def/test_int_parse_error.csv deleted file mode 100644 index f24926c..0000000 --- a/tests/orbitales_def/test_int_parse_error.csv +++ /dev/null @@ -1,2 +0,0 @@ -i,j,n -0,1,k diff --git a/tests/orbitales_def/test_parse_2x2.csv b/tests/orbitales_def/test_parse_2x2.csv deleted file mode 100644 index eb2b5e3..0000000 --- a/tests/orbitales_def/test_parse_2x2.csv +++ /dev/null @@ -1,2 +0,0 @@ -i,j,n -0,1,1 diff --git a/tests/orbitales_def/test_parse_3x3.csv b/tests/orbitales_def/test_parse_3x3.csv deleted file mode 100644 index ca3e18a..0000000 --- a/tests/orbitales_def/test_parse_3x3.csv +++ /dev/null @@ -1,4 +0,0 @@ -i,j,n -0,1,1 -0,2,2 -1,2,3 diff --git a/tests/parse_orbitale_test.rs b/tests/parse_orbitale_test.rs deleted file mode 100644 index 9f77418..0000000 --- a/tests/parse_orbitale_test.rs +++ /dev/null @@ -1,51 +0,0 @@ -use std::path::Path; - -use impurity::parse::orbitale::parse_orbitale_def; - -#[test] -#[ignore] -fn test_parse_2x2() { - let size = 2; - let fp = Path::new("tests/orbitales_def/test_parse_2x2.csv"); - let fij = parse_orbitale_def(&fp.to_path_buf(), size).unwrap(); - assert_eq!(fij.len(), 4); - assert_eq!(fij[1], -fij[2]); - assert!(fij[1] <= 0.0); -} - -#[test] -#[ignore] -fn test_parse_3x3() { - let size = 3; - let fp = Path::new("tests/orbitales_def/test_parse_3x3.csv"); - let fij = parse_orbitale_def(&fp.to_path_buf(), size).unwrap(); - assert_eq!(fij.len(), 9); - assert_eq!(fij[1], -fij[3]); - assert_eq!(fij[6], -fij[2]); - assert_eq!(fij[5], -fij[7]); - assert!(fij[1] <= 0.0); - assert!(fij[2] <= 0.0); - assert!(fij[5] <= 0.0); -} - -#[test] -fn test_int_parse_error() { - let size = 2; - let fp = Path::new("tests/orbitales_def/test_int_parse_error.csv"); - let fij = parse_orbitale_def(&fp.to_path_buf(), size); - match fij { - Ok(_) => panic!("Should have errored."), - Err(_) => println!("Error as expected."), - } -} - -#[test] -fn test_csv_parse_error() { - let size = 2; - let fp = Path::new("tests/orbitales_def/test_csv_parse_error.csv"); - let fij = parse_orbitale_def(&fp.to_path_buf(), size); - match fij { - Ok(_) => panic!("Should have errored."), - Err(_) => println!("Error as expected."), - } -} diff --git a/tests/vmc_2sites.rs b/tests/vmc_2sites.rs index 66e6be8..0eb4024 100644 --- a/tests/vmc_2sites.rs +++ b/tests/vmc_2sites.rs @@ -357,9 +357,11 @@ fn comupte_energy_from_all_states() { clean_update_frequency: CLEAN_UPDATE_FREQUENCY, nmcwarmup: NMCWARMUP, nmcsample: NMCSAMP, + nwarmupchains: 1, tolerance_sherman_morrison: TOL_SHERMAN, tolerance_singularity: TOL_SINGULARITY, pair_wavefunction: false, + _opt_iter: 0, }; let mut fij = [ 0.0, 0.0, 0.0, 0.0, @@ -423,19 +425,19 @@ fn comupte_energy_from_all_states() { energy_individual_state(&states_names[i], ¶meters), 1e-12); } close(mean_energy, analytic(¶meters), 1e-12); - let mut otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; - let mut expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut visited: Vec = vec![0; NMCSAMP]; + let otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited: Vec = vec![0; NMCSAMP]; let mut der = DerivativeOperator { - o_tilde: &mut otilde, - expval_o: &mut expvalo, - ho: &mut expval_ho, + o_tilde: otilde.into_boxed_slice(), + expval_o: expvalo.into_boxed_slice(), + ho: expval_ho.into_boxed_slice(), n: (NFIJ + NVIJ + NGI) as i32, nsamp: NMCSAMP as f64, nsamp_int: 1, mu: -1, - visited: &mut visited, + visited: visited.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: 0.0, @@ -466,7 +468,7 @@ fn comupte_energy_from_all_states() { // Test derivatives let exp_val = analytic_derivatives_expval(¶meters); - print_der(der.expval_o, &exp_val, sys.ngi+sys.nvij+sys.nfij); + print_der(&der.expval_o, &exp_val, sys.ngi+sys.nvij+sys.nfij); let psi = norm(¶meters); println!("Norm: {:10.4e}", psi); for i in 0..sys.ngi+sys.nvij+sys.nfij { @@ -474,7 +476,7 @@ fn comupte_energy_from_all_states() { } let exp_val_ho = analytic_ho_expval(¶meters); - print_der(der.ho, &exp_val_ho, sys.ngi+sys.nvij+sys.nfij); + print_der(&der.ho, &exp_val_ho, sys.ngi+sys.nvij+sys.nfij); let psi = norm(¶meters); println!("Norm: {:10.4e}", psi); for i in 0..sys.ngi+sys.nvij+sys.nfij { diff --git a/tests/vmc_2sites_analytic_conv.rs b/tests/vmc_2sites_analytic_conv.rs index 7ea254f..460f10d 100644 --- a/tests/vmc_2sites_analytic_conv.rs +++ b/tests/vmc_2sites_analytic_conv.rs @@ -326,9 +326,11 @@ fn comupte_energy_from_all_states() { clean_update_frequency: CLEAN_UPDATE_FREQUENCY, nmcwarmup: NMCWARMUP, nmcsample: NMCSAMP, + nwarmupchains: 1, tolerance_sherman_morrison: TOL_SHERMAN, tolerance_singularity: TOL_SINGULARITY, pair_wavefunction: false, + _opt_iter: 0, }; let mut fij = [ 0.0, 0.0, 0.0, 0.0, @@ -393,36 +395,36 @@ fn comupte_energy_from_all_states() { } close(mean_energy, analytic(¶meters), 1e-12); mean_energy = mean_energy / norm(¶meters); - let mut otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; - let mut expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut visited: Vec = vec![0; NMCSAMP]; + let otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited: Vec = vec![0; NMCSAMP]; let mut der = DerivativeOperator { - o_tilde: &mut otilde, - expval_o: &mut expvalo, - ho: &mut expval_ho, + o_tilde: otilde.into_boxed_slice(), + expval_o: expvalo.into_boxed_slice(), + ho: expval_ho.into_boxed_slice(), n: (NFIJ + NVIJ + NGI) as i32, nsamp: NMCSAMP as f64, nsamp_int: 1, mu: -1, - visited: &mut visited, + visited: visited.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: 0.0, }; - let mut otilde_pair: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; - let mut expvalo_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut expval_ho_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut visited_pair: Vec = vec![0; NMCSAMP]; + let otilde_pair: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let expvalo_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited_pair: Vec = vec![0; NMCSAMP]; let mut der_pair = DerivativeOperator { - o_tilde: &mut otilde_pair, - expval_o: &mut expvalo_pair, - ho: &mut expval_ho_pair, + o_tilde: otilde_pair.into_boxed_slice(), + expval_o: expvalo_pair.into_boxed_slice(), + ho: expval_ho_pair.into_boxed_slice(), n: (sys.size*sys.size + NVIJ + NGI) as i32, nsamp: NMCSAMP as f64, nsamp_int: 1, mu: -1, - visited: &mut visited_pair, + visited: visited_pair.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: 0.0, @@ -448,7 +450,7 @@ fn comupte_energy_from_all_states() { // Test derivatives let exp_val = analytic_derivatives_expval(¶meters); println!("Checking "); - print_der(der_pair.expval_o, &exp_val, sys.ngi+sys.nvij+sys.nfij); + print_der(&der_pair.expval_o, &exp_val, sys.ngi+sys.nvij+sys.nfij); let psi = norm(¶meters); println!("Norm: {:10.4e}", psi); for i in 0..sys.ngi+sys.nvij+sys.nfij { @@ -458,7 +460,7 @@ fn comupte_energy_from_all_states() { let exp_val_ho = analytic_ho_expval(¶meters); println!("Checking "); - print_der(der_pair.ho, &exp_val_ho, sys.ngi+sys.nvij+sys.nfij); + print_der(&der_pair.ho, &exp_val_ho, sys.ngi+sys.nvij+sys.nfij); let psi = norm(¶meters); println!("Norm: {:10.4e}", psi); for i in 0..sys.ngi+sys.nvij+sys.nfij { @@ -478,8 +480,8 @@ fn comupte_energy_from_all_states() { unsafe { let incx = 1; let incy = 1; - daxpy(der_pair.n, -mean_energy, der_pair.expval_o, incx, der_pair.ho, incy); - dcopy(der_pair.n, der_pair.ho, incx, &mut b, incy); + daxpy(der_pair.n, -mean_energy, &der_pair.expval_o, incx, &mut der_pair.ho, incy); + dcopy(der_pair.n, &der_pair.ho, incx, &mut b, incy); } println!("x0 = {:?}", x0); } diff --git a/tests/vmc_2sites_exact_sum.rs b/tests/vmc_2sites_exact_sum.rs index 048d67a..29eb7b6 100644 --- a/tests/vmc_2sites_exact_sum.rs +++ b/tests/vmc_2sites_exact_sum.rs @@ -322,9 +322,11 @@ fn comupte_energy_from_all_states() { clean_update_frequency: CLEAN_UPDATE_FREQUENCY, nmcwarmup: NMCWARMUP, nmcsample: NMCSAMP, + nwarmupchains: 1, tolerance_sherman_morrison: TOL_SHERMAN, tolerance_singularity: TOL_SINGULARITY, pair_wavefunction: false, + _opt_iter: 0, }; let mut fij = [ 0.0, 0.0, 0.0, 0.0, @@ -389,36 +391,36 @@ fn comupte_energy_from_all_states() { } close(mean_energy, analytic(¶meters), 1e-12); mean_energy = mean_energy / norm(¶meters); - let mut otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; - let mut expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut visited: Vec = vec![0; NMCSAMP]; + let otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited: Vec = vec![0; NMCSAMP]; let mut der = DerivativeOperator { - o_tilde: &mut otilde, - expval_o: &mut expvalo, - ho: &mut expval_ho, + o_tilde: otilde.into_boxed_slice(), + expval_o: expvalo.into_boxed_slice(), + ho: expval_ho.into_boxed_slice(), n: (NFIJ + NVIJ + NGI) as i32, nsamp: NMCSAMP as f64, nsamp_int: 1, mu: -1, - visited: &mut visited, + visited: visited.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: 0.0, }; - let mut otilde_pair: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; - let mut expvalo_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut expval_ho_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut visited_pair: Vec = vec![0; NMCSAMP]; + let otilde_pair: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let expvalo_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho_pair: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited_pair: Vec = vec![0; NMCSAMP]; let mut der_pair = DerivativeOperator { - o_tilde: &mut otilde_pair, - expval_o: &mut expvalo_pair, - ho: &mut expval_ho_pair, + o_tilde: otilde_pair.into_boxed_slice(), + expval_o: expvalo_pair.into_boxed_slice(), + ho: expval_ho_pair.into_boxed_slice(), n: (sys.size*sys.size + NVIJ + NGI) as i32, nsamp: NMCSAMP as f64, nsamp_int: 1, mu: -1, - visited: &mut visited_pair, + visited: visited_pair.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: 0.0, @@ -444,7 +446,7 @@ fn comupte_energy_from_all_states() { // Test derivatives let exp_val = analytic_derivatives_expval(¶meters); println!("Checking "); - print_der(der_pair.expval_o, &exp_val, sys.ngi+sys.nvij+sys.nfij); + print_der(&der_pair.expval_o, &exp_val, sys.ngi+sys.nvij+sys.nfij); let psi = norm(¶meters); println!("Norm: {:10.4e}", psi); for i in 0..sys.ngi+sys.nvij+sys.nfij { @@ -454,7 +456,7 @@ fn comupte_energy_from_all_states() { let exp_val_ho = analytic_ho_expval(¶meters); println!("Checking "); - print_der(der_pair.ho, &exp_val_ho, sys.ngi+sys.nvij+sys.nfij); + print_der(&der_pair.ho, &exp_val_ho, sys.ngi+sys.nvij+sys.nfij); let psi = norm(¶meters); println!("Norm: {:10.4e}", psi); for i in 0..sys.ngi+sys.nvij+sys.nfij { diff --git a/tests/vmc_clean_updates.rs b/tests/vmc_clean_updates.rs index 855c7b2..4e67ab4 100644 --- a/tests/vmc_clean_updates.rs +++ b/tests/vmc_clean_updates.rs @@ -86,25 +86,27 @@ fn monte_carlo_first_iteration() { clean_update_frequency: CLEAN_UPDATE_FREQUENCY, nmcwarmup: NMCWARMUP, nmcsample: NMCSAMP, + nwarmupchains: 1, tolerance_sherman_morrison: TOLERENCE_SHERMAN_MORRISSON, tolerance_singularity: TOLERANCE_SINGULARITY, pair_wavefunction: false, + _opt_iter: 0, }; log_system_parameters(&sys); - let mut otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; - let mut expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; - let mut visited: Vec = vec![0; NMCSAMP]; + let otilde: Vec = vec![0.0; (NFIJ + NVIJ + NGI) * NMCSAMP]; + let expvalo: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let expval_ho: Vec = vec![0.0; NFIJ + NVIJ + NGI]; + let visited: Vec = vec![0; NMCSAMP]; let mut der = DerivativeOperator { - o_tilde: &mut otilde, - expval_o: &mut expvalo, - ho: &mut expval_ho, + o_tilde: otilde.into_boxed_slice(), + expval_o: expvalo.into_boxed_slice(), + ho: expval_ho.into_boxed_slice(), n: (NFIJ + NVIJ + NGI) as i32, nsamp: NMCSAMP as f64, nsamp_int: 1, mu: -1, - visited: &mut visited, + visited: visited.into_boxed_slice(), pfaff_off: NGI + NVIJ, jas_off: NGI, epsilon: 0.0,