From 00518ece5421f42c67fa725059ba4129b0ca829b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Juli=C3=A1n=20D=2E=20Ot=C3=A1lvaro?=
 <juliandavid347@gmail.com>
Date: Thu, 12 Feb 2026 09:03:44 +0000
Subject: [PATCH 1/7] feat: add BestDose optimization functionality

Add BestDose dose optimization feature ported from Pmetrics_rust bestdose branch:

- R/PM_bestdose.R: PM_bestdose and PM_bestdose_problem R6 classes for
  Bayesian dose optimization with support for concentration and AUC targets
- src/rust/src/bestdose_executor.rs: Rust backend for BestDose optimization
  using pmcore's BestDoseProblem with ODE model support
- Updated lib.rs with bestdose, bestdose_prepare, bestdose_optimize exports
- Updated extendr-wrappers.R with R-side wrapper functions
- Updated NAMESPACE with PM_bestdose, PM_bestdose_problem, bestdose exports
- Bumped pmcore dependency from 0.21.1 to 0.22.1 (required for bestdose)
- Added libloading dependency for dynamic model loading
- Added bestdose example data (past, prior, target CSVs) and test script
- Fixed executor.rs mutability issue for pmcore 0.22.1 compatibility
---
 .gitignore                                   |   2 +
 Cargo.lock                                   | 713 +++++++------------
 NAMESPACE                                    |   3 +
 R/PM_bestdose.R                              | 381 ++++++++++
 R/extendr-wrappers.R                         |   8 +
 inst/Examples/Rscript/bestdose_simple_test.R |  50 ++
 inst/Examples/src/bestdose_past.csv          |   9 +
 inst/Examples/src/bestdose_prior.csv         |  47 ++
 inst/Examples/src/bestdose_target.csv        |   9 +
 man/PM_bestdose.Rd                           | 172 +++++
 man/PM_bestdose_problem.Rd                   |  88 +++
 man/bestdose.Rd                              |  23 +
 src/rust/Cargo.toml                          |   3 +-
 src/rust/src/bestdose_executor.rs            | 426 +++++++++++
 src/rust/src/executor.rs                     |   2 +-
 src/rust/src/lib.rs                          | 110 +++
 16 files changed, 1608 insertions(+), 438 deletions(-)
 create mode 100644 R/PM_bestdose.R
 create mode 100644 inst/Examples/Rscript/bestdose_simple_test.R
 create mode 100644 inst/Examples/src/bestdose_past.csv
 create mode 100644 inst/Examples/src/bestdose_prior.csv
 create mode 100644 inst/Examples/src/bestdose_target.csv
 create mode 100644 man/PM_bestdose.Rd
 create mode 100644 man/PM_bestdose_problem.Rd
 create mode 100644 man/bestdose.Rd
 create mode 100644 src/rust/src/bestdose_executor.rs

diff --git a/.gitignore b/.gitignore
index ed5327c28..36556e27a 100755
--- a/.gitignore
+++ b/.gitignore
@@ -63,6 +63,7 @@ inst/Learn/_freeze/
 data-raw/Runs/
 tests/testthat/Runs/
 Examples
+!inst/Examples/
 Experimental
 other/
 1
@@ -71,6 +72,7 @@ other/
 
 # Misc project outputs
 *_test.R
+!inst/Examples/Rscript/*_test.R
 docs/
 docs
 src/rust/vendor
diff --git a/Cargo.lock b/Cargo.lock
index 3c31d49bd..bff62dcb9 100755
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -4,21 +4,21 @@ version = 4
 
 [[package]]
 name = "ahash"
-version = "0.8.11"
+version = "0.8.12"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "e89da841a80418a9b391ebaea17f5c112ffaaa96f621d2c285b5174da76b9011"
+checksum = "5a15f179cd60c4584b8a8c596927aadc462e27f2ca70c04e0071964a73ba7a75"
 dependencies = [
  "cfg-if",
  "once_cell",
  "version_check",
- "zerocopy 0.7.35",
+ "zerocopy",
 ]
 
 [[package]]
 name = "aho-corasick"
-version = "1.1.3"
+version = "1.1.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "8e60d3430d3a69478ad0993f19238d2df97c507009a52b3c10addcd7f6bcb916"
+checksum = "ddd31a130427c27518df266943a5308ed92d4b226cc639f5a8f1002816174301"
 dependencies = [
  "memchr",
 ]
@@ -31,9 +31,9 @@ checksum = "683d7910e743518b0e34f1186f92494becacb047c7b6bf616c96772180fef923"
 
 [[package]]
 name = "anyhow"
-version = "1.0.100"
+version = "1.0.101"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "a23eb6b1614318a8071c9b2521f36b424b2c83db5eb3a0fead4a6c0809af6e61"
+checksum = "5f0e0fee31ef5ed1ba1316088939cea399010ed7731dba877ed44aeb407a75ea"
 
 [[package]]
 name = "approx"
@@ -56,7 +56,7 @@ dependencies = [
  "paste",
  "rand 0.9.2",
  "rand_xoshiro",
- "thiserror 2.0.17",
+ "thiserror 2.0.18",
  "web-time",
 ]
 
@@ -72,7 +72,7 @@ dependencies = [
  "num-integer",
  "num-traits",
  "rand 0.9.2",
- "thiserror 2.0.17",
+ "thiserror 2.0.18",
 ]
 
 [[package]]
@@ -87,15 +87,15 @@ dependencies = [
 
 [[package]]
 name = "autocfg"
-version = "1.4.0"
+version = "1.5.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "ace50bade8e6234aa140d9a2f552bbee1db4d353f69b8217bc503490fc1a9f26"
+checksum = "c08606f8c3cbf4ce6ec8e28fb0014a2c086708fe954eaa885384a6165172e7e8"
 
 [[package]]
 name = "bitflags"
-version = "2.9.0"
+version = "2.10.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "5c8214115b7bf84099f1309324e63141d4c5d7cc26862f97a0a857dbefe165bd"
+checksum = "812e12b5285cc515a9c72a5c1d3b6d46a19dac5acfef5265968c166106e31dd3"
 
 [[package]]
 name = "block-buffer"
@@ -108,28 +108,28 @@ dependencies = [
 
 [[package]]
 name = "bumpalo"
-version = "3.17.0"
+version = "3.19.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1628fb46dfa0b37568d12e5edd512553eccf6a22a78e8bde00bb4aed84d5bdbf"
+checksum = "5dd9dc738b7a8311c7ade152424974d8115f2cdad61e8dab8dac9f2362298510"
 
 [[package]]
 name = "bytemuck"
-version = "1.23.2"
+version = "1.25.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "3995eaeebcdf32f91f980d360f78732ddc061097ab4e39991ae7a6ace9194677"
+checksum = "c8efb64bd706a16a1bdde310ae86b351e4d21550d98d056f22f8a7f7a2183fec"
 dependencies = [
  "bytemuck_derive",
 ]
 
 [[package]]
 name = "bytemuck_derive"
-version = "1.10.1"
+version = "1.10.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "4f154e572231cb6ba2bd1176980827e3d5dc04cc183a75dea38109fbdd672d29"
+checksum = "f9abbd1bc6865053c427f7198e6af43bfdedc55ab791faed4fbd361d789575ff"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.106",
+ "syn 2.0.115",
 ]
 
 [[package]]
@@ -147,9 +147,9 @@ dependencies = [
  "ahash",
  "cached_proc_macro",
  "cached_proc_macro_types",
- "hashbrown",
+ "hashbrown 0.15.5",
  "once_cell",
- "thiserror 2.0.17",
+ "thiserror 2.0.18",
  "web-time",
 ]
 
@@ -162,7 +162,7 @@ dependencies = [
  "darling",
  "proc-macro2",
  "quote",
- "syn 2.0.106",
+ "syn 2.0.115",
 ]
 
 [[package]]
@@ -173,18 +173,19 @@ checksum = "ade8366b8bd5ba243f0a58f036cc0ca8a2f069cff1a2351ef1cac6b083e16fc0"
 
 [[package]]
 name = "cc"
-version = "1.2.23"
+version = "1.2.55"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "5f4ac86a9e5bc1e2b3449ab9d7d3a6a405e3d1bb28d7b9be8614f55846ae3766"
+checksum = "47b26a0954ae34af09b50f0de26458fa95369a0d478d8236d3f93082b219bd29"
 dependencies = [
+ "find-msvc-tools",
  "shlex",
 ]
 
 [[package]]
 name = "cfg-if"
-version = "1.0.0"
+version = "1.0.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd"
+checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801"
 
 [[package]]
 name = "cpufeatures"
@@ -253,15 +254,15 @@ checksum = "d0a5c400df2834b80a4c3327b3aad3a4c4cd4de0629063962b03235697506a28"
 
 [[package]]
 name = "crunchy"
-version = "0.2.3"
+version = "0.2.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "43da5946c66ffcc7745f48db692ffbb10a83bfe0afd96235c5c2a4fb23994929"
+checksum = "460fbee9c2c2f33933d720630a6a0bac33ba7053db5344fac858d4b8952d77d5"
 
 [[package]]
 name = "crypto-common"
-version = "0.1.6"
+version = "0.1.7"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1bfb12502f3fc46cca1bb51ac28df9d618d813cdc3d2f25b9fe775a34af26bb3"
+checksum = "78c8292055d1c1df0cce5d180393dc8cce0abec0a7102adb6c7b1eef6016d60a"
 dependencies = [
  "generic-array",
  "typenum",
@@ -269,21 +270,21 @@ dependencies = [
 
 [[package]]
 name = "csv"
-version = "1.3.1"
+version = "1.4.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "acdc4883a9c96732e4733212c01447ebd805833b7275a73ca3ee080fd77afdaf"
+checksum = "52cd9d68cf7efc6ddfaaee42e7288d3a99d613d4b50f76ce9827ae0c6e14f938"
 dependencies = [
  "csv-core",
  "itoa",
  "ryu",
- "serde",
+ "serde_core",
 ]
 
 [[package]]
 name = "csv-core"
-version = "0.1.12"
+version = "0.1.13"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "7d02f3b0da4c6504f86e9cd789d8dbafab48c2321be74e9987593de5a894d93d"
+checksum = "704a3c26996a80471189265814dbc2c257598b96b8a7feae2d31ace646bb9782"
 dependencies = [
  "memchr",
 ]
@@ -309,7 +310,7 @@ dependencies = [
  "proc-macro2",
  "quote",
  "strsim",
- "syn 2.0.106",
+ "syn 2.0.115",
 ]
 
 [[package]]
@@ -320,7 +321,7 @@ checksum = "fc34b93ccb385b40dc71c6fceac4b2ad23662c7eeb248cf10d529b7e055b6ead"
 dependencies = [
  "darling_core",
  "quote",
- "syn 2.0.106",
+ "syn 2.0.115",
 ]
 
 [[package]]
@@ -331,9 +332,9 @@ checksum = "930c7171c8df9fb1782bdf9b918ed9ed2d33d1d22300abb754f9085bc48bf8e8"
 
 [[package]]
 name = "deranged"
-version = "0.4.1"
+version = "0.5.6"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "28cfac68e08048ae1883171632c2aef3ebc555621ae56fbccce1cbf22dd7f058"
+checksum = "cc3dc5ad92c2e2d1c193bbbbdf2ea477cb81331de4f3103f267ca18368b988c4"
 dependencies = [
  "powerfmt",
 ]
@@ -351,7 +352,7 @@ dependencies = [
  "num-traits",
  "petgraph",
  "serde",
- "thiserror 2.0.17",
+ "thiserror 2.0.18",
 ]
 
 [[package]]
@@ -366,13 +367,20 @@ dependencies = [
 
 [[package]]
 name = "dyn-stack"
-version = "0.13.0"
+version = "0.13.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "490bd48eb68fffcfed519b4edbfd82c69cbe741d175b84f0e0cbe8c57cbe0bdd"
+checksum = "1c4713e43e2886ba72b8271aa66c93d722116acf7a75555cce11dcde84388fe8"
 dependencies = [
  "bytemuck",
+ "dyn-stack-macros",
 ]
 
+[[package]]
+name = "dyn-stack-macros"
+version = "0.1.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e1d926b4d407d372f141f93bb444696142c29d32962ccbd3531117cf3aa0bfa9"
+
 [[package]]
 name = "either"
 version = "1.15.0"
@@ -388,7 +396,7 @@ dependencies = [
  "heck",
  "proc-macro2",
  "quote",
- "syn 2.0.106",
+ "syn 2.0.115",
 ]
 
 [[package]]
@@ -417,7 +425,7 @@ checksum = "3bf679796c0322556351f287a51b49e48f7c4986e727b5dd78c972d30e2e16cc"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.106",
+ "syn 2.0.115",
 ]
 
 [[package]]
@@ -428,7 +436,7 @@ checksum = "44f23cf4b44bfce11a86ace86f8a73ffdec849c9fd00a386a53d278bd9e81fb3"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.106",
+ "syn 2.0.115",
 ]
 
 [[package]]
@@ -439,25 +447,31 @@ checksum = "877a4ace8713b0bcf2a4e7eec82529c029f1d0619886d18145fea96c3ffe5c0f"
 
 [[package]]
 name = "extendr-api"
-version = "0.7.1"
+version = "0.8.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "67505d96c7faa49d20e749dba7ba2447db52c40a788fd88cc2b6bef02c02277a"
+checksum = "ea54977c6e37236839ffcbc20b5dcea58aa32ae43fbef54a81e1011dc6b19061"
 dependencies = [
+ "extendr-ffi",
  "extendr-macros",
- "libR-sys",
  "once_cell",
  "paste",
 ]
 
+[[package]]
+name = "extendr-ffi"
+version = "0.8.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c76777174a82bdb3e66872f580687d3d0143eed1df9b9cd72b321b9596a23ca7"
+
 [[package]]
 name = "extendr-macros"
-version = "0.7.1"
+version = "0.8.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "81b58838056f294411d0b2c35ac1a2b24c507d6828b75f2c1e74f00ee9b99267"
+checksum = "661cc4ae29de9c4dafe16cfcbda1dbb9f31bd2568f96ebad232cc1f9bcc8b04d"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.106",
+ "syn 2.0.115",
 ]
 
 [[package]]
@@ -507,7 +521,7 @@ checksum = "2cc4b8cd876795d3b19ddfd59b03faa303c0b8adb9af6e188e81fc647c485bb9"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.106",
+ "syn 2.0.115",
 ]
 
 [[package]]
@@ -528,6 +542,12 @@ dependencies = [
  "reborrow",
 ]
 
+[[package]]
+name = "find-msvc-tools"
+version = "0.1.9"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5baebc0774151f905a1a2cc41989300b1e6fbb29aff0ceffa1064fdd3088d582"
+
 [[package]]
 name = "fixedbitset"
 version = "0.5.7"
@@ -673,16 +693,17 @@ checksum = "5881e4c3c2433fe4905bb19cfd2b5d49d4248274862b68c27c33d9ba4e13f9ec"
 
 [[package]]
 name = "generator"
-version = "0.8.5"
+version = "0.8.8"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "d18470a76cb7f8ff746cf1f7470914f900252ec36bbc40b569d74b1258446827"
+checksum = "52f04ae4152da20c76fe800fa48659201d5cf627c5149ca0b707b69d7eef6cf9"
 dependencies = [
  "cc",
  "cfg-if",
  "libc",
  "log",
  "rustversion",
- "windows",
+ "windows-link",
+ "windows-result",
 ]
 
 [[package]]
@@ -697,25 +718,25 @@ dependencies = [
 
 [[package]]
 name = "getrandom"
-version = "0.2.15"
+version = "0.2.17"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "c4567c8db10ae91089c99af84c68c38da3ec2f087c3f82960bcdbf3656b6f4d7"
+checksum = "ff2abc00be7fca6ebc474524697ae276ad847ad0a6b3faa4bcb027e9a4614ad0"
 dependencies = [
  "cfg-if",
  "libc",
- "wasi 0.11.0+wasi-snapshot-preview1",
+ "wasi",
 ]
 
 [[package]]
 name = "getrandom"
-version = "0.3.2"
+version = "0.3.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "73fea8450eea4bac3940448fb7ae50d91f034f941199fcd9d909a5a07aa455f0"
+checksum = "899def5c37c4fd7b2664648c28120ecec138e4d395b459e5ca34f9cce2dd77fd"
 dependencies = [
  "cfg-if",
  "libc",
  "r-efi",
- "wasi 0.14.2+wasi-0.2.4",
+ "wasip2",
 ]
 
 [[package]]
@@ -810,33 +831,40 @@ checksum = "8babf46d4c1c9d92deac9f7be466f76dfc4482b6452fc5024b5e8daf6ffeb3ee"
 
 [[package]]
 name = "glam"
-version = "0.30.8"
+version = "0.30.10"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "e12d847aeb25f41be4c0ec9587d624e9cd631bc007a8fd7ce3f5851e064c6460"
+checksum = "19fc433e8437a212d1b6f1e68c7824af3aed907da60afa994e7f542d18d12aa9"
 
 [[package]]
 name = "half"
-version = "2.5.0"
+version = "2.7.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "7db2ff139bba50379da6aa0766b52fdcb62cb5b263009b09ed58ba604e14bbd1"
+checksum = "6ea2d84b969582b4b1864a92dc5d27cd2b77b622a8d79306834f1be5ba20d84b"
 dependencies = [
  "bytemuck",
  "cfg-if",
  "crunchy",
  "num-traits",
+ "zerocopy",
 ]
 
 [[package]]
 name = "hashbrown"
-version = "0.15.2"
+version = "0.15.5"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "bf151400ff0baff5465007dd2f3e717f3fe502074ca563069ce3a6629d07b289"
+checksum = "9229cfe53dfd69f0609a49f65461bd93001ea1ef889cd5529dd176593f5338a1"
 dependencies = [
  "allocator-api2",
  "equivalent",
  "foldhash",
 ]
 
+[[package]]
+name = "hashbrown"
+version = "0.16.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "841d1cc9bed7f9236f321df977030373f4a4163ae1a7dbfe1a51a2c1a51d9100"
+
 [[package]]
 name = "heck"
 version = "0.5.0"
@@ -845,9 +873,9 @@ checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea"
 
 [[package]]
 name = "hermit-abi"
-version = "0.3.9"
+version = "0.5.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "d231dfb89cfffdbc30e7fc41579ed6066ad03abda9e567ccafae602b97ec5024"
+checksum = "fc0fef456e4baa96da950455cd02c081ca953b141298e41db3fc7e36b1da849c"
 
 [[package]]
 name = "ident_case"
@@ -857,12 +885,12 @@ checksum = "b9e0384b61958566e926dc50660321d12159025e767c18e043daf26b70104c39"
 
 [[package]]
 name = "indexmap"
-version = "2.8.0"
+version = "2.13.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "3954d50fe15b02142bf25d3b8bdadb634ec3948f103d04ffe3031bc8fe9d7058"
+checksum = "7714e70437a7dc3ac8eb7e6f8df75fd8eb422675fc7678aff7364301092b1017"
 dependencies = [
  "equivalent",
- "hashbrown",
+ "hashbrown 0.16.1",
 ]
 
 [[package]]
@@ -878,15 +906,15 @@ dependencies = [
 
 [[package]]
 name = "itoa"
-version = "1.0.15"
+version = "1.0.17"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "4a5f13b858c8d314ee3e8f639011f7ccefe71f97f96e50151fb991f267928e2c"
+checksum = "92ecc6618181def0457392ccd0ee51198e065e016d1d527a7ac1b6dc7c1f09d2"
 
 [[package]]
 name = "js-sys"
-version = "0.3.77"
+version = "0.3.85"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1cfaf33c695fc6e08064efbc1f72ec937429614f25eef83af942d0e227c3a28f"
+checksum = "8c942ebf8e95485ca0d52d97da7c5a2c387d0e7f0ba4c35e93bfcaee045955b3"
 dependencies = [
  "once_cell",
  "wasm-bindgen",
@@ -898,39 +926,33 @@ version = "1.5.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "bbd2bcb4c963f2ddae06a2efc7e9f3591312473c50c6685e1f298068316e66fe"
 
-[[package]]
-name = "libR-sys"
-version = "0.7.1"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "06ac9752bc1e83f5a354a62b9e81bd8db4468b1008e29f262441e7f0e91e6bb3"
-
 [[package]]
 name = "libc"
-version = "0.2.171"
+version = "0.2.181"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "c19937216e9d3aa9956d9bb8dfc0b0c8beb6058fc4f7a4dc4d850edf86a237d6"
+checksum = "459427e2af2b9c839b132acb702a1c654d95e10f8c326bfc2ad11310e458b1c5"
 
 [[package]]
 name = "libloading"
-version = "0.8.6"
+version = "0.8.9"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "fc2f4eb4bc735547cfed7c0a4922cbd04a4655978c09b54f1f7b228750664c34"
+checksum = "d7c4b02199fee7c5d21a5ae7d8cfa79a6ef5bb2fc834d6e9058e89c825efdc55"
 dependencies = [
  "cfg-if",
- "windows-targets",
+ "windows-link",
 ]
 
 [[package]]
 name = "libm"
-version = "0.2.15"
+version = "0.2.16"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "f9fbbcab51052fe104eb5e5d351cf728d30a5be1fe14d9be8a3b097481fb97de"
+checksum = "b6d2cec3eae94f9f509c767b45932f1ada8350c4bdb85af2fcab4a3c14807981"
 
 [[package]]
 name = "log"
-version = "0.4.27"
+version = "0.4.29"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "13dc2df351e3202783a1fe0d44375f7295ffb4049267b0f3018346dc122a1d94"
+checksum = "5e5032e24019045c762d3c0f28f5b6b8bbf38563a65908389bf7978758920897"
 
 [[package]]
 name = "loom"
@@ -956,9 +978,9 @@ dependencies = [
 
 [[package]]
 name = "matrixmultiply"
-version = "0.3.9"
+version = "0.3.10"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "9380b911e3e96d10c1f415da0876389aaf1b56759054eeb0de7df940c456ba1a"
+checksum = "a06de3016e9fae57a36fd14dba131fccf49f74b40b7fbdb472f96e361ec71a08"
 dependencies = [
  "autocfg",
  "rawpointer",
@@ -966,9 +988,9 @@ dependencies = [
 
 [[package]]
 name = "memchr"
-version = "2.7.4"
+version = "2.8.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "78ca9ab1a0babb1e7d5695e3530886289c18cf2f87ec19a575a0abdce112e3a3"
+checksum = "f8ca58f447f06ed17d5fc4043ce1b10dd205e060fb3ce5b979b8ed8e59ff3f79"
 
 [[package]]
 name = "nalgebra"
@@ -1009,7 +1031,7 @@ dependencies = [
  "glam 0.27.0",
  "glam 0.28.0",
  "glam 0.29.3",
- "glam 0.30.8",
+ "glam 0.30.10",
  "matrixmultiply",
  "nalgebra-macros",
  "num-complex",
@@ -1027,7 +1049,7 @@ checksum = "973e7178a678cfd059ccec50887658d482ce16b0aa9da3888ddeab5cd5eb4889"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.106",
+ "syn 2.0.115",
 ]
 
 [[package]]
@@ -1145,7 +1167,7 @@ version = "0.50.3"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "7957b9740744892f114936ab4a57b3f487491bbeafaf8083688b16841a4240e5"
 dependencies = [
- "windows-sys 0.59.0",
+ "windows-sys 0.61.2",
 ]
 
 [[package]]
@@ -1171,9 +1193,9 @@ dependencies = [
 
 [[package]]
 name = "num-conv"
-version = "0.1.0"
+version = "0.2.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "51d515d32fb182ee37cda2ccdcb92950d6a3c2893aa280e540671c2cd0f3b1d9"
+checksum = "cf97ec579c3c42f953ef76dbf8d55ac91fb219dde70e49aa4a6b7d74e9919050"
 
 [[package]]
 name = "num-integer"
@@ -1207,9 +1229,9 @@ dependencies = [
 
 [[package]]
 name = "num_cpus"
-version = "1.16.0"
+version = "1.17.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "4161fcb6d602d4d2081af7c3a45852d875a03dd337a6bfdd6e06407b61342a43"
+checksum = "91df4bbde75afed763b708b7eee1e8e7651e02d97f6d5dd763e89367e957b23b"
 dependencies = [
  "hermit-abi",
  "libc",
@@ -1229,20 +1251,19 @@ checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a"
 
 [[package]]
 name = "pest"
-version = "2.8.0"
+version = "2.8.6"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "198db74531d58c70a361c42201efde7e2591e976d518caf7662a47dc5720e7b6"
+checksum = "e0848c601009d37dfa3430c4666e147e49cdcf1b92ecd3e63657d8a5f19da662"
 dependencies = [
  "memchr",
- "thiserror 2.0.17",
  "ucd-trie",
 ]
 
 [[package]]
 name = "pest_derive"
-version = "2.8.0"
+version = "2.8.6"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "d725d9cfd79e87dccc9341a2ef39d1b6f6353d68c4b33c177febbe1a402c97c5"
+checksum = "11f486f1ea21e6c10ed15d5a7c77165d0ee443402f0780849d1768e7d9d6fe77"
 dependencies = [
  "pest",
  "pest_generator",
@@ -1250,24 +1271,23 @@ dependencies = [
 
 [[package]]
 name = "pest_generator"
-version = "2.8.0"
+version = "2.8.6"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "db7d01726be8ab66ab32f9df467ae8b1148906685bbe75c82d1e65d7f5b3f841"
+checksum = "8040c4647b13b210a963c1ed407c1ff4fdfa01c31d6d2a098218702e6664f94f"
 dependencies = [
  "pest",
  "pest_meta",
  "proc-macro2",
  "quote",
- "syn 2.0.106",
+ "syn 2.0.115",
 ]
 
 [[package]]
 name = "pest_meta"
-version = "2.8.0"
+version = "2.8.6"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "7f9f832470494906d1fca5329f8ab5791cc60beb230c74815dff541cbd2b5ca0"
+checksum = "89815c69d36021a140146f26659a81d6c2afa33d216d736dd4be5381a7362220"
 dependencies = [
- "once_cell",
  "pest",
  "sha2",
 ]
@@ -1279,16 +1299,16 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "8701b58ea97060d5e5b155d383a69952a60943f0e6dfe30b04c287beb0b27455"
 dependencies = [
  "fixedbitset",
- "hashbrown",
+ "hashbrown 0.15.5",
  "indexmap",
  "serde",
 ]
 
 [[package]]
 name = "pharmsol"
-version = "0.20.0"
+version = "0.21.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "98b8e2ab3a0e91cd4b20c28544cb3676e8df31aa490cf5680ec0531259b5fa4e"
+checksum = "2fc25564d039d0cd5701013aa3785a339b14cf0b51409d7b817320bc360dc944"
 dependencies = [
  "argmin",
  "argmin-math",
@@ -1304,7 +1324,7 @@ dependencies = [
  "serde",
  "serde_json",
  "statrs",
- "thiserror 2.0.17",
+ "thiserror 2.0.18",
  "tracing",
 ]
 
@@ -1320,6 +1340,7 @@ version = "0.1.0"
 dependencies = [
  "anyhow",
  "extendr-api",
+ "libloading",
  "pmcore",
  "rayon",
  "tracing",
@@ -1328,13 +1349,12 @@ dependencies = [
 
 [[package]]
 name = "pmcore"
-version = "0.21.1"
+version = "0.22.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "703e83f4a6a919cc60b85936d560840947b1b07a2d8ccfa7c87144d1722b6d63"
+checksum = "3866100507aa3bcba475381af3102d84b5e503bce82cb82f56cf0fa46cc1e408"
 dependencies = [
  "anyhow",
  "argmin",
- "argmin-math",
  "csv",
  "faer",
  "faer-ext",
@@ -1351,15 +1371,15 @@ dependencies = [
 
 [[package]]
 name = "portable-atomic"
-version = "1.11.0"
+version = "1.13.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "350e9b48cbc6b0e028b0473b114454c6316e57336ee184ceab6e53f72c178b3e"
+checksum = "c33a9471896f1c69cecef8d20cbe2f7accd12527ce60845ff44c153bb2a21b49"
 
 [[package]]
 name = "portable-atomic-util"
-version = "0.2.4"
+version = "0.2.5"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "d8a2f0d8d040d7848a709caf78912debcc3f33ee4b3cac47d73d1e1069e83507"
+checksum = "7a9db96d7fa8782dd8c15ce32ffe8680bbd1e978a43bf51a34d39483540495f5"
 dependencies = [
  "portable-atomic",
 ]
@@ -1376,14 +1396,14 @@ version = "0.2.21"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "85eae3c4ed2f50dcfe72643da4befc30deadb458a9b590d720cde2f2b1e97da9"
 dependencies = [
- "zerocopy 0.8.24",
+ "zerocopy",
 ]
 
 [[package]]
 name = "private-gemm-x86"
-version = "0.1.18"
+version = "0.1.20"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "0b8138b380908e85071bdd6b2841a38b0858ef09848b754a15219d0b9ca90928"
+checksum = "0af8c3e5087969c323f667ccb4b789fa0954f5aa650550e38e81cf9108be21b5"
 dependencies = [
  "crossbeam",
  "defer",
@@ -1397,9 +1417,9 @@ dependencies = [
 
 [[package]]
 name = "proc-macro2"
-version = "1.0.101"
+version = "1.0.106"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "89ae43fd86e4158d6db51ad8e2b80f313af9cc74f5c0e03ccb87de09998732de"
+checksum = "8fd00f0bb2e90d81d1044c2b32617f68fcb9fa3bb7640c23e9c748e53fb30934"
 dependencies = [
  "unicode-ident",
 ]
@@ -1433,9 +1453,9 @@ dependencies = [
 
 [[package]]
 name = "qd"
-version = "0.7.4"
+version = "0.7.7"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "73940173cf92cd24f3650f5f388946524026712a6ca170762340acf5fb3fde0f"
+checksum = "ff8bb755b6008c3b41bf8a0866c8dd4e1245a2f011ceaa22a13ee55c538493e2"
 dependencies = [
  "bytemuck",
  "libm",
@@ -1445,18 +1465,18 @@ dependencies = [
 
 [[package]]
 name = "quote"
-version = "1.0.40"
+version = "1.0.44"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1885c039570dc00dcb4ff087a89e185fd56bae234ddc7f056a945bf36467248d"
+checksum = "21b2ebcf727b7760c461f091f9f0f539b77b8e87f2fd88131e7f1b433b3cece4"
 dependencies = [
  "proc-macro2",
 ]
 
 [[package]]
 name = "r-efi"
-version = "5.2.0"
+version = "5.3.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "74765f6d916ee2faa39bc8e68e4f3ed8949b48cccdac59983d287a7cb71ce9c5"
+checksum = "69cdb34c158ceb288df11e18b4bd39de994f6657d83847bdffdbd7f346754b0f"
 
 [[package]]
 name = "rand"
@@ -1476,7 +1496,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "6db2770f06117d490610c7488547d543617b21bfa07796d7a12f6f1bd53850d1"
 dependencies = [
  "rand_chacha 0.9.0",
- "rand_core 0.9.3",
+ "rand_core 0.9.5",
 ]
 
 [[package]]
@@ -1496,7 +1516,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "d3022b5f1df60f26e1ffddd6c66e8aa15de382ae63b3a0c1bfc0e4d3e3f325cb"
 dependencies = [
  "ppv-lite86",
- "rand_core 0.9.3",
+ "rand_core 0.9.5",
 ]
 
 [[package]]
@@ -1505,16 +1525,16 @@ version = "0.6.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c"
 dependencies = [
- "getrandom 0.2.15",
+ "getrandom 0.2.17",
 ]
 
 [[package]]
 name = "rand_core"
-version = "0.9.3"
+version = "0.9.5"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "99d9a13982dcf210057a8a78572b2217b667c3beacbf3a0d8b454f6f82837d38"
+checksum = "76afc826de14238e6e8c374ddcc1fa19e374fd8dd986b0d2af0d02377261d83c"
 dependencies = [
- "getrandom 0.3.2",
+ "getrandom 0.3.4",
 ]
 
 [[package]]
@@ -1543,14 +1563,14 @@ version = "0.7.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "f703f4665700daf5512dcca5f43afa6af89f09db47fb56be587f80636bda2d41"
 dependencies = [
- "rand_core 0.9.3",
+ "rand_core 0.9.5",
 ]
 
 [[package]]
 name = "raw-cpuid"
-version = "11.5.0"
+version = "11.6.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "c6df7ab838ed27997ba19a4664507e6f82b41fe6e20be42929332156e5e85146"
+checksum = "498cd0dc59d73224351ee52a95fee0f1a617a2eae0e7d9d720cc622c73a54186"
 dependencies = [
  "bitflags",
 ]
@@ -1589,9 +1609,9 @@ checksum = "03251193000f4bd3b042892be858ee50e8b3719f2b08e5833ac4353724632430"
 
 [[package]]
 name = "regex-automata"
-version = "0.4.9"
+version = "0.4.14"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "809e8dc61f6de73b46c85f4c96486310fe304c434cfa43669d7b40f711150908"
+checksum = "6e1dd4122fc1595e8162618945476892eefca7b88c52820e74af6262213cae8f"
 dependencies = [
  "aho-corasick",
  "memchr",
@@ -1600,21 +1620,21 @@ dependencies = [
 
 [[package]]
 name = "regex-syntax"
-version = "0.8.5"
+version = "0.8.9"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "2b15c43186be67a4fd63bee50d0303afffcef381492ebe2c5d87f324e1b8815c"
+checksum = "a96887878f22d7bad8a3b6dc5b7440e0ada9a245242924394987b21cf2210a4c"
 
 [[package]]
 name = "rustversion"
-version = "1.0.20"
+version = "1.0.22"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "eded382c5f5f786b989652c49544c4877d9f015cc22e145a5ea8ea66c2921cd2"
+checksum = "b39cdef0fa800fc44525c84ccb54a029961a8215f9619753635a9c0d2538d46d"
 
 [[package]]
 name = "ryu"
-version = "1.0.20"
+version = "1.0.23"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "28d3b2b1366ec20994f1fd18c3c594f05c5dd4bc44d8bb0c1c632c8d6829481f"
+checksum = "9774ba4a74de5f7b1c1451ed6cd5285a32eddb5cccb8cc655a4e50009e06477f"
 
 [[package]]
 name = "safe_arch"
@@ -1648,41 +1668,52 @@ checksum = "1bc711410fbe7399f390ca1c3b60ad0f53f80e95c5eb935e52268a0e2cd49acc"
 
 [[package]]
 name = "serde"
-version = "1.0.219"
+version = "1.0.228"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9a8e94ea7f378bd32cbbd37198a4a91436180c5bb472411e48b5ec2e2124ae9e"
+dependencies = [
+ "serde_core",
+ "serde_derive",
+]
+
+[[package]]
+name = "serde_core"
+version = "1.0.228"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "5f0e2c6ed6606019b4e29e69dbaba95b11854410e5347d525002456dbbb786b6"
+checksum = "41d385c7d4ca58e59fc732af25c3983b67ac852c1a25000afe1175de458b67ad"
 dependencies = [
  "serde_derive",
 ]
 
 [[package]]
 name = "serde_derive"
-version = "1.0.219"
+version = "1.0.228"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "5b0276cf7f2c73365f7157c8123c21cd9a50fbbd844757af28ca1f5925fc2a00"
+checksum = "d540f220d3187173da220f885ab66608367b6574e925011a9353e4badda91d79"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.106",
+ "syn 2.0.115",
 ]
 
 [[package]]
 name = "serde_json"
-version = "1.0.140"
+version = "1.0.149"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "20068b6e96dc6c9bd23e01df8827e6c7e1f2fddd43c21810382803c136b99373"
+checksum = "83fc039473c5595ace860d8c4fafa220ff474b3fc6bfdb4293327f1a37e94d86"
 dependencies = [
  "itoa",
  "memchr",
- "ryu",
  "serde",
+ "serde_core",
+ "zmij",
 ]
 
 [[package]]
 name = "sha2"
-version = "0.10.8"
+version = "0.10.9"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "793db75ad2bcafc3ffa7c68b215fee268f537982cd901d132f89c6343f3a3dc8"
+checksum = "a7507d819769d01a365ab707794a4084392c824f54a7a6a7862f8c3d0892b283"
 dependencies = [
  "cfg-if",
  "cpufeatures",
@@ -1706,9 +1737,9 @@ checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64"
 
 [[package]]
 name = "simba"
-version = "0.9.0"
+version = "0.9.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "b3a386a501cd104797982c15ae17aafe8b9261315b5d07e3ec803f2ea26be0fa"
+checksum = "c99284beb21666094ba2b75bbceda012e610f5479dfcc2d6e2426f53197ffd95"
 dependencies = [
  "approx",
  "num-complex",
@@ -1719,9 +1750,9 @@ dependencies = [
 
 [[package]]
 name = "smallvec"
-version = "1.14.0"
+version = "1.15.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "7fcf8323ef1faaee30a44a340193b1ac6814fd9b7b4e88e9d4519a3e4abe1cfd"
+checksum = "67b1b7a3b5fe4f1376887184045fcf45c69e92af734b7aaddc05fb777b6fbd03"
 
 [[package]]
 name = "sobol_burley"
@@ -1731,9 +1762,9 @@ checksum = "09f37cae1d97c4078377153ede7a26f7813b689ad5c6b76ff45dc52e53afe1d1"
 
 [[package]]
 name = "spindle"
-version = "0.2.5"
+version = "0.2.6"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "f794dedb367e82477aa6bbf83ea9bbce9bc074b3caacaa82fc4ba398ec9b701d"
+checksum = "673aaca3d8aa5387a6eba861fbf984af5348d9df5d940c25c6366b19556fdf64"
 dependencies = [
  "atomic-wait",
  "crossbeam",
@@ -1773,9 +1804,9 @@ dependencies = [
 
 [[package]]
 name = "syn"
-version = "2.0.106"
+version = "2.0.115"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "ede7c438028d4436d71104916910f5bb611972c5cfd7f89b8300a8186e6fada6"
+checksum = "6e614ed320ac28113fa64972c4262d5dbc89deacdfd00c34a3e4cea073243c12"
 dependencies = [
  "proc-macro2",
  "quote",
@@ -1807,11 +1838,11 @@ dependencies = [
 
 [[package]]
 name = "thiserror"
-version = "2.0.17"
+version = "2.0.18"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "f63587ca0f12b72a0600bcba1d40081f830876000bb46dd2337a3051618f4fc8"
+checksum = "4288b5bcbc7920c07a1149a35cf9590a2aa808e0bc1eafaade0b80947865fbc4"
 dependencies = [
- "thiserror-impl 2.0.17",
+ "thiserror-impl 2.0.18",
 ]
 
 [[package]]
@@ -1822,56 +1853,55 @@ checksum = "4fee6c4efc90059e10f81e6d42c60a18f76588c3d74cb83a0b242a2b6c7504c1"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.106",
+ "syn 2.0.115",
 ]
 
 [[package]]
 name = "thiserror-impl"
-version = "2.0.17"
+version = "2.0.18"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "3ff15c8ecd7de3849db632e14d18d2571fa09dfc5ed93479bc4485c7a517c913"
+checksum = "ebc4ee7f67670e9b64d05fa4253e753e016c6c95ff35b89b7941d6b856dec1d5"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.106",
+ "syn 2.0.115",
 ]
 
 [[package]]
 name = "thread_local"
-version = "1.1.8"
+version = "1.1.9"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "8b9ef9bad013ada3808854ceac7b46812a6465ba368859a37e2100283d2d719c"
+checksum = "f60246a4944f24f6e018aa17cdeffb7818b76356965d03b07d6a9886e8962185"
 dependencies = [
  "cfg-if",
- "once_cell",
 ]
 
 [[package]]
 name = "time"
-version = "0.3.41"
+version = "0.3.47"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "8a7619e19bc266e0f9c5e6686659d394bc57973859340060a69221e57dbc0c40"
+checksum = "743bd48c283afc0388f9b8827b976905fb217ad9e647fae3a379a9283c4def2c"
 dependencies = [
  "deranged",
  "itoa",
  "num-conv",
  "powerfmt",
- "serde",
+ "serde_core",
  "time-core",
  "time-macros",
 ]
 
 [[package]]
 name = "time-core"
-version = "0.1.4"
+version = "0.1.8"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "c9e9a38711f559d9e3ce1cdb06dd7c5b8ea546bc90052da6d06bb76da74bb07c"
+checksum = "7694e1cfe791f8d31026952abf09c69ca6f6fa4e1a1229e18988f06a04a12dca"
 
 [[package]]
 name = "time-macros"
-version = "0.2.22"
+version = "0.2.27"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "3526739392ec93fd8b359c8e98514cb3e8e021beb4e5f597b00a0221f8ed8a49"
+checksum = "2e70e4c5a0e0a8a4823ad65dfe1a6930e4f4d756dcd9dd7939022b5e8c501215"
 dependencies = [
  "num-conv",
  "time-core",
@@ -1879,9 +1909,9 @@ dependencies = [
 
 [[package]]
 name = "tracing"
-version = "0.1.41"
+version = "0.1.44"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "784e0ac535deb450455cbfa28a6f0df145ea1bb7ae51b821cf5e7927fdcfbdd0"
+checksum = "63e71662fa4b2a2c3a26f570f037eb95bb1f85397f3cd8076caed2f026a6d100"
 dependencies = [
  "pin-project-lite",
  "tracing-attributes",
@@ -1890,20 +1920,20 @@ dependencies = [
 
 [[package]]
 name = "tracing-attributes"
-version = "0.1.28"
+version = "0.1.31"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "395ae124c09f9e6918a2310af6038fba074bcf474ac352496d5910dd59a2226d"
+checksum = "7490cfa5ec963746568740651ac6781f701c9c5ea257c58e057f3ba8cf69e8da"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.106",
+ "syn 2.0.115",
 ]
 
 [[package]]
 name = "tracing-core"
-version = "0.1.33"
+version = "0.1.36"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "e672c95779cf947c5311f83787af4fa8fffd12fb27e4993211a84bdfd9610f9c"
+checksum = "db97caf9d906fbde555dd62fa95ddba9eecfd14cb388e4f491a66d74cd5fb79a"
 dependencies = [
  "once_cell",
  "valuable",
@@ -1922,9 +1952,9 @@ dependencies = [
 
 [[package]]
 name = "tracing-subscriber"
-version = "0.3.20"
+version = "0.3.22"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "2054a14f5307d601f88daf0553e1cbf472acc4f2c51afab632431cdcd72124d5"
+checksum = "2f30143827ddab0d256fd843b7a66d164e9f271cfa0dde49142c5ca0ca291f1e"
 dependencies = [
  "matchers",
  "nu-ansi-term",
@@ -1941,9 +1971,9 @@ dependencies = [
 
 [[package]]
 name = "typenum"
-version = "1.18.0"
+version = "1.19.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1dccffe3ce07af9386bfd29e80c0ab1a8205a2fc34e4bcd40364df902cfa8f3f"
+checksum = "562d481066bde0658276a35467c4af00bdc6ee726305698a55b86e61d7ad82bb"
 
 [[package]]
 name = "ucd-trie"
@@ -1953,9 +1983,9 @@ checksum = "2896d95c02a80c6d6a5d6e953d479f5ddf2dfdb6a244441010e373ac0fb88971"
 
 [[package]]
 name = "unicode-ident"
-version = "1.0.18"
+version = "1.0.23"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "5a5f39404a5da50712a4c1eecf25e90dd62b613502b7e925fd4e4d19b5c96512"
+checksum = "537dd038a89878be9b64dd4bd1b260315c1bb94f4d784956b81e27a088d9a09e"
 
 [[package]]
 name = "valuable"
@@ -1981,49 +2011,37 @@ dependencies = [
 
 [[package]]
 name = "wasi"
-version = "0.11.0+wasi-snapshot-preview1"
+version = "0.11.1+wasi-snapshot-preview1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "9c8d87e72b64a3b4db28d11ce29237c246188f4f51057d65a7eab63b7987e423"
+checksum = "ccf3ec651a847eb01de73ccad15eb7d99f80485de043efb2f370cd654f4ea44b"
 
 [[package]]
-name = "wasi"
-version = "0.14.2+wasi-0.2.4"
+name = "wasip2"
+version = "1.0.2+wasi-0.2.9"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "9683f9a5a998d873c0d21fcbe3c083009670149a8fab228644b8bd36b2c48cb3"
+checksum = "9517f9239f02c069db75e65f174b3da828fe5f5b945c4dd26bd25d89c03ebcf5"
 dependencies = [
- "wit-bindgen-rt",
+ "wit-bindgen",
 ]
 
 [[package]]
 name = "wasm-bindgen"
-version = "0.2.100"
+version = "0.2.108"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1edc8929d7499fc4e8f0be2262a241556cfc54a0bea223790e71446f2aab1ef5"
+checksum = "64024a30ec1e37399cf85a7ffefebdb72205ca1c972291c51512360d90bd8566"
 dependencies = [
  "cfg-if",
  "once_cell",
+ "rustversion",
  "wasm-bindgen-macro",
-]
-
-[[package]]
-name = "wasm-bindgen-backend"
-version = "0.2.100"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "2f0a0651a5c2bc21487bde11ee802ccaf4c51935d0d3d42a6101f98161700bc6"
-dependencies = [
- "bumpalo",
- "log",
- "proc-macro2",
- "quote",
- "syn 2.0.106",
  "wasm-bindgen-shared",
 ]
 
 [[package]]
 name = "wasm-bindgen-macro"
-version = "0.2.100"
+version = "0.2.108"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "7fe63fc6d09ed3792bd0897b314f53de8e16568c2b3f7982f468c0bf9bd0b407"
+checksum = "008b239d9c740232e71bd39e8ef6429d27097518b6b30bdf9086833bd5b6d608"
 dependencies = [
  "quote",
  "wasm-bindgen-macro-support",
@@ -2031,22 +2049,22 @@ dependencies = [
 
 [[package]]
 name = "wasm-bindgen-macro-support"
-version = "0.2.100"
+version = "0.2.108"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "8ae87ea40c9f689fc23f209965b6fb8a99ad69aeeb0231408be24920604395de"
+checksum = "5256bae2d58f54820e6490f9839c49780dff84c65aeab9e772f15d5f0e913a55"
 dependencies = [
+ "bumpalo",
  "proc-macro2",
  "quote",
- "syn 2.0.106",
- "wasm-bindgen-backend",
+ "syn 2.0.115",
  "wasm-bindgen-shared",
 ]
 
 [[package]]
 name = "wasm-bindgen-shared"
-version = "0.2.100"
+version = "0.2.108"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1a05d73b933a847d6cccdda8f838a22ff101ad9bf93e33684f39c1f5f0eece3d"
+checksum = "1f01b580c9ac74c8d8f0c0e4afb04eeef2acf145458e52c03845ee9cd23e3d12"
 dependencies = [
  "unicode-ident",
 ]
@@ -2063,9 +2081,9 @@ dependencies = [
 
 [[package]]
 name = "wide"
-version = "0.7.32"
+version = "0.7.33"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "41b5576b9a81633f3e8df296ce0063042a73507636cbe956c61133dd7034ab22"
+checksum = "0ce5da8ecb62bcd8ec8b7ea19f69a51275e91299be594ea5cc6ef7819e16cd03"
 dependencies = [
  "bytemuck",
  "safe_arch",
@@ -2073,111 +2091,24 @@ dependencies = [
 
 [[package]]
 name = "winapi-util"
-version = "0.1.9"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "cf221c93e13a30d793f7645a0e7762c55d169dbb0a49671918a2319d289b10bb"
-dependencies = [
- "windows-sys 0.59.0",
-]
-
-[[package]]
-name = "windows"
-version = "0.61.1"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "c5ee8f3d025738cb02bad7868bbb5f8a6327501e870bf51f1b455b0a2454a419"
-dependencies = [
- "windows-collections",
- "windows-core",
- "windows-future",
- "windows-link",
- "windows-numerics",
-]
-
-[[package]]
-name = "windows-collections"
-version = "0.2.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "3beeceb5e5cfd9eb1d76b381630e82c4241ccd0d27f1a39ed41b2760b255c5e8"
-dependencies = [
- "windows-core",
-]
-
-[[package]]
-name = "windows-core"
-version = "0.61.1"
+version = "0.1.11"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "46ec44dc15085cea82cf9c78f85a9114c463a369786585ad2882d1ff0b0acf40"
+checksum = "c2a7b1c03c876122aa43f3020e6c3c3ee5c05081c9a00739faf7503aeba10d22"
 dependencies = [
- "windows-implement",
- "windows-interface",
- "windows-link",
- "windows-result",
- "windows-strings",
-]
-
-[[package]]
-name = "windows-future"
-version = "0.2.1"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "fc6a41e98427b19fe4b73c550f060b59fa592d7d686537eebf9385621bfbad8e"
-dependencies = [
- "windows-core",
- "windows-link",
- "windows-threading",
-]
-
-[[package]]
-name = "windows-implement"
-version = "0.60.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "a47fddd13af08290e67f4acabf4b459f647552718f683a7b415d290ac744a836"
-dependencies = [
- "proc-macro2",
- "quote",
- "syn 2.0.106",
-]
-
-[[package]]
-name = "windows-interface"
-version = "0.59.1"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "bd9211b69f8dcdfa817bfd14bf1c97c9188afa36f4750130fcdf3f400eca9fa8"
-dependencies = [
- "proc-macro2",
- "quote",
- "syn 2.0.106",
+ "windows-sys 0.61.2",
 ]
 
 [[package]]
 name = "windows-link"
-version = "0.1.1"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "76840935b766e1b0a05c0066835fb9ec80071d4c09a16f6bd5f7e655e3c14c38"
-
-[[package]]
-name = "windows-numerics"
-version = "0.2.0"
+version = "0.2.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "9150af68066c4c5c07ddc0ce30421554771e528bde427614c61038bc2c92c2b1"
-dependencies = [
- "windows-core",
- "windows-link",
-]
+checksum = "f0805222e57f7521d6a62e36fa9163bc891acd422f971defe97d64e70d0a4fe5"
 
 [[package]]
 name = "windows-result"
-version = "0.3.3"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "4b895b5356fc36103d0f64dd1e94dfa7ac5633f1c9dd6e80fe9ec4adef69e09d"
-dependencies = [
- "windows-link",
-]
-
-[[package]]
-name = "windows-strings"
 version = "0.4.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "2a7ab927b2637c19b3dbe0965e75d8f2d30bdd697a1516191cad2ec4df8fb28a"
+checksum = "7781fa89eaf60850ac3d2da7af8e5242a5ea78d1a11c49bf2910bb5a73853eb5"
 dependencies = [
  "windows-link",
 ]
@@ -2188,45 +2119,20 @@ version = "0.42.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "5a3e1820f08b8513f676f7ab6c1f99ff312fb97b553d30ff4dd86f9f15728aa7"
 dependencies = [
- "windows_aarch64_gnullvm 0.42.2",
- "windows_aarch64_msvc 0.42.2",
- "windows_i686_gnu 0.42.2",
- "windows_i686_msvc 0.42.2",
- "windows_x86_64_gnu 0.42.2",
- "windows_x86_64_gnullvm 0.42.2",
- "windows_x86_64_msvc 0.42.2",
+ "windows_aarch64_gnullvm",
+ "windows_aarch64_msvc",
+ "windows_i686_gnu",
+ "windows_i686_msvc",
+ "windows_x86_64_gnu",
+ "windows_x86_64_gnullvm",
+ "windows_x86_64_msvc",
 ]
 
 [[package]]
 name = "windows-sys"
-version = "0.59.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1e38bc4d79ed67fd075bcc251a1c39b32a1776bbe92e5bef1f0bf1f8c531853b"
-dependencies = [
- "windows-targets",
-]
-
-[[package]]
-name = "windows-targets"
-version = "0.52.6"
+version = "0.61.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "9b724f72796e036ab90c1021d4780d4d3d648aca59e491e6b98e725b84e99973"
-dependencies = [
- "windows_aarch64_gnullvm 0.52.6",
- "windows_aarch64_msvc 0.52.6",
- "windows_i686_gnu 0.52.6",
- "windows_i686_gnullvm",
- "windows_i686_msvc 0.52.6",
- "windows_x86_64_gnu 0.52.6",
- "windows_x86_64_gnullvm 0.52.6",
- "windows_x86_64_msvc 0.52.6",
-]
-
-[[package]]
-name = "windows-threading"
-version = "0.1.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "b66463ad2e0ea3bbf808b7f1d371311c80e115c0b71d60efc142cafbcfb057a6"
+checksum = "ae137229bcbd6cdf0f7b80a31df61766145077ddf49416a728b02cb3921ff3fc"
 dependencies = [
  "windows-link",
 ]
@@ -2237,78 +2143,36 @@ version = "0.42.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "597a5118570b68bc08d8d59125332c54f1ba9d9adeedeef5b99b02ba2b0698f8"
 
-[[package]]
-name = "windows_aarch64_gnullvm"
-version = "0.52.6"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "32a4622180e7a0ec044bb555404c800bc9fd9ec262ec147edd5989ccd0c02cd3"
-
 [[package]]
 name = "windows_aarch64_msvc"
 version = "0.42.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "e08e8864a60f06ef0d0ff4ba04124db8b0fb3be5776a5cd47641e942e58c4d43"
 
-[[package]]
-name = "windows_aarch64_msvc"
-version = "0.52.6"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "09ec2a7bb152e2252b53fa7803150007879548bc709c039df7627cabbd05d469"
-
 [[package]]
 name = "windows_i686_gnu"
 version = "0.42.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "c61d927d8da41da96a81f029489353e68739737d3beca43145c8afec9a31a84f"
 
-[[package]]
-name = "windows_i686_gnu"
-version = "0.52.6"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "8e9b5ad5ab802e97eb8e295ac6720e509ee4c243f69d781394014ebfe8bbfa0b"
-
-[[package]]
-name = "windows_i686_gnullvm"
-version = "0.52.6"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "0eee52d38c090b3caa76c563b86c3a4bd71ef1a819287c19d586d7334ae8ed66"
-
 [[package]]
 name = "windows_i686_msvc"
 version = "0.42.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "44d840b6ec649f480a41c8d80f9c65108b92d89345dd94027bfe06ac444d1060"
 
-[[package]]
-name = "windows_i686_msvc"
-version = "0.52.6"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "240948bc05c5e7c6dabba28bf89d89ffce3e303022809e73deaefe4f6ec56c66"
-
 [[package]]
 name = "windows_x86_64_gnu"
 version = "0.42.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "8de912b8b8feb55c064867cf047dda097f92d51efad5b491dfb98f6bbb70cb36"
 
-[[package]]
-name = "windows_x86_64_gnu"
-version = "0.52.6"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "147a5c80aabfbf0c7d901cb5895d1de30ef2907eb21fbbab29ca94c5b08b1a78"
-
 [[package]]
 name = "windows_x86_64_gnullvm"
 version = "0.42.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "26d41b46a36d453748aedef1486d5c7a85db22e56aff34643984ea85514e94a3"
 
-[[package]]
-name = "windows_x86_64_gnullvm"
-version = "0.52.6"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "24d5b23dc417412679681396f2b49f3de8c1473deb516bd34410872eff51ed0d"
-
 [[package]]
 name = "windows_x86_64_msvc"
 version = "0.42.2"
@@ -2316,56 +2180,33 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "9aec5da331524158c6d1a4ac0ab1541149c0b9505fde06423b02f5ef0106b9f0"
 
 [[package]]
-name = "windows_x86_64_msvc"
-version = "0.52.6"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "589f6da84c646204747d1270a2a5661ea66ed1cced2631d546fdfb155959f9ec"
-
-[[package]]
-name = "wit-bindgen-rt"
-version = "0.39.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "6f42320e61fe2cfd34354ecb597f86f413484a798ba44a8ca1165c58d42da6c1"
-dependencies = [
- "bitflags",
-]
-
-[[package]]
-name = "zerocopy"
-version = "0.7.35"
+name = "wit-bindgen"
+version = "0.51.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1b9b4fd18abc82b8136838da5d50bae7bdea537c574d8dc1a34ed098d6c166f0"
-dependencies = [
- "zerocopy-derive 0.7.35",
-]
+checksum = "d7249219f66ced02969388cf2bb044a09756a083d0fab1e566056b04d9fbcaa5"
 
 [[package]]
 name = "zerocopy"
-version = "0.8.24"
+version = "0.8.39"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "2586fea28e186957ef732a5f8b3be2da217d65c5969d4b1e17f973ebbe876879"
+checksum = "db6d35d663eadb6c932438e763b262fe1a70987f9ae936e60158176d710cae4a"
 dependencies = [
- "zerocopy-derive 0.8.24",
+ "zerocopy-derive",
 ]
 
 [[package]]
 name = "zerocopy-derive"
-version = "0.7.35"
+version = "0.8.39"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "fa4f8080344d4671fb4e831a13ad1e68092748387dfc4f55e356242fae12ce3e"
+checksum = "4122cd3169e94605190e77839c9a40d40ed048d305bfdc146e7df40ab0f3e517"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.106",
+ "syn 2.0.115",
 ]
 
 [[package]]
-name = "zerocopy-derive"
-version = "0.8.24"
+name = "zmij"
+version = "1.0.21"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "a996a8f63c5c4448cd959ac1bab0aaa3306ccfd060472f85943ee0750f0169be"
-dependencies = [
- "proc-macro2",
- "quote",
- "syn 2.0.106",
-]
+checksum = "b8848ee67ecc8aedbaf3e4122217aff892639231befc6a1b58d29fff4c2cabaa"
diff --git a/NAMESPACE b/NAMESPACE
index 56f0a9aba..e0d36c1d1 100755
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -45,6 +45,8 @@ export(NPreport)
 export(NPrun)
 export(PMFortranConfig)
 export(PM_batch)
+export(PM_bestdose)
+export(PM_bestdose_problem)
 export(PM_build)
 export(PM_compare)
 export(PM_cov)
@@ -95,6 +97,7 @@ export(add_shapes)
 export(add_smooth)
 export(additive)
 export(all_is_numeric)
+export(bestdose)
 export(build_model)
 export(build_plot)
 export(clear_build)
diff --git a/R/PM_bestdose.R b/R/PM_bestdose.R
new file mode 100644
index 000000000..684c04e02
--- /dev/null
+++ b/R/PM_bestdose.R
@@ -0,0 +1,381 @@
+bestdose_parse_prior <- function(prior) {
+    if (inherits(prior, "PM_result")) {
+        theta_path <- file.path(prior$rundir, "outputs", "theta.csv")
+        if (!file.exists(theta_path)) {
+            cli::cli_abort("theta.csv not found in PM_result outputs")
+        }
+        theta_path
+    } else if (inherits(prior, "PM_final")) {
+        temp_path <- tempfile(fileext = ".csv")
+        bestdose_write_prior_csv(prior, temp_path)
+        temp_path
+    } else if (is.character(prior)) {
+        if (!file.exists(prior)) {
+            cli::cli_abort("Prior file not found: {prior}")
+        }
+        prior
+    } else {
+        cli::cli_abort("prior must be PM_result, PM_final, or path to theta.csv")
+    }
+}
+
+bestdose_write_prior_csv <- function(prior, path) {
+    df <- as.data.frame(prior$popPoints)
+    df$prob <- prior$popProb
+    write.csv(df, path, row.names = FALSE, quote = FALSE)
+}
+
+bestdose_parse_model <- function(model) {
+    if (inherits(model, "PM_model")) {
+        compiled_path <- model$binary_path
+        if (is.null(compiled_path) || !file.exists(compiled_path)) {
+            cli::cli_abort("Model must be compiled first. Use model$compile()")
+        }
+
+        kind <- if (!is.null(model$model_list$analytical) && model$model_list$analytical) {
+            "analytical"
+        } else {
+            "ode"
+        }
+
+        list(path = compiled_path, kind = kind, model = model)
+    } else if (is.character(model)) {
+        if (!file.exists(model)) {
+            cli::cli_abort("Model file not found: {model}")
+        }
+        kind <- if (grepl("analytical", model, ignore.case = TRUE)) {
+            "analytical"
+        } else {
+            "ode"
+        }
+        list(path = model, kind = kind, model = NULL)
+    } else {
+        cli::cli_abort("model must be PM_model or path to compiled model")
+    }
+}
+
+bestdose_parse_data <- function(data) {
+    if (inherits(data, "PM_data")) {
+        temp_path <- tempfile(fileext = ".csv")
+        write.csv(data$standard_data, temp_path, row.names = FALSE, quote = FALSE)
+        temp_path
+    } else if (is.character(data)) {
+        if (!file.exists(data)) {
+            cli::cli_abort("Data file not found: {data}")
+        }
+        data
+    } else {
+        cli::cli_abort("data must be PM_data or path to CSV file")
+    }
+}
+
+bestdose_default_settings <- function(prior, model) {
+    if (inherits(prior, "PM_result")) {
+        return(prior$settings)
+    }
+
+    param_ranges <- lapply(model$model_list$pri, function(x) {
+        c(x$min, x$max)
+    })
+    names(param_ranges) <- tolower(names(param_ranges))
+
+    list(
+        algorithm = "NPAG",
+        ranges = param_ranges,
+        error_models = list(
+            list(
+                initial = 0.0,
+                type = "additive",
+                coeff = c(0.0, 0.2, 0.0, 0.0)
+            )
+        ),
+        max_cycles = 500,
+        points = 2028,
+        seed = 22,
+        prior = "prior.csv",
+        idelta = 0.25,
+        tad = 0.0
+    )
+}
+
+#' @title
+#' Object to contain BestDose optimization results
+#'
+#' @description
+#' `r lifecycle::badge("experimental")`
+#'
+#' This object is created after a successful BestDose optimization run.
+#' BestDose finds optimal dosing regimens to achieve target drug concentrations
+#' or AUC values using Bayesian optimization.
+#'
+#' @export
+PM_bestdose <- R6::R6Class(
+    "PM_bestdose",
+    public = list(
+        result = NULL,
+        problem = NULL,
+        bias_weight = NULL,
+        initialize = function(prior = NULL,
+                              model = NULL,
+                              past_data = NULL,
+                              target = NULL,
+                              dose_range = list(min = 0, max = 1000),
+                              bias_weight = 0.5,
+                              target_type = "concentration",
+                              time_offset = NULL,
+                              settings = NULL,
+                              result = NULL,
+                              problem_obj = NULL,
+                              bias_override = NULL) {
+            if (!is.null(result)) {
+                private$.set_result(result, problem_obj, bias_override)
+                return(invisible(self))
+            }
+
+            if (is.null(target)) {
+                cli::cli_abort("target must be supplied when computing a new BestDose result")
+            }
+
+            problem <- PM_bestdose_problem$new(
+                prior = prior,
+                model = model,
+                past_data = past_data,
+                target = target,
+                dose_range = dose_range,
+                bias_weight = bias_weight,
+                target_type = target_type,
+                time_offset = time_offset,
+                settings = settings
+            )
+
+            raw <- problem$optimize_raw(bias_weight = bias_weight)
+            private$.set_result(raw$result, problem, raw$bias_weight)
+            invisible(self)
+        },
+
+        #' @description
+        #' Print summary of BestDose results
+        print = function() {
+            cat("BestDose Optimization Results\n")
+            cat("==============================\n\n")
+            cat(sprintf("Optimal doses: [%.2f, %.2f] mg\n", self$get_doses()[1], self$get_doses()[2]))
+            cat(sprintf("Objective function: %.10f\n", self$get_objf()))
+            cat(sprintf("ln(Objective): %.4f\n", log(self$get_objf())))
+            cat(sprintf("Method: %s\n", self$get_method()))
+            cat(sprintf("Status: %s\n", self$get_status()))
+            if (!is.null(self$bias_weight)) {
+                cat(sprintf("Bias weight (lambda): %.2f\n", self$bias_weight))
+            }
+            cat(sprintf("\nNumber of predictions: %d\n", nrow(self$result$predictions)))
+            if (!is.null(self$result$auc_predictions)) {
+                cat(sprintf("Number of AUC predictions: %d\n", nrow(self$result$auc_predictions)))
+            }
+            invisible(self)
+        },
+
+        #' @description
+        #' Get optimal dose values
+        #' @return Numeric vector of optimal doses
+        get_doses = function() {
+            self$result$doses
+        },
+
+        #' @description
+        #' Get concentration-time predictions
+        #' @return Data frame with predictions
+        get_predictions = function() {
+            self$result$predictions
+        },
+
+        #' @description
+        #' Get AUC predictions (if available)
+        #' @return Data frame with AUC predictions or NULL
+        get_auc_predictions = function() {
+            self$result$auc_predictions
+        },
+
+        #' @description
+        #' Get objective function value
+        #' @return Numeric objective function value
+        get_objf = function() {
+            self$result$objf
+        },
+
+        #' @description
+        #' Get optimization status
+        #' @return Character string with status
+        get_status = function() {
+            self$result$status
+        },
+
+        #' @description
+        #' Get optimization method used
+        #' @return Character string: "posterior" or "uniform"
+        get_method = function() {
+            self$result$method
+        },
+
+        #' @description
+        #' Save results to RDS file
+        #' @param filename Path to save file. Default: "bestdose_result.rds"
+        save = function(filename = "bestdose_result.rds") {
+            saveRDS(self, filename)
+            cli::cli_alert_success("Results saved to {filename}")
+            invisible(self)
+        }
+    ),
+    private = list(
+        .set_result = function(result, problem, bias_weight) {
+            self$result <- result
+            self$problem <- problem
+            self$bias_weight <- bias_weight
+        }
+    )
+)
+
+#' @title
+#' Prepare a reusable BestDose optimization problem
+#'
+#' @description
+#' `r lifecycle::badge("experimental")`
+#'
+#' Use `PM_bestdose_problem` to mirror the Rust workflow: compute posterior
+#' support points once, inspect them in R, and solve for multiple bias weights
+#' without repeating the expensive initialization step.
+#'
+#' @export
+PM_bestdose_problem <- R6::R6Class(
+    "PM_bestdose_problem",
+    public = list(
+        handle = NULL,
+        theta = NULL,
+        theta_dim = NULL,
+        param_names = NULL,
+        posterior_weights = NULL,
+        population_weights = NULL,
+        bias_weight = NULL,
+        target_type = NULL,
+        dose_range = NULL,
+        model_info = NULL,
+        settings = NULL,
+        initialize = function(prior,
+                              model,
+                              past_data = NULL,
+                              target,
+                              dose_range = list(min = 0, max = 1000),
+                              bias_weight = 0.5,
+                              target_type = "concentration",
+                              time_offset = NULL,
+                              settings = NULL) {
+            if (!target_type %in% c("concentration", "auc_from_zero", "auc_from_last_dose")) {
+                cli::cli_abort("target_type must be one of: concentration, auc_from_zero, auc_from_last_dose")
+            }
+
+            if (bias_weight < 0 || bias_weight > 1) {
+                cli::cli_abort("bias_weight must be between 0 and 1")
+            }
+
+            if (is.null(dose_range$min) || is.null(dose_range$max)) {
+                cli::cli_abort("dose_range must have both 'min' and 'max' elements")
+            }
+
+            if (dose_range$min >= dose_range$max) {
+                cli::cli_abort("dose_range$min must be less than dose_range$max")
+            }
+
+            prior_path <- bestdose_parse_prior(prior)
+            model_info <- bestdose_parse_model(model)
+            past_data_path <- if (!is.null(past_data)) bestdose_parse_data(past_data) else NULL
+            target_data_path <- bestdose_parse_data(target)
+
+            if (is.null(settings)) {
+                model_for_settings <- if (!is.null(model_info$model)) model_info$model else model
+                settings <- bestdose_default_settings(prior, model_for_settings)
+            }
+
+            prep <- bestdose_prepare(
+                model_path = model_info$path,
+                prior_path = prior_path,
+                past_data_path = past_data_path,
+                target_data_path = target_data_path,
+                time_offset = time_offset,
+                dose_min = dose_range$min,
+                dose_max = dose_range$max,
+                bias_weight = bias_weight,
+                target_type = target_type,
+                params = settings,
+                kind = model_info$kind
+            )
+
+            if (is.character(prep)) {
+                cli::cli_abort(prep)
+            }
+
+            dim <- as.integer(prep$theta_dim)
+            theta_matrix <- matrix(prep$theta_values, nrow = dim[1], ncol = dim[2])
+            colnames(theta_matrix) <- prep$param_names
+
+            self$handle <- prep$handle
+            self$theta <- theta_matrix
+            self$theta_dim <- dim
+            self$param_names <- prep$param_names
+            self$posterior_weights <- prep$posterior_weights
+            self$population_weights <- prep$population_weights
+            self$bias_weight <- prep$bias_weight
+            self$target_type <- prep$target_type
+            self$dose_range <- dose_range
+            self$model_info <- model_info
+            self$settings <- settings
+
+            cli::cli_alert_success("BestDose problem prepared with %d support points", dim[1])
+        },
+        finalize = function() {
+            self$handle <- NULL
+        },
+
+        #' @description
+        #' Run optimization and return raw list (doses, objf, predictions)
+        optimize_raw = function(bias_weight = NULL) {
+            private$.run_optimize(bias_weight)
+        },
+
+        #' @description
+        #' Run optimization and return a `PM_bestdose` result object
+        optimize = function(bias_weight = NULL) {
+            raw <- self$optimize_raw(bias_weight)
+            PM_bestdose$new(
+                result = raw$result,
+                problem_obj = self,
+                bias_override = raw$bias_weight
+            )
+        }
+    ),
+    private = list(
+        .run_optimize = function(bias_weight) {
+            if (is.null(self$handle)) {
+                cli::cli_abort("BestDose problem handle has been released")
+            }
+
+            bw <- if (is.null(bias_weight)) self$bias_weight else bias_weight
+
+            if (bw < 0 || bw > 1) {
+                cli::cli_abort("bias_weight must be between 0 and 1")
+            }
+
+            res <- bestdose_optimize(self$handle, bw)
+            if (is.character(res)) {
+                cli::cli_abort(res)
+            }
+
+            list(result = res, bias_weight = bw)
+        }
+    )
+)
+
+#' @export
+PM_bestdose$load <- function(filename = "bestdose_result.rds") {
+    if (!file.exists(filename)) {
+        cli::cli_abort("File not found: {filename}")
+    }
+    readRDS(filename)
+}
diff --git a/R/extendr-wrappers.R b/R/extendr-wrappers.R
index c7c3b3552..9e701cb55 100755
--- a/R/extendr-wrappers.R
+++ b/R/extendr-wrappers.R
@@ -76,5 +76,13 @@ temporary_path <- function() .Call(wrap__temporary_path)
 #' @export
 setup_logs <- function() .Call(wrap__setup_logs)
 
+#' Run BestDose optimization to find optimal doses
+#'@export
+bestdose <- function(model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind) .Call(wrap__bestdose, model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind)
+
+bestdose_prepare <- function(model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind) .Call(wrap__bestdose_prepare, model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind)
+
+bestdose_optimize <- function(handle, bias_weight) .Call(wrap__bestdose_optimize, handle, bias_weight)
+
 
 # nolint end
diff --git a/inst/Examples/Rscript/bestdose_simple_test.R b/inst/Examples/Rscript/bestdose_simple_test.R
new file mode 100644
index 000000000..7cab124fa
--- /dev/null
+++ b/inst/Examples/Rscript/bestdose_simple_test.R
@@ -0,0 +1,50 @@
+library(Pmetrics)
+
+setwd("inst/Examples/Runs")
+
+mod_onecomp <- PM_model$new(
+    pri = list(
+        ke = ab(0.001, 3.0),
+        v = ab(25.0, 250.0)
+    ),
+    eqn = function() {
+        dx[1] <- -ke * X[1] + B[1]
+    },
+    out = function() {
+        Y[1] <- X[1] / v
+    },
+    err = list(
+        additive(1, c(0, 0.20, 0, 0))
+    )
+)
+
+
+past_file <- "../src/bestdose_past.csv"
+target_file <- "../src/bestdose_target.csv"
+prior_file <- "../src/bestdose_prior.csv"
+
+
+# Prepare the problem once (posterior + handle to optimized model)
+problem <- PM_bestdose_problem$new(
+    prior = prior_file,
+    model = mod_onecomp,
+    past_data = past_file,
+    target = target_file,
+    dose_range = list(min = 0, max = 300),
+    bias_weight = 0.0,
+    target_type = "concentration" # "concentration", "auc_from_zero", "auc_from_last_dose"
+)
+
+cat("\nPosterior support points:\n")
+print(head(problem$theta))
+
+# Reuse the same problem for different bias weights
+bias_weights <- seq(0, 1, by = 0.25)
+results <- lapply(bias_weights, function(lambda) {
+    problem$optimize(bias_weight = lambda)
+})
+
+for (i in seq_along(results)) {
+    cat("\n=== Bias weight:", bias_weights[i], "===\n")
+    results[[i]]$print()
+}
diff --git a/inst/Examples/src/bestdose_past.csv b/inst/Examples/src/bestdose_past.csv
new file mode 100644
index 000000000..f31bbc5e9
--- /dev/null
+++ b/inst/Examples/src/bestdose_past.csv
@@ -0,0 +1,9 @@
+"id","evid","time","dur","dose","addl","ii","input","out","outeq","c0","c1","c2","c3"
+1,1,0,0,150,0,0,1,.,.,.,.,.,.
+1,0,2,.,.,.,.,.,0.759050697604428,1,.,.,.,.
+1,0,4,.,.,.,.,.,0.384085169721793,1,.,.,.,.
+1,0,6,.,.,.,.,.,0.194349887386702,1,.,.,.,.
+1,1,12,0,75,0,0,1,.,.,.,.,.,.
+1,0,14,.,.,.,.,.,0.392266577540038,1,.,.,.,.
+1,0,16,.,.,.,.,.,0.198489739204705,1,.,.,.,.
+1,0,18,.,.,.,.,.,0.100437250648841,1,.,.,.,.
diff --git a/inst/Examples/src/bestdose_prior.csv b/inst/Examples/src/bestdose_prior.csv
new file mode 100644
index 000000000..0559dea12
--- /dev/null
+++ b/inst/Examples/src/bestdose_prior.csv
@@ -0,0 +1,47 @@
+ke,v,prob
+0.08736658442020416,104.1576635837555,0.06411826509758521
+0.1305751524925232,97.43967413902283,0.03921568628597478
+0.3540655417442322,86.90321147441864,0.03930537754808413
+0.2908282660484314,113.02810430526733,0.05369391213529181
+0.10566180405616761,140.3564077615738,0.01960784315457728
+0.9795492498397828,221.82963848114014,0.019607843139348134
+0.3214192095041275,68.54407012462616,0.019607843027666164
+0.04363290762901306,85.5460512638092,0.01960721617788871
+0.3493796042442322,75.69715678691864,0.019607073213094923
+0.09591740503311158,70.4919171333313,0.019606975069368423
+0.3318302191734314,92.54958868026733,0.019607488410309255
+0.3091094713449478,92.16600060462952,0.01969419594801783
+0.3319534166574478,118.26951622962952,0.019631823530513178
+0.06250557525157929,107.17031538486481,0.019670607599452994
+0.06885235497951508,76.57066643238068,0.01970716827239974
+0.104416601395607,91.60142242908478,0.019644091781941257
+0.11730292952060699,78.54966461658478,0.019597452684180432
+0.06250557525157929,95.70058882236481,0.019558918237699178
+0.013994081234931946,178.82485330104828,0.019739518261060587
+0.3406021231412888,99.99475717544556,0.019495129816772868
+0.07190197534561157,67.0641827583313,0.019509009848450342
+0.02050767450332642,140.61931788921356,0.019494247900270656
+0.283145404958725,125.85132718086243,0.021689589949674865
+0.28171865940093993,140.98521947860718,0.020065175979272123
+0.2902425238609314,94.48318243026733,0.02049255770567427
+0.09111055810451509,93.09410393238068,0.019720100726323925
+0.2925854926109314,128.45290899276733,0.019127402830562693
+0.2925854926109314,105.60134649276733,0.02235654024528764
+0.3019573676109314,127.61794805526733,0.01954540376390492
+0.044804392004013066,99.6524965763092,0.021082697423336363
+0.28490263152122497,92.62867093086243,0.018708872021706954
+0.08112040762901307,98.0704653263092,0.021075404370655606
+0.036757444095611574,109.3835186958313,0.018114822039674315
+0.09826037378311157,111.2731671333313,0.021727337750962878
+0.09298869409561158,113.6022686958313,0.01713443002609149
+0.3132096666574478,115.89646935462952,0.00474084568678399
+0.07658791284561157,113.3385968208313,0.002945117861059519
+0.08678084223270416,87.5463354587555,0.02933991222454541
+0.08736658442020416,128.0199682712555,0.019952065962974042
+0.05615650746822357,130.38859486579895,0.000015108929159013186
+0.056742249655723574,131.13566517829895,0.012605671648896763
+0.08736658442020416,87.5463354587555,0.007152585492468555
+0.07658791284561157,113.3825421333313,0.03210569535127533
+0.3132096666574478,115.94041466712952,0.03478763091025865
+0.05615650746822357,130.43254017829895,0.04616814785864839
+0.08736658442020416,127.9760229587555,0.00001919610085464287
diff --git a/inst/Examples/src/bestdose_target.csv b/inst/Examples/src/bestdose_target.csv
new file mode 100644
index 000000000..29d8788a6
--- /dev/null
+++ b/inst/Examples/src/bestdose_target.csv
@@ -0,0 +1,9 @@
+"id","evid","time","dur","dose","addl","ii","input","out","outeq","c0","c1","c2","c3"
+1,1,0,0,0,0,0,1,.,.,.,.,.,.
+1,0,2,.,.,.,.,.,0.759050697604428,1,.,.,.,.
+1,0,4,.,.,.,.,.,0.384085169721793,1,.,.,.,.
+1,0,6,.,.,.,.,.,0.194349887386702,1,.,.,.,.
+1,1,12,0,0,0,0,1,.,.,.,.,.,.
+1,0,14,.,.,.,.,.,0.392266577540038,1,.,.,.,.
+1,0,16,.,.,.,.,.,0.198489739204705,1,.,.,.,.
+1,0,18,.,.,.,.,.,0.100437250648841,1,.,.,.,.
diff --git a/man/PM_bestdose.Rd b/man/PM_bestdose.Rd
new file mode 100644
index 000000000..458570517
--- /dev/null
+++ b/man/PM_bestdose.Rd
@@ -0,0 +1,172 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PM_bestdose.R
+\name{PM_bestdose}
+\alias{PM_bestdose}
+\title{Object to contain BestDose optimization results}
+\description{
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}}
+
+This object is created after a successful BestDose optimization run.
+BestDose finds optimal dosing regimens to achieve target drug concentrations
+or AUC values using Bayesian optimization.
+}
+\section{Methods}{
+\subsection{Public methods}{
+\itemize{
+\item \href{#method-PM_bestdose-new}{\code{PM_bestdose$new()}}
+\item \href{#method-PM_bestdose-print}{\code{PM_bestdose$print()}}
+\item \href{#method-PM_bestdose-get_doses}{\code{PM_bestdose$get_doses()}}
+\item \href{#method-PM_bestdose-get_predictions}{\code{PM_bestdose$get_predictions()}}
+\item \href{#method-PM_bestdose-get_auc_predictions}{\code{PM_bestdose$get_auc_predictions()}}
+\item \href{#method-PM_bestdose-get_objf}{\code{PM_bestdose$get_objf()}}
+\item \href{#method-PM_bestdose-get_status}{\code{PM_bestdose$get_status()}}
+\item \href{#method-PM_bestdose-get_method}{\code{PM_bestdose$get_method()}}
+\item \href{#method-PM_bestdose-save}{\code{PM_bestdose$save()}}
+\item \href{#method-PM_bestdose-clone}{\code{PM_bestdose$clone()}}
+}
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-PM_bestdose-new"></a>}}
+\if{latex}{\out{\hypertarget{method-PM_bestdose-new}{}}}
+\subsection{Method \code{new()}}{
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{PM_bestdose$new(
+  prior = NULL,
+  model = NULL,
+  past_data = NULL,
+  target = NULL,
+  dose_range = list(min = 0, max = 1000),
+  bias_weight = 0.5,
+  target_type = "concentration",
+  time_offset = NULL,
+  settings = NULL,
+  result = NULL,
+  problem_obj = NULL,
+  bias_override = NULL
+)}\if{html}{\out{</div>}}
+}
+
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-PM_bestdose-print"></a>}}
+\if{latex}{\out{\hypertarget{method-PM_bestdose-print}{}}}
+\subsection{Method \code{print()}}{
+Print summary of BestDose results
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{PM_bestdose$print()}\if{html}{\out{</div>}}
+}
+
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-PM_bestdose-get_doses"></a>}}
+\if{latex}{\out{\hypertarget{method-PM_bestdose-get_doses}{}}}
+\subsection{Method \code{get_doses()}}{
+Get optimal dose values
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{PM_bestdose$get_doses()}\if{html}{\out{</div>}}
+}
+
+\subsection{Returns}{
+Numeric vector of optimal doses
+}
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-PM_bestdose-get_predictions"></a>}}
+\if{latex}{\out{\hypertarget{method-PM_bestdose-get_predictions}{}}}
+\subsection{Method \code{get_predictions()}}{
+Get concentration-time predictions
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{PM_bestdose$get_predictions()}\if{html}{\out{</div>}}
+}
+
+\subsection{Returns}{
+Data frame with predictions
+}
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-PM_bestdose-get_auc_predictions"></a>}}
+\if{latex}{\out{\hypertarget{method-PM_bestdose-get_auc_predictions}{}}}
+\subsection{Method \code{get_auc_predictions()}}{
+Get AUC predictions (if available)
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{PM_bestdose$get_auc_predictions()}\if{html}{\out{</div>}}
+}
+
+\subsection{Returns}{
+Data frame with AUC predictions or NULL
+}
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-PM_bestdose-get_objf"></a>}}
+\if{latex}{\out{\hypertarget{method-PM_bestdose-get_objf}{}}}
+\subsection{Method \code{get_objf()}}{
+Get objective function value
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{PM_bestdose$get_objf()}\if{html}{\out{</div>}}
+}
+
+\subsection{Returns}{
+Numeric objective function value
+}
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-PM_bestdose-get_status"></a>}}
+\if{latex}{\out{\hypertarget{method-PM_bestdose-get_status}{}}}
+\subsection{Method \code{get_status()}}{
+Get optimization status
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{PM_bestdose$get_status()}\if{html}{\out{</div>}}
+}
+
+\subsection{Returns}{
+Character string with status
+}
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-PM_bestdose-get_method"></a>}}
+\if{latex}{\out{\hypertarget{method-PM_bestdose-get_method}{}}}
+\subsection{Method \code{get_method()}}{
+Get optimization method used
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{PM_bestdose$get_method()}\if{html}{\out{</div>}}
+}
+
+\subsection{Returns}{
+Character string: "posterior" or "uniform"
+}
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-PM_bestdose-save"></a>}}
+\if{latex}{\out{\hypertarget{method-PM_bestdose-save}{}}}
+\subsection{Method \code{save()}}{
+Save results to RDS file
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{PM_bestdose$save(filename = "bestdose_result.rds")}\if{html}{\out{</div>}}
+}
+
+\subsection{Arguments}{
+\if{html}{\out{<div class="arguments">}}
+\describe{
+\item{\code{filename}}{Path to save file. Default: "bestdose_result.rds"}
+}
+\if{html}{\out{</div>}}
+}
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-PM_bestdose-clone"></a>}}
+\if{latex}{\out{\hypertarget{method-PM_bestdose-clone}{}}}
+\subsection{Method \code{clone()}}{
+The objects of this class are cloneable with this method.
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{PM_bestdose$clone(deep = FALSE)}\if{html}{\out{</div>}}
+}
+
+\subsection{Arguments}{
+\if{html}{\out{<div class="arguments">}}
+\describe{
+\item{\code{deep}}{Whether to make a deep clone.}
+}
+\if{html}{\out{</div>}}
+}
+}
+}
diff --git a/man/PM_bestdose_problem.Rd b/man/PM_bestdose_problem.Rd
new file mode 100644
index 000000000..1806c8c2f
--- /dev/null
+++ b/man/PM_bestdose_problem.Rd
@@ -0,0 +1,88 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PM_bestdose.R
+\name{PM_bestdose_problem}
+\alias{PM_bestdose_problem}
+\title{Prepare a reusable BestDose optimization problem}
+\description{
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}}
+
+Use \code{PM_bestdose_problem} to mirror the Rust workflow: compute posterior
+support points once, inspect them in R, and solve for multiple bias weights
+without repeating the expensive initialization step.
+}
+\section{Methods}{
+\subsection{Public methods}{
+\itemize{
+\item \href{#method-PM_bestdose_problem-new}{\code{PM_bestdose_problem$new()}}
+\item \href{#method-PM_bestdose_problem-finalize}{\code{PM_bestdose_problem$finalize()}}
+\item \href{#method-PM_bestdose_problem-optimize_raw}{\code{PM_bestdose_problem$optimize_raw()}}
+\item \href{#method-PM_bestdose_problem-optimize}{\code{PM_bestdose_problem$optimize()}}
+\item \href{#method-PM_bestdose_problem-clone}{\code{PM_bestdose_problem$clone()}}
+}
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-PM_bestdose_problem-new"></a>}}
+\if{latex}{\out{\hypertarget{method-PM_bestdose_problem-new}{}}}
+\subsection{Method \code{new()}}{
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{PM_bestdose_problem$new(
+  prior,
+  model,
+  past_data = NULL,
+  target,
+  dose_range = list(min = 0, max = 1000),
+  bias_weight = 0.5,
+  target_type = "concentration",
+  time_offset = NULL,
+  settings = NULL
+)}\if{html}{\out{</div>}}
+}
+
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-PM_bestdose_problem-finalize"></a>}}
+\if{latex}{\out{\hypertarget{method-PM_bestdose_problem-finalize}{}}}
+\subsection{Method \code{finalize()}}{
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{PM_bestdose_problem$finalize()}\if{html}{\out{</div>}}
+}
+
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-PM_bestdose_problem-optimize_raw"></a>}}
+\if{latex}{\out{\hypertarget{method-PM_bestdose_problem-optimize_raw}{}}}
+\subsection{Method \code{optimize_raw()}}{
+Run optimization and return raw list (doses, objf, predictions)
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{PM_bestdose_problem$optimize_raw(bias_weight = NULL)}\if{html}{\out{</div>}}
+}
+
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-PM_bestdose_problem-optimize"></a>}}
+\if{latex}{\out{\hypertarget{method-PM_bestdose_problem-optimize}{}}}
+\subsection{Method \code{optimize()}}{
+Run optimization and return a \code{PM_bestdose} result object
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{PM_bestdose_problem$optimize(bias_weight = NULL)}\if{html}{\out{</div>}}
+}
+
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-PM_bestdose_problem-clone"></a>}}
+\if{latex}{\out{\hypertarget{method-PM_bestdose_problem-clone}{}}}
+\subsection{Method \code{clone()}}{
+The objects of this class are cloneable with this method.
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{PM_bestdose_problem$clone(deep = FALSE)}\if{html}{\out{</div>}}
+}
+
+\subsection{Arguments}{
+\if{html}{\out{<div class="arguments">}}
+\describe{
+\item{\code{deep}}{Whether to make a deep clone.}
+}
+\if{html}{\out{</div>}}
+}
+}
+}
diff --git a/man/bestdose.Rd b/man/bestdose.Rd
new file mode 100644
index 000000000..93fc8e78b
--- /dev/null
+++ b/man/bestdose.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/extendr-wrappers.R
+\name{bestdose}
+\alias{bestdose}
+\title{Run BestDose optimization to find optimal doses}
+\usage{
+bestdose(
+  model_path,
+  prior_path,
+  past_data_path,
+  target_data_path,
+  time_offset,
+  dose_min,
+  dose_max,
+  bias_weight,
+  target_type,
+  params,
+  kind
+)
+}
+\description{
+Run BestDose optimization to find optimal doses
+}
diff --git a/src/rust/Cargo.toml b/src/rust/Cargo.toml
index 2ee9e7262..431000536 100755
--- a/src/rust/Cargo.toml
+++ b/src/rust/Cargo.toml
@@ -9,8 +9,9 @@ name = 'pm_rs'
 
 [dependencies]
 extendr-api = '*'
-pmcore = {version ="=0.21.1",  features = ["exa"]}
+pmcore = { version = "=0.22.1", features = ["exa"] }
 # pmcore = { path = "../../../PMcore", features = ["exa"] }
+libloading = "0.8"
 
 rayon = "1.10.0"
 anyhow = "1.0.97"
diff --git a/src/rust/src/bestdose_executor.rs b/src/rust/src/bestdose_executor.rs
new file mode 100644
index 000000000..2d83731f6
--- /dev/null
+++ b/src/rust/src/bestdose_executor.rs
@@ -0,0 +1,426 @@
+use crate::{logs::RFormatLayer, settings::settings};
+use extendr_api::prelude::*;
+use pmcore::bestdose::{BestDoseProblem, BestDoseResult, DoseRange, Target};
+use pmcore::prelude::{data, ODE};
+use pmcore::routines::initialization::parse_prior;
+use std::path::PathBuf;
+
+/// Helper to parse target type from string
+pub(crate) fn parse_target_type(target_str: &str) -> std::result::Result<Target, String> {
+    match target_str.to_lowercase().as_str() {
+        "concentration" => Ok(Target::Concentration),
+        "auc_from_zero" | "auc" => Ok(Target::AUCFromZero),
+        "auc_from_last_dose" | "auc_interval" => Ok(Target::AUCFromLastDose),
+        _ => Err(format!(
+            "Invalid target type: {}. Must be 'concentration', 'auc_from_zero', or 'auc_from_last_dose'",
+            target_str
+        )),
+    }
+}
+
+/// R-compatible prediction row for BestDose output
+#[derive(Debug, IntoDataFrameRow)]
+pub struct BestDosePredictionRow {
+    id: String,
+    time: f64,
+    observed: f64,
+    pop_mean: f64,
+    pop_median: f64,
+    post_mean: f64,
+    post_median: f64,
+    outeq: usize,
+}
+
+impl BestDosePredictionRow {
+    pub fn from_np_prediction(
+        pred: &pmcore::routines::output::predictions::NPPredictionRow,
+        id: &str,
+    ) -> Self {
+        Self {
+            id: id.to_string(),
+            time: pred.time(),
+            observed: pred.obs().unwrap_or(0.0),
+            pop_mean: pred.pop_mean(),
+            pop_median: pred.pop_median(),
+            post_mean: pred.post_mean(),
+            post_median: pred.post_median(),
+            outeq: pred.outeq(),
+        }
+    }
+}
+
+/// R-compatible AUC prediction row
+#[derive(Debug, IntoDataFrameRow)]
+pub struct BestDoseAucRow {
+    time: f64,
+    auc: f64,
+}
+
+/// Convert BestDoseResult to R-compatible list structure
+pub(crate) fn convert_bestdose_result_to_r(
+    result: BestDoseResult,
+) -> std::result::Result<Robj, String> {
+    // Extract doses
+    let doses: Vec<f64> = result.doses();
+
+    // Objective function
+    let objf = result.objf();
+
+    // Status
+    let status_str = format!("{:?}", result.status());
+
+    // Predictions as data frame
+    let pred_rows: Vec<BestDosePredictionRow> = result
+        .predictions()
+        .predictions()
+        .iter()
+        .map(|p| BestDosePredictionRow::from_np_prediction(p, "subject_1"))
+        .collect();
+    let pred_df = pred_rows
+        .into_dataframe()
+        .map_err(|e| format!("Failed to create predictions dataframe: {:?}", e))?;
+
+    // AUC predictions (if available)
+    let auc_val = if let Some(auc_preds) = result.auc_predictions() {
+        let auc_rows: Vec<BestDoseAucRow> = auc_preds
+            .iter()
+            .map(|(time, auc)| BestDoseAucRow {
+                time: *time,
+                auc: *auc,
+            })
+            .collect();
+        let auc_df = auc_rows
+            .into_dataframe()
+            .map_err(|e| format!("Failed to create AUC dataframe: {:?}", e))?;
+        Robj::from(auc_df)
+    } else {
+        Robj::from(()) // NULL for no AUC
+    };
+
+    // Optimization method
+    let method_str = format!("{}", result.optimization_method());
+
+    // Build the list using list! macro
+    let output = list!(
+        doses = doses,
+        objf = objf,
+        status = status_str,
+        predictions = pred_df,
+        auc_predictions = auc_val,
+        method = method_str
+    );
+
+    Ok(output.into())
+}
+/// Opaque handle that keeps the dynamic model library alive while reusing the
+/// prepared `BestDoseProblem` for multiple optimization runs.
+pub struct BestDoseProblemHandle {
+    problem: BestDoseProblem,
+    #[allow(dead_code)]
+    library: libloading::Library,
+}
+
+impl BestDoseProblemHandle {
+    #[allow(clippy::too_many_arguments)]
+    pub fn new(
+        model_path: PathBuf,
+        prior_path: PathBuf,
+        past_data_path: Option<PathBuf>,
+        target_data_path: PathBuf,
+        time_offset: Option<f64>,
+        dose_min: f64,
+        dose_max: f64,
+        bias_weight: f64,
+        target_type: &str,
+        params: List,
+    ) -> std::result::Result<Self, String> {
+        let (library, (eq, meta)) =
+            unsafe { pmcore::prelude::pharmsol::exa::load::load::<ODE>(model_path) };
+
+        let settings = settings(params, meta.get_params(), "/tmp/bestdose")
+            .map_err(|e| format!("Failed to parse settings: {}", e))?;
+
+        let (population_theta, prior_weights) =
+            parse_prior(&prior_path.to_str().unwrap().to_string(), &settings)
+                .map_err(|e| format!("Failed to parse prior: {}", e))?;
+
+        let population_weights = prior_weights
+            .ok_or_else(|| "Prior file must contain a 'prob' column with weights".to_string())?;
+
+        let past_data = if let Some(path) = past_data_path {
+            let data = data::read_pmetrics(path.to_str().unwrap())
+                .map_err(|e| format!("Failed to read past data: {}", e))?;
+            let subjects = data.subjects();
+            if subjects.is_empty() {
+                return Err("Past data file contains no subjects".to_string());
+            }
+            Some(subjects[0].clone())
+        } else {
+            None
+        };
+
+        let target_data = {
+            let data = data::read_pmetrics(target_data_path.to_str().unwrap())
+                .map_err(|e| format!("Failed to read target data: {}", e))?;
+            let subjects = data.subjects();
+            if subjects.is_empty() {
+                return Err("Target data file contains no subjects".to_string());
+            }
+            subjects[0].clone()
+        };
+
+        let target_enum = parse_target_type(target_type)?;
+        let doserange = DoseRange::new(dose_min, dose_max);
+
+        let problem = BestDoseProblem::new(
+            &population_theta,
+            &population_weights,
+            past_data,
+            target_data,
+            time_offset,
+            eq,
+            doserange,
+            bias_weight,
+            settings,
+            target_enum,
+        )
+        .map_err(|e| format!("Failed to create BestDose problem: {}", e))?;
+
+        Ok(Self { problem, library })
+    }
+
+    pub fn optimize(
+        &self,
+        bias_weight: Option<f64>,
+    ) -> std::result::Result<BestDoseResult, String> {
+        let configured_problem = match bias_weight {
+            Some(weight) => self.problem.clone().with_bias_weight(weight),
+            None => self.problem.clone(),
+        };
+
+        configured_problem
+            .optimize()
+            .map_err(|e| format!("Optimization failed: {}", e))
+    }
+
+    pub fn problem(&self) -> &BestDoseProblem {
+        &self.problem
+    }
+}
+
+pub(crate) fn bestdose_ode(
+    model_path: PathBuf,
+    prior_path: PathBuf,
+    past_data_path: Option<PathBuf>,
+    target_data_path: PathBuf,
+    time_offset: Option<f64>,
+    dose_min: f64,
+    dose_max: f64,
+    bias_weight: f64,
+    target_type: &str,
+    params: List,
+) -> std::result::Result<BestDoseResult, String> {
+    let handle = BestDoseProblemHandle::new(
+        model_path,
+        prior_path,
+        past_data_path,
+        target_data_path,
+        time_offset,
+        dose_min,
+        dose_max,
+        bias_weight,
+        target_type,
+        params,
+    )?;
+
+    handle.optimize(None)
+}
+
+/// Execute bestdose optimization for analytical models (placeholder - not yet supported)
+pub(crate) fn bestdose_analytical(
+    _model_path: PathBuf,
+    _prior_path: PathBuf,
+    _past_data_path: Option<PathBuf>,
+    _target_data_path: PathBuf,
+    _time_offset: Option<f64>,
+    _dose_min: f64,
+    _dose_max: f64,
+    _bias_weight: f64,
+    _target_type: &str,
+    _params: List,
+) -> std::result::Result<BestDoseResult, String> {
+    Err("BestDose for analytical models is not yet supported".to_string())
+}
+
+pub(crate) struct PosteriorSummary {
+    theta_values: Vec<f64>,
+    theta_dim: (i32, i32),
+    param_names: Vec<String>,
+    posterior_weights: Vec<f64>,
+    population_weights: Vec<f64>,
+    bias_weight: f64,
+    target_type: Target,
+}
+
+fn summarize_problem(problem: &BestDoseProblem) -> PosteriorSummary {
+    let theta = problem.posterior_theta();
+    let matrix = theta.matrix();
+    let nrows = matrix.nrows() as i32;
+    let ncols = matrix.ncols() as i32;
+    let mut theta_values = vec![0.0; (nrows * ncols) as usize];
+
+    for col in 0..ncols as usize {
+        for row in 0..nrows as usize {
+            theta_values[row + col * nrows as usize] = *matrix.get(row, col);
+        }
+    }
+
+    PosteriorSummary {
+        theta_values,
+        theta_dim: (nrows, ncols),
+        param_names: theta.param_names(),
+        posterior_weights: problem.posterior_weights().to_vec(),
+        population_weights: problem.population_weights().to_vec(),
+        bias_weight: problem.bias_weight(),
+        target_type: problem.target_type(),
+    }
+}
+
+fn vec_to_doubles(values: Vec<f64>, label: &str) -> std::result::Result<Doubles, String> {
+    Doubles::try_from(values)
+        .map_err(|e| format!("Failed to convert {} to doubles: {:?}", label, e))
+}
+
+fn dims_to_integers(dim: (i32, i32)) -> std::result::Result<Integers, String> {
+    Integers::try_from(vec![dim.0, dim.1])
+        .map_err(|e| format!("Failed to convert dims to integers: {:?}", e))
+}
+
+fn names_to_strings(names: &[String]) -> Strings {
+    Strings::from_values(names.iter().map(|s| s.as_str()))
+}
+
+pub(crate) fn prepare_bestdose_problem(
+    model_path: PathBuf,
+    prior_path: PathBuf,
+    past_data_path: Option<PathBuf>,
+    target_data_path: PathBuf,
+    time_offset: Option<f64>,
+    dose_min: f64,
+    dose_max: f64,
+    bias_weight: f64,
+    target_type: &str,
+    params: List,
+) -> std::result::Result<(BestDoseProblemHandle, PosteriorSummary), String> {
+    let handle = BestDoseProblemHandle::new(
+        model_path,
+        prior_path,
+        past_data_path,
+        target_data_path,
+        time_offset,
+        dose_min,
+        dose_max,
+        bias_weight,
+        target_type,
+        params,
+    )?;
+
+    let summary = summarize_problem(handle.problem());
+    Ok((handle, summary))
+}
+
+pub(crate) fn bestdose_prepare_internal(
+    model_path: &str,
+    prior_path: &str,
+    past_data_path: Nullable<String>,
+    target_data_path: &str,
+    time_offset: Nullable<f64>,
+    dose_min: f64,
+    dose_max: f64,
+    bias_weight: f64,
+    target_type: &str,
+    params: List,
+    kind: &str,
+) -> Robj {
+    RFormatLayer::reset_global_timer();
+    let _ = crate::setup_logs();
+
+    let past_path = past_data_path.into_option().map(PathBuf::from);
+    let time_offset = time_offset.into_option();
+
+    let preparation = match kind {
+        "ode" => prepare_bestdose_problem(
+            PathBuf::from(model_path),
+            PathBuf::from(prior_path),
+            past_path,
+            PathBuf::from(target_data_path),
+            time_offset,
+            dose_min,
+            dose_max,
+            bias_weight,
+            target_type,
+            params.clone(),
+        ),
+        "analytical" => Err("BestDose for analytical models is not yet supported".to_string()),
+        other => Err(format!("{} is not a supported model type", other)),
+    };
+
+    match preparation {
+        Ok((handle, summary)) => {
+            let theta_values = match vec_to_doubles(summary.theta_values, "theta_values") {
+                Ok(values) => values,
+                Err(e) => return Robj::from(e),
+            };
+            let theta_dim = match dims_to_integers(summary.theta_dim) {
+                Ok(dim) => dim,
+                Err(e) => return Robj::from(e),
+            };
+            let posterior_weights =
+                match vec_to_doubles(summary.posterior_weights, "posterior_weights") {
+                    Ok(values) => values,
+                    Err(e) => return Robj::from(e),
+                };
+            let population_weights =
+                match vec_to_doubles(summary.population_weights, "population_weights") {
+                    Ok(values) => values,
+                    Err(e) => return Robj::from(e),
+                };
+            let param_names = names_to_strings(&summary.param_names);
+            let handle_ptr = ExternalPtr::new(handle);
+
+            let output = list!(
+                handle = handle_ptr,
+                theta_values = theta_values,
+                theta_dim = theta_dim,
+                param_names = param_names,
+                posterior_weights = posterior_weights,
+                population_weights = population_weights,
+                bias_weight = summary.bias_weight,
+                target_type = format!("{:?}", summary.target_type),
+                nspp = summary.theta_dim.0,
+                n_parameters = summary.theta_dim.1
+            );
+
+            output.into()
+        }
+        Err(e) => Robj::from(format!("BestDose prepare failed: {}", e)),
+    }
+}
+
+pub(crate) fn bestdose_optimize_internal(
+    handle: ExternalPtr<BestDoseProblemHandle>,
+    bias_weight: Nullable<f64>,
+) -> Robj {
+    let weight = bias_weight.into_option();
+
+    match handle.try_addr() {
+        Ok(inner) => match inner.optimize(weight) {
+            Ok(result) => match convert_bestdose_result_to_r(result) {
+                Ok(robj) => robj,
+                Err(e) => Robj::from(format!("Failed to convert result: {}", e)),
+            },
+            Err(e) => Robj::from(format!("BestDose optimization failed: {}", e)),
+        },
+        Err(e) => Robj::from(format!("Invalid BestDose handle: {}", e)),
+    }
+}
diff --git a/src/rust/src/executor.rs b/src/rust/src/executor.rs
index 2866673a7..2f22bbabe 100755
--- a/src/rust/src/executor.rs
+++ b/src/rust/src/executor.rs
@@ -42,7 +42,7 @@ pub(crate) fn fit<E: Equation + Send>(
     let data = data::read_pmetrics(data.to_str().unwrap()).expect("Failed to read data");
     //dbg!(&data);
     let mut algorithm = dispatch_algorithm(settings, eq, data)?;
-    let result = algorithm.fit()?;
+    let mut result = algorithm.fit()?;
     result.write_outputs()?;
     Ok(())
 }
diff --git a/src/rust/src/lib.rs b/src/rust/src/lib.rs
index d280ceff1..a1b2475c7 100755
--- a/src/rust/src/lib.rs
+++ b/src/rust/src/lib.rs
@@ -1,4 +1,5 @@
 // mod build;
+mod bestdose_executor;
 mod executor;
 mod logs;
 mod settings;
@@ -8,9 +9,11 @@ use anyhow::Result;
 use extendr_api::prelude::*;
 use pmcore::prelude::{data::read_pmetrics, pharmsol::exa::build, Analytical, ODE};
 use simulation::SimulationRow;
+use std::path::PathBuf;
 use std::process::Command;
 use tracing_subscriber::layer::SubscriberExt;
 
+use crate::bestdose_executor::BestDoseProblemHandle;
 use crate::logs::RFormatLayer;
 
 fn validate_paths(data_path: &str, model_path: &str) {
@@ -291,6 +294,110 @@ fn setup_logs() -> anyhow::Result<()> {
 // Macro to generate exports.
 // This ensures exported functions are registered with R.
 // See corresponding C code in `entrypoint.c`.
+/// Run BestDose optimization to find optimal doses
+///@export
+#[extendr]
+fn bestdose(
+    model_path: &str,
+    prior_path: &str,
+    past_data_path: Nullable<String>,
+    target_data_path: &str,
+    time_offset: Nullable<f64>,
+    dose_min: f64,
+    dose_max: f64,
+    bias_weight: f64,
+    target_type: &str,
+    params: List,
+    kind: &str,
+) -> Robj {
+    RFormatLayer::reset_global_timer();
+    let _ = setup_logs();
+
+    println!("Starting BestDose optimization...");
+
+    let past_path = match past_data_path.into_option() {
+        Some(p) => Some(PathBuf::from(p)),
+        None => None,
+    };
+
+    let time_offset_opt = time_offset.into_option();
+
+    let result = match kind {
+        "ode" => bestdose_executor::bestdose_ode(
+            model_path.into(),
+            prior_path.into(),
+            past_path,
+            target_data_path.into(),
+            time_offset_opt,
+            dose_min,
+            dose_max,
+            bias_weight,
+            target_type,
+            params.clone(),
+        ),
+        "analytical" => bestdose_executor::bestdose_analytical(
+            model_path.into(),
+            prior_path.into(),
+            past_path,
+            target_data_path.into(),
+            time_offset_opt,
+            dose_min,
+            dose_max,
+            bias_weight,
+            target_type,
+            params.clone(),
+        ),
+        _ => {
+            return Robj::from(format!("{} is not a supported model type", kind));
+        }
+    };
+
+    match result {
+        Ok(bd_result) => match bestdose_executor::convert_bestdose_result_to_r(bd_result) {
+            Ok(r) => r,
+            Err(e) => Robj::from(format!("Failed to convert result: {}", e)),
+        },
+        Err(e) => Robj::from(format!("BestDose failed: {}", e)),
+    }
+}
+
+#[extendr]
+fn bestdose_prepare(
+    model_path: &str,
+    prior_path: &str,
+    past_data_path: Nullable<String>,
+    target_data_path: &str,
+    time_offset: Nullable<f64>,
+    dose_min: f64,
+    dose_max: f64,
+    bias_weight: f64,
+    target_type: &str,
+    params: List,
+    kind: &str,
+) -> Robj {
+    bestdose_executor::bestdose_prepare_internal(
+        model_path,
+        prior_path,
+        past_data_path,
+        target_data_path,
+        time_offset,
+        dose_min,
+        dose_max,
+        bias_weight,
+        target_type,
+        params,
+        kind,
+    )
+}
+
+#[extendr]
+fn bestdose_optimize(
+    handle: ExternalPtr<BestDoseProblemHandle>,
+    bias_weight: Nullable<f64>,
+) -> Robj {
+    bestdose_executor::bestdose_optimize_internal(handle, bias_weight)
+}
+
 extendr_module! {
     mod Pmetrics;
     fn simulate_one;
@@ -302,6 +409,9 @@ extendr_module! {
     fn model_parameters;
     fn temporary_path;
     fn setup_logs;
+    fn bestdose;
+    fn bestdose_prepare;
+    fn bestdose_optimize;
 }
 
 // To generate the exported function in R, run the following command:

From 04f741b964001cad6c1aa86ba83d1a0df6982c46 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Juli=C3=A1n=20D=2E=20Ot=C3=A1lvaro?=
 <juliandavid347@gmail.com>
Date: Thu, 12 Feb 2026 09:08:22 +0000
Subject: [PATCH 2/7] fix: use devtools::load_all() in bestdose example script

---
 inst/Examples/Rscript/bestdose_simple_test.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/inst/Examples/Rscript/bestdose_simple_test.R b/inst/Examples/Rscript/bestdose_simple_test.R
index 7cab124fa..58ded6da4 100644
--- a/inst/Examples/Rscript/bestdose_simple_test.R
+++ b/inst/Examples/Rscript/bestdose_simple_test.R
@@ -1,4 +1,4 @@
-library(Pmetrics)
+devtools::load_all()
 
 setwd("inst/Examples/Runs")
 

From 44b0d6aa1ce4c8a48cdfd7e54eccc6a92288490b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Juli=C3=A1n=20D=2E=20Ot=C3=A1lvaro?=
 <juliandavid347@gmail.com>
Date: Wed, 4 Mar 2026 12:25:31 +0000
Subject: [PATCH 3/7] fix: max_cycles are passed correctly, also error model

---
 R/PM_bestdose.R                              | 21 +++++++-------------
 inst/Examples/Rscript/bestdose_simple_test.R |  3 ++-
 2 files changed, 9 insertions(+), 15 deletions(-)

diff --git a/R/PM_bestdose.R b/R/PM_bestdose.R
index 684c04e02..baf87b38e 100644
--- a/R/PM_bestdose.R
+++ b/R/PM_bestdose.R
@@ -69,11 +69,7 @@ bestdose_parse_data <- function(data) {
     }
 }
 
-bestdose_default_settings <- function(prior, model) {
-    if (inherits(prior, "PM_result")) {
-        return(prior$settings)
-    }
-
+bestdose_default_settings <- function(prior, model, max_cycles = 500) {
     param_ranges <- lapply(model$model_list$pri, function(x) {
         c(x$min, x$max)
     })
@@ -82,14 +78,8 @@ bestdose_default_settings <- function(prior, model) {
     list(
         algorithm = "NPAG",
         ranges = param_ranges,
-        error_models = list(
-            list(
-                initial = 0.0,
-                type = "additive",
-                coeff = c(0.0, 0.2, 0.0, 0.0)
-            )
-        ),
-        max_cycles = 500,
+        error_models = lapply(model$model_list$err, function(x) x$flatten()),
+        max_cycles = max_cycles,
         points = 2028,
         seed = 22,
         prior = "prior.csv",
@@ -123,6 +113,7 @@ PM_bestdose <- R6::R6Class(
                               bias_weight = 0.5,
                               target_type = "concentration",
                               time_offset = NULL,
+                              max_cycles = 500,
                               settings = NULL,
                               result = NULL,
                               problem_obj = NULL,
@@ -145,6 +136,7 @@ PM_bestdose <- R6::R6Class(
                 bias_weight = bias_weight,
                 target_type = target_type,
                 time_offset = time_offset,
+                max_cycles = max_cycles,
                 settings = settings
             )
 
@@ -266,6 +258,7 @@ PM_bestdose_problem <- R6::R6Class(
                               bias_weight = 0.5,
                               target_type = "concentration",
                               time_offset = NULL,
+                              max_cycles = 500,
                               settings = NULL) {
             if (!target_type %in% c("concentration", "auc_from_zero", "auc_from_last_dose")) {
                 cli::cli_abort("target_type must be one of: concentration, auc_from_zero, auc_from_last_dose")
@@ -290,7 +283,7 @@ PM_bestdose_problem <- R6::R6Class(
 
             if (is.null(settings)) {
                 model_for_settings <- if (!is.null(model_info$model)) model_info$model else model
-                settings <- bestdose_default_settings(prior, model_for_settings)
+                settings <- bestdose_default_settings(prior, model_for_settings, max_cycles = max_cycles)
             }
 
             prep <- bestdose_prepare(
diff --git a/inst/Examples/Rscript/bestdose_simple_test.R b/inst/Examples/Rscript/bestdose_simple_test.R
index 58ded6da4..6867ce7d5 100644
--- a/inst/Examples/Rscript/bestdose_simple_test.R
+++ b/inst/Examples/Rscript/bestdose_simple_test.R
@@ -32,7 +32,8 @@ problem <- PM_bestdose_problem$new(
     target = target_file,
     dose_range = list(min = 0, max = 300),
     bias_weight = 0.0,
-    target_type = "concentration" # "concentration", "auc_from_zero", "auc_from_last_dose"
+    target_type = "concentration", # "concentration", "auc_from_zero", "auc_from_last_dose"
+    max_cycles = 500
 )
 
 cat("\nPosterior support points:\n")

From 9fc29ddd84d9856d35d9450cdf069b0029c7b3da Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Juli=C3=A1n=20D=2E=20Ot=C3=A1lvaro?=
 <juliandavid347@gmail.com>
Date: Fri, 6 Mar 2026 15:21:02 +0000
Subject: [PATCH 4/7] feat: new bestdose API

---
 Cargo.lock                                   |  16 +-
 NAMESPACE                                    |   2 +-
 R/PM_bestdose.R                              | 123 ++++++++-------
 R/extendr-wrappers.R                         |   4 +-
 inst/Examples/Rscript/bestdose_simple_test.R |  19 +--
 src/rust/Cargo.toml                          |   6 +-
 src/rust/src/bestdose_executor.rs            | 150 +++++++------------
 src/rust/src/lib.rs                          |  33 ++--
 8 files changed, 156 insertions(+), 197 deletions(-)

diff --git a/Cargo.lock b/Cargo.lock
index bff62dcb9..103699091 100755
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -341,9 +341,9 @@ dependencies = [
 
 [[package]]
 name = "diffsol"
-version = "0.7.0"
+version = "0.8.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "fcb280582d6bf3f01ad1bc197c38f0d88793e22e300ce96866381e237f5fc684"
+checksum = "31fd20d3173a897a4a016c7f320b82e81a4bae05f3616c1acd0d9da97d080ee7"
 dependencies = [
  "faer",
  "faer-traits",
@@ -934,9 +934,9 @@ checksum = "459427e2af2b9c839b132acb702a1c654d95e10f8c326bfc2ad11310e458b1c5"
 
 [[package]]
 name = "libloading"
-version = "0.8.9"
+version = "0.9.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "d7c4b02199fee7c5d21a5ae7d8cfa79a6ef5bb2fc834d6e9058e89c825efdc55"
+checksum = "754ca22de805bb5744484a5b151a9e1a8e837d5dc232c2d7d8c2e3492edc8b60"
 dependencies = [
  "cfg-if",
  "windows-link",
@@ -1306,9 +1306,9 @@ dependencies = [
 
 [[package]]
 name = "pharmsol"
-version = "0.21.0"
+version = "0.22.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "2fc25564d039d0cd5701013aa3785a339b14cf0b51409d7b817320bc360dc944"
+checksum = "8fbfd50c78d8848657bd38effdbaf404c8308e01d6956ba345a774f98ad7d938"
 dependencies = [
  "argmin",
  "argmin-math",
@@ -1349,9 +1349,7 @@ dependencies = [
 
 [[package]]
 name = "pmcore"
-version = "0.22.1"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "3866100507aa3bcba475381af3102d84b5e503bce82cb82f56cf0fa46cc1e408"
+version = "0.22.2"
 dependencies = [
  "anyhow",
  "argmin",
diff --git a/NAMESPACE b/NAMESPACE
index e0d36c1d1..62480dfb4 100755
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -46,7 +46,7 @@ export(NPrun)
 export(PMFortranConfig)
 export(PM_batch)
 export(PM_bestdose)
-export(PM_bestdose_problem)
+export(PM_bestdose_posterior)
 export(PM_build)
 export(PM_compare)
 export(PM_cov)
diff --git a/R/PM_bestdose.R b/R/PM_bestdose.R
index baf87b38e..72ca569a6 100644
--- a/R/PM_bestdose.R
+++ b/R/PM_bestdose.R
@@ -103,7 +103,7 @@ PM_bestdose <- R6::R6Class(
     "PM_bestdose",
     public = list(
         result = NULL,
-        problem = NULL,
+        posterior = NULL,
         bias_weight = NULL,
         initialize = function(prior = NULL,
                               model = NULL,
@@ -116,10 +116,10 @@ PM_bestdose <- R6::R6Class(
                               max_cycles = 500,
                               settings = NULL,
                               result = NULL,
-                              problem_obj = NULL,
+                              posterior_obj = NULL,
                               bias_override = NULL) {
             if (!is.null(result)) {
-                private$.set_result(result, problem_obj, bias_override)
+                private$.set_result(result, posterior_obj, bias_override)
                 return(invisible(self))
             }
 
@@ -127,21 +127,22 @@ PM_bestdose <- R6::R6Class(
                 cli::cli_abort("target must be supplied when computing a new BestDose result")
             }
 
-            problem <- PM_bestdose_problem$new(
+            posterior <- PM_bestdose_posterior$new(
                 prior = prior,
                 model = model,
                 past_data = past_data,
+                max_cycles = max_cycles,
+                settings = settings
+            )
+
+            raw <- posterior$optimize_raw(
                 target = target,
                 dose_range = dose_range,
                 bias_weight = bias_weight,
                 target_type = target_type,
-                time_offset = time_offset,
-                max_cycles = max_cycles,
-                settings = settings
+                time_offset = time_offset
             )
-
-            raw <- problem$optimize_raw(bias_weight = bias_weight)
-            private$.set_result(raw$result, problem, raw$bias_weight)
+            private$.set_result(raw$result, posterior, raw$bias_weight)
             invisible(self)
         },
 
@@ -217,27 +218,27 @@ PM_bestdose <- R6::R6Class(
         }
     ),
     private = list(
-        .set_result = function(result, problem, bias_weight) {
+        .set_result = function(result, posterior, bias_weight) {
             self$result <- result
-            self$problem <- problem
+            self$posterior <- posterior
             self$bias_weight <- bias_weight
         }
     )
 )
 
 #' @title
-#' Prepare a reusable BestDose optimization problem
+#' Compute a reusable BestDose posterior
 #'
 #' @description
 #' `r lifecycle::badge("experimental")`
 #'
-#' Use `PM_bestdose_problem` to mirror the Rust workflow: compute posterior
-#' support points once, inspect them in R, and solve for multiple bias weights
-#' without repeating the expensive initialization step.
+#' Use `PM_bestdose_posterior` to compute the Bayesian posterior once from
+#' prior population data and patient history, then call `$optimize()` multiple
+#' times with different targets, dose ranges, or bias weights.
 #'
 #' @export
-PM_bestdose_problem <- R6::R6Class(
-    "PM_bestdose_problem",
+PM_bestdose_posterior <- R6::R6Class(
+    "PM_bestdose_posterior",
     public = list(
         handle = NULL,
         theta = NULL,
@@ -245,41 +246,16 @@ PM_bestdose_problem <- R6::R6Class(
         param_names = NULL,
         posterior_weights = NULL,
         population_weights = NULL,
-        bias_weight = NULL,
-        target_type = NULL,
-        dose_range = NULL,
         model_info = NULL,
         settings = NULL,
         initialize = function(prior,
                               model,
                               past_data = NULL,
-                              target,
-                              dose_range = list(min = 0, max = 1000),
-                              bias_weight = 0.5,
-                              target_type = "concentration",
-                              time_offset = NULL,
                               max_cycles = 500,
                               settings = NULL) {
-            if (!target_type %in% c("concentration", "auc_from_zero", "auc_from_last_dose")) {
-                cli::cli_abort("target_type must be one of: concentration, auc_from_zero, auc_from_last_dose")
-            }
-
-            if (bias_weight < 0 || bias_weight > 1) {
-                cli::cli_abort("bias_weight must be between 0 and 1")
-            }
-
-            if (is.null(dose_range$min) || is.null(dose_range$max)) {
-                cli::cli_abort("dose_range must have both 'min' and 'max' elements")
-            }
-
-            if (dose_range$min >= dose_range$max) {
-                cli::cli_abort("dose_range$min must be less than dose_range$max")
-            }
-
             prior_path <- bestdose_parse_prior(prior)
             model_info <- bestdose_parse_model(model)
             past_data_path <- if (!is.null(past_data)) bestdose_parse_data(past_data) else NULL
-            target_data_path <- bestdose_parse_data(target)
 
             if (is.null(settings)) {
                 model_for_settings <- if (!is.null(model_info$model)) model_info$model else model
@@ -290,12 +266,6 @@ PM_bestdose_problem <- R6::R6Class(
                 model_path = model_info$path,
                 prior_path = prior_path,
                 past_data_path = past_data_path,
-                target_data_path = target_data_path,
-                time_offset = time_offset,
-                dose_min = dose_range$min,
-                dose_max = dose_range$max,
-                bias_weight = bias_weight,
-                target_type = target_type,
                 params = settings,
                 kind = model_info$kind
             )
@@ -314,13 +284,10 @@ PM_bestdose_problem <- R6::R6Class(
             self$param_names <- prep$param_names
             self$posterior_weights <- prep$posterior_weights
             self$population_weights <- prep$population_weights
-            self$bias_weight <- prep$bias_weight
-            self$target_type <- prep$target_type
-            self$dose_range <- dose_range
             self$model_info <- model_info
             self$settings <- settings
 
-            cli::cli_alert_success("BestDose problem prepared with %d support points", dim[1])
+            cli::cli_alert_success("BestDose posterior computed with {dim[1]} support points")
         },
         finalize = function() {
             self$handle <- NULL
@@ -328,39 +295,67 @@ PM_bestdose_problem <- R6::R6Class(
 
         #' @description
         #' Run optimization and return raw list (doses, objf, predictions)
-        optimize_raw = function(bias_weight = NULL) {
-            private$.run_optimize(bias_weight)
+        optimize_raw = function(target,
+                                dose_range = list(min = 0, max = 1000),
+                                bias_weight = 0.5,
+                                target_type = "concentration",
+                                time_offset = NULL) {
+            private$.run_optimize(target, dose_range, bias_weight, target_type, time_offset)
         },
 
         #' @description
         #' Run optimization and return a `PM_bestdose` result object
-        optimize = function(bias_weight = NULL) {
-            raw <- self$optimize_raw(bias_weight)
+        optimize = function(target,
+                            dose_range = list(min = 0, max = 1000),
+                            bias_weight = 0.5,
+                            target_type = "concentration",
+                            time_offset = NULL) {
+            raw <- self$optimize_raw(target, dose_range, bias_weight, target_type, time_offset)
             PM_bestdose$new(
                 result = raw$result,
-                problem_obj = self,
+                posterior_obj = self,
                 bias_override = raw$bias_weight
             )
         }
     ),
     private = list(
-        .run_optimize = function(bias_weight) {
+        .run_optimize = function(target, dose_range, bias_weight, target_type, time_offset) {
             if (is.null(self$handle)) {
-                cli::cli_abort("BestDose problem handle has been released")
+                cli::cli_abort("BestDose posterior handle has been released")
             }
 
-            bw <- if (is.null(bias_weight)) self$bias_weight else bias_weight
+            if (!target_type %in% c("concentration", "auc_from_zero", "auc_from_last_dose")) {
+                cli::cli_abort("target_type must be one of: concentration, auc_from_zero, auc_from_last_dose")
+            }
 
-            if (bw < 0 || bw > 1) {
+            if (bias_weight < 0 || bias_weight > 1) {
                 cli::cli_abort("bias_weight must be between 0 and 1")
             }
 
-            res <- bestdose_optimize(self$handle, bw)
+            if (is.null(dose_range$min) || is.null(dose_range$max)) {
+                cli::cli_abort("dose_range must have both 'min' and 'max' elements")
+            }
+
+            if (dose_range$min >= dose_range$max) {
+                cli::cli_abort("dose_range$min must be less than dose_range$max")
+            }
+
+            target_data_path <- bestdose_parse_data(target)
+
+            res <- bestdose_optimize(
+                self$handle,
+                target_data_path,
+                time_offset,
+                dose_range$min,
+                dose_range$max,
+                bias_weight,
+                target_type
+            )
             if (is.character(res)) {
                 cli::cli_abort(res)
             }
 
-            list(result = res, bias_weight = bw)
+            list(result = res, bias_weight = bias_weight)
         }
     )
 )
diff --git a/R/extendr-wrappers.R b/R/extendr-wrappers.R
index 9e701cb55..9f6077b39 100755
--- a/R/extendr-wrappers.R
+++ b/R/extendr-wrappers.R
@@ -80,9 +80,9 @@ setup_logs <- function() .Call(wrap__setup_logs)
 #'@export
 bestdose <- function(model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind) .Call(wrap__bestdose, model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind)
 
-bestdose_prepare <- function(model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind) .Call(wrap__bestdose_prepare, model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind)
+bestdose_prepare <- function(model_path, prior_path, past_data_path, params, kind) .Call(wrap__bestdose_prepare, model_path, prior_path, past_data_path, params, kind)
 
-bestdose_optimize <- function(handle, bias_weight) .Call(wrap__bestdose_optimize, handle, bias_weight)
+bestdose_optimize <- function(handle, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type) .Call(wrap__bestdose_optimize, handle, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type)
 
 
 # nolint end
diff --git a/inst/Examples/Rscript/bestdose_simple_test.R b/inst/Examples/Rscript/bestdose_simple_test.R
index 6867ce7d5..f05a7b15d 100644
--- a/inst/Examples/Rscript/bestdose_simple_test.R
+++ b/inst/Examples/Rscript/bestdose_simple_test.R
@@ -24,25 +24,26 @@ target_file <- "../src/bestdose_target.csv"
 prior_file <- "../src/bestdose_prior.csv"
 
 
-# Prepare the problem once (posterior + handle to optimized model)
-problem <- PM_bestdose_problem$new(
+# Step 1: Compute the posterior once (expensive step)
+posterior <- PM_bestdose_posterior$new(
     prior = prior_file,
     model = mod_onecomp,
     past_data = past_file,
-    target = target_file,
-    dose_range = list(min = 0, max = 300),
-    bias_weight = 0.0,
-    target_type = "concentration", # "concentration", "auc_from_zero", "auc_from_last_dose"
     max_cycles = 500
 )
 
 cat("\nPosterior support points:\n")
-print(head(problem$theta))
+print(head(posterior$theta))
 
-# Reuse the same problem for different bias weights
+# Step 2: Reuse the posterior for different bias weights
 bias_weights <- seq(0, 1, by = 0.25)
 results <- lapply(bias_weights, function(lambda) {
-    problem$optimize(bias_weight = lambda)
+    posterior$optimize(
+        target = target_file,
+        dose_range = list(min = 0, max = 300),
+        bias_weight = lambda,
+        target_type = "concentration"
+    )
 })
 
 for (i in seq_along(results)) {
diff --git a/src/rust/Cargo.toml b/src/rust/Cargo.toml
index 431000536..dfbb32e67 100755
--- a/src/rust/Cargo.toml
+++ b/src/rust/Cargo.toml
@@ -9,9 +9,9 @@ name = 'pm_rs'
 
 [dependencies]
 extendr-api = '*'
-pmcore = { version = "=0.22.1", features = ["exa"] }
-# pmcore = { path = "../../../PMcore", features = ["exa"] }
-libloading = "0.8"
+pmcore = { path = "../../../PMcore", features = ["exa"] }
+# pmcore = { version = "=0.22.1", features = ["exa"] }
+libloading = "0.9"
 
 rayon = "1.10.0"
 anyhow = "1.0.97"
diff --git a/src/rust/src/bestdose_executor.rs b/src/rust/src/bestdose_executor.rs
index 2d83731f6..aafee6871 100644
--- a/src/rust/src/bestdose_executor.rs
+++ b/src/rust/src/bestdose_executor.rs
@@ -1,6 +1,6 @@
 use crate::{logs::RFormatLayer, settings::settings};
 use extendr_api::prelude::*;
-use pmcore::bestdose::{BestDoseProblem, BestDoseResult, DoseRange, Target};
+use pmcore::bestdose::{BestDosePosterior, BestDoseResult, DoseRange, Target};
 use pmcore::prelude::{data, ODE};
 use pmcore::routines::initialization::parse_prior;
 use std::path::PathBuf;
@@ -113,25 +113,18 @@ pub(crate) fn convert_bestdose_result_to_r(
     Ok(output.into())
 }
 /// Opaque handle that keeps the dynamic model library alive while reusing the
-/// prepared `BestDoseProblem` for multiple optimization runs.
-pub struct BestDoseProblemHandle {
-    problem: BestDoseProblem,
+/// prepared `BestDosePosterior` for multiple optimization runs.
+pub struct BestDosePosteriorHandle {
+    posterior: BestDosePosterior,
     #[allow(dead_code)]
     library: libloading::Library,
 }
 
-impl BestDoseProblemHandle {
-    #[allow(clippy::too_many_arguments)]
+impl BestDosePosteriorHandle {
     pub fn new(
         model_path: PathBuf,
         prior_path: PathBuf,
         past_data_path: Option<PathBuf>,
-        target_data_path: PathBuf,
-        time_offset: Option<f64>,
-        dose_min: f64,
-        dose_max: f64,
-        bias_weight: f64,
-        target_type: &str,
         params: List,
     ) -> std::result::Result<Self, String> {
         let (library, (eq, meta)) =
@@ -159,52 +152,48 @@ impl BestDoseProblemHandle {
             None
         };
 
-        let target_data = {
-            let data = data::read_pmetrics(target_data_path.to_str().unwrap())
-                .map_err(|e| format!("Failed to read target data: {}", e))?;
-            let subjects = data.subjects();
-            if subjects.is_empty() {
-                return Err("Target data file contains no subjects".to_string());
-            }
-            subjects[0].clone()
-        };
-
-        let target_enum = parse_target_type(target_type)?;
-        let doserange = DoseRange::new(dose_min, dose_max);
-
-        let problem = BestDoseProblem::new(
+        let posterior = BestDosePosterior::compute(
             &population_theta,
             &population_weights,
             past_data,
-            target_data,
-            time_offset,
             eq,
-            doserange,
-            bias_weight,
             settings,
-            target_enum,
         )
-        .map_err(|e| format!("Failed to create BestDose problem: {}", e))?;
+        .map_err(|e| format!("Failed to compute BestDose posterior: {}", e))?;
 
-        Ok(Self { problem, library })
+        Ok(Self { posterior, library })
     }
 
+    #[allow(clippy::too_many_arguments)]
     pub fn optimize(
         &self,
-        bias_weight: Option<f64>,
+        target_data_path: PathBuf,
+        time_offset: Option<f64>,
+        dose_min: f64,
+        dose_max: f64,
+        bias_weight: f64,
+        target_type: &str,
     ) -> std::result::Result<BestDoseResult, String> {
-        let configured_problem = match bias_weight {
-            Some(weight) => self.problem.clone().with_bias_weight(weight),
-            None => self.problem.clone(),
+        let target_data = {
+            let data = data::read_pmetrics(target_data_path.to_str().unwrap())
+                .map_err(|e| format!("Failed to read target data: {}", e))?;
+            let subjects = data.subjects();
+            if subjects.is_empty() {
+                return Err("Target data file contains no subjects".to_string());
+            }
+            subjects[0].clone()
         };
 
-        configured_problem
-            .optimize()
+        let target_enum = parse_target_type(target_type)?;
+        let dose_range = DoseRange::new(dose_min, dose_max);
+
+        self.posterior
+            .optimize(target_data, time_offset, dose_range, bias_weight, target_enum)
             .map_err(|e| format!("Optimization failed: {}", e))
     }
 
-    pub fn problem(&self) -> &BestDoseProblem {
-        &self.problem
+    pub fn posterior(&self) -> &BestDosePosterior {
+        &self.posterior
     }
 }
 
@@ -220,20 +209,14 @@ pub(crate) fn bestdose_ode(
     target_type: &str,
     params: List,
 ) -> std::result::Result<BestDoseResult, String> {
-    let handle = BestDoseProblemHandle::new(
+    let handle = BestDosePosteriorHandle::new(
         model_path,
         prior_path,
         past_data_path,
-        target_data_path,
-        time_offset,
-        dose_min,
-        dose_max,
-        bias_weight,
-        target_type,
         params,
     )?;
 
-    handle.optimize(None)
+    handle.optimize(target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type)
 }
 
 /// Execute bestdose optimization for analytical models (placeholder - not yet supported)
@@ -258,12 +241,10 @@ pub(crate) struct PosteriorSummary {
     param_names: Vec<String>,
     posterior_weights: Vec<f64>,
     population_weights: Vec<f64>,
-    bias_weight: f64,
-    target_type: Target,
 }
 
-fn summarize_problem(problem: &BestDoseProblem) -> PosteriorSummary {
-    let theta = problem.posterior_theta();
+fn summarize_handle(handle: &BestDosePosteriorHandle) -> PosteriorSummary {
+    let theta = handle.posterior().theta();
     let matrix = theta.matrix();
     let nrows = matrix.nrows() as i32;
     let ncols = matrix.ncols() as i32;
@@ -279,10 +260,8 @@ fn summarize_problem(problem: &BestDoseProblem) -> PosteriorSummary {
         theta_values,
         theta_dim: (nrows, ncols),
         param_names: theta.param_names(),
-        posterior_weights: problem.posterior_weights().to_vec(),
-        population_weights: problem.population_weights().to_vec(),
-        bias_weight: problem.bias_weight(),
-        target_type: problem.target_type(),
+        posterior_weights: handle.posterior().posterior_weights().to_vec(),
+        population_weights: handle.posterior().population_weights().to_vec(),
     }
 }
 
@@ -300,32 +279,20 @@ fn names_to_strings(names: &[String]) -> Strings {
     Strings::from_values(names.iter().map(|s| s.as_str()))
 }
 
-pub(crate) fn prepare_bestdose_problem(
+pub(crate) fn prepare_bestdose_posterior(
     model_path: PathBuf,
     prior_path: PathBuf,
     past_data_path: Option<PathBuf>,
-    target_data_path: PathBuf,
-    time_offset: Option<f64>,
-    dose_min: f64,
-    dose_max: f64,
-    bias_weight: f64,
-    target_type: &str,
     params: List,
-) -> std::result::Result<(BestDoseProblemHandle, PosteriorSummary), String> {
-    let handle = BestDoseProblemHandle::new(
+) -> std::result::Result<(BestDosePosteriorHandle, PosteriorSummary), String> {
+    let handle = BestDosePosteriorHandle::new(
         model_path,
         prior_path,
         past_data_path,
-        target_data_path,
-        time_offset,
-        dose_min,
-        dose_max,
-        bias_weight,
-        target_type,
         params,
     )?;
 
-    let summary = summarize_problem(handle.problem());
+    let summary = summarize_handle(&handle);
     Ok((handle, summary))
 }
 
@@ -333,12 +300,6 @@ pub(crate) fn bestdose_prepare_internal(
     model_path: &str,
     prior_path: &str,
     past_data_path: Nullable<String>,
-    target_data_path: &str,
-    time_offset: Nullable<f64>,
-    dose_min: f64,
-    dose_max: f64,
-    bias_weight: f64,
-    target_type: &str,
     params: List,
     kind: &str,
 ) -> Robj {
@@ -346,19 +307,12 @@ pub(crate) fn bestdose_prepare_internal(
     let _ = crate::setup_logs();
 
     let past_path = past_data_path.into_option().map(PathBuf::from);
-    let time_offset = time_offset.into_option();
 
     let preparation = match kind {
-        "ode" => prepare_bestdose_problem(
+        "ode" => prepare_bestdose_posterior(
             PathBuf::from(model_path),
             PathBuf::from(prior_path),
             past_path,
-            PathBuf::from(target_data_path),
-            time_offset,
-            dose_min,
-            dose_max,
-            bias_weight,
-            target_type,
             params.clone(),
         ),
         "analytical" => Err("BestDose for analytical models is not yet supported".to_string()),
@@ -395,8 +349,6 @@ pub(crate) fn bestdose_prepare_internal(
                 param_names = param_names,
                 posterior_weights = posterior_weights,
                 population_weights = population_weights,
-                bias_weight = summary.bias_weight,
-                target_type = format!("{:?}", summary.target_type),
                 nspp = summary.theta_dim.0,
                 n_parameters = summary.theta_dim.1
             );
@@ -408,13 +360,25 @@ pub(crate) fn bestdose_prepare_internal(
 }
 
 pub(crate) fn bestdose_optimize_internal(
-    handle: ExternalPtr<BestDoseProblemHandle>,
-    bias_weight: Nullable<f64>,
+    handle: ExternalPtr<BestDosePosteriorHandle>,
+    target_data_path: &str,
+    time_offset: Nullable<f64>,
+    dose_min: f64,
+    dose_max: f64,
+    bias_weight: f64,
+    target_type: &str,
 ) -> Robj {
-    let weight = bias_weight.into_option();
+    let time_offset = time_offset.into_option();
 
     match handle.try_addr() {
-        Ok(inner) => match inner.optimize(weight) {
+        Ok(inner) => match inner.optimize(
+            PathBuf::from(target_data_path),
+            time_offset,
+            dose_min,
+            dose_max,
+            bias_weight,
+            target_type,
+        ) {
             Ok(result) => match convert_bestdose_result_to_r(result) {
                 Ok(robj) => robj,
                 Err(e) => Robj::from(format!("Failed to convert result: {}", e)),
diff --git a/src/rust/src/lib.rs b/src/rust/src/lib.rs
index a1b2475c7..1f63e9916 100755
--- a/src/rust/src/lib.rs
+++ b/src/rust/src/lib.rs
@@ -13,7 +13,7 @@ use std::path::PathBuf;
 use std::process::Command;
 use tracing_subscriber::layer::SubscriberExt;
 
-use crate::bestdose_executor::BestDoseProblemHandle;
+use crate::bestdose_executor::BestDosePosteriorHandle;
 use crate::logs::RFormatLayer;
 
 fn validate_paths(data_path: &str, model_path: &str) {
@@ -366,12 +366,6 @@ fn bestdose_prepare(
     model_path: &str,
     prior_path: &str,
     past_data_path: Nullable<String>,
-    target_data_path: &str,
-    time_offset: Nullable<f64>,
-    dose_min: f64,
-    dose_max: f64,
-    bias_weight: f64,
-    target_type: &str,
     params: List,
     kind: &str,
 ) -> Robj {
@@ -379,12 +373,6 @@ fn bestdose_prepare(
         model_path,
         prior_path,
         past_data_path,
-        target_data_path,
-        time_offset,
-        dose_min,
-        dose_max,
-        bias_weight,
-        target_type,
         params,
         kind,
     )
@@ -392,10 +380,23 @@ fn bestdose_prepare(
 
 #[extendr]
 fn bestdose_optimize(
-    handle: ExternalPtr<BestDoseProblemHandle>,
-    bias_weight: Nullable<f64>,
+    handle: ExternalPtr<BestDosePosteriorHandle>,
+    target_data_path: &str,
+    time_offset: Nullable<f64>,
+    dose_min: f64,
+    dose_max: f64,
+    bias_weight: f64,
+    target_type: &str,
 ) -> Robj {
-    bestdose_executor::bestdose_optimize_internal(handle, bias_weight)
+    bestdose_executor::bestdose_optimize_internal(
+        handle,
+        target_data_path,
+        time_offset,
+        dose_min,
+        dose_max,
+        bias_weight,
+        target_type,
+    )
 }
 
 extendr_module! {

From f909cc4d7397bb37f2fe4c91c1f7fbd0d9bd785d Mon Sep 17 00:00:00 2001
From: Michael Neely <mneely@usc.edu>
Date: Fri, 6 Mar 2026 09:50:17 -0800
Subject: [PATCH 5/7] chore: remove old functions

---
 R/PM_bestdose.R | 22 +++++++---------------
 1 file changed, 7 insertions(+), 15 deletions(-)

diff --git a/R/PM_bestdose.R b/R/PM_bestdose.R
index 72ca569a6..bb18eef6b 100644
--- a/R/PM_bestdose.R
+++ b/R/PM_bestdose.R
@@ -293,15 +293,6 @@ PM_bestdose_posterior <- R6::R6Class(
             self$handle <- NULL
         },
 
-        #' @description
-        #' Run optimization and return raw list (doses, objf, predictions)
-        optimize_raw = function(target,
-                                dose_range = list(min = 0, max = 1000),
-                                bias_weight = 0.5,
-                                target_type = "concentration",
-                                time_offset = NULL) {
-            private$.run_optimize(target, dose_range, bias_weight, target_type, time_offset)
-        },
 
         #' @description
         #' Run optimization and return a `PM_bestdose` result object
@@ -310,7 +301,7 @@ PM_bestdose_posterior <- R6::R6Class(
                             bias_weight = 0.5,
                             target_type = "concentration",
                             time_offset = NULL) {
-            raw <- self$optimize_raw(target, dose_range, bias_weight, target_type, time_offset)
+            raw <- private$.run_optimize(target, dose_range, bias_weight, target_type, time_offset)
             PM_bestdose$new(
                 result = raw$result,
                 posterior_obj = self,
@@ -321,23 +312,24 @@ PM_bestdose_posterior <- R6::R6Class(
     private = list(
         .run_optimize = function(target, dose_range, bias_weight, target_type, time_offset) {
             if (is.null(self$handle)) {
-                cli::cli_abort("BestDose posterior handle has been released")
+                cli::cli_abort(c( x = "PM_bestdose_posterior object is not properly initialized",
+                "i" = "Re-run `PM_bestdose$new()`."   ))
             }
 
             if (!target_type %in% c("concentration", "auc_from_zero", "auc_from_last_dose")) {
-                cli::cli_abort("target_type must be one of: concentration, auc_from_zero, auc_from_last_dose")
+                cli::cli_abort(" {.arg target_type} must be one of: concentration, auc_from_zero, auc_from_last_dose")
             }
 
             if (bias_weight < 0 || bias_weight > 1) {
-                cli::cli_abort("bias_weight must be between 0 and 1")
+                cli::cli_abort("{.arg bias_weight} must be between 0 and 1")
             }
 
             if (is.null(dose_range$min) || is.null(dose_range$max)) {
-                cli::cli_abort("dose_range must have both 'min' and 'max' elements")
+                cli::cli_abort("{.arg dose_range} must have both 'min' and 'max' elements")
             }
 
             if (dose_range$min >= dose_range$max) {
-                cli::cli_abort("dose_range$min must be less than dose_range$max")
+                cli::cli_abort("{.arg dose_range$min} must be less than {.arg dose_range$max}")
             }
 
             target_data_path <- bestdose_parse_data(target)

From 7c8e43344ff0ff0ce3190affd7ef7eb71c3b0d32 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Juli=C3=A1n=20D=2E=20Ot=C3=A1lvaro?=
 <juliandavid347@gmail.com>
Date: Fri, 6 Mar 2026 18:13:27 +0000
Subject: [PATCH 6/7] fix: using PMcore on gh

---
 Cargo.lock          | 1 +
 src/rust/Cargo.toml | 3 ++-
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/Cargo.lock b/Cargo.lock
index 103699091..37e4fcd97 100755
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -1350,6 +1350,7 @@ dependencies = [
 [[package]]
 name = "pmcore"
 version = "0.22.2"
+source = "git+https://github.com/LAPKB/PMcore?branch=feat%2Fnew-bestdose-api#1e4c6ca2fa34a91f4d9ac95fb898b67fb8e68e7f"
 dependencies = [
  "anyhow",
  "argmin",
diff --git a/src/rust/Cargo.toml b/src/rust/Cargo.toml
index dfbb32e67..ba495ddc1 100755
--- a/src/rust/Cargo.toml
+++ b/src/rust/Cargo.toml
@@ -9,7 +9,8 @@ name = 'pm_rs'
 
 [dependencies]
 extendr-api = '*'
-pmcore = { path = "../../../PMcore", features = ["exa"] }
+# pmcore = { path = "../../../PMcore", features = ["exa"] }
+pmcore = { git = "https://github.com/LAPKB/PMcore", branch = "feat/new-bestdose-api", features = ["exa"] }
 # pmcore = { version = "=0.22.1", features = ["exa"] }
 libloading = "0.9"
 

From 2ebf493344f3d2ffb2b3fae8714deb26a2af5128 Mon Sep 17 00:00:00 2001
From: Michael Neely <mneely@usc.edu>
Date: Thu, 9 Apr 2026 06:53:10 -0700
Subject: [PATCH 7/7] fix: update pkgdown

---
 DESCRIPTION  | 4 ++--
 _pkgdown.yml | 3 +++
 2 files changed, 5 insertions(+), 2 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index eb0d6da70..a6e3184d3 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -16,8 +16,8 @@ Authors@R: c(
   person("Robert", "Leary", email = "", role = "ctb")
     )
 Version: 3.0.7
-URL: https://lapkb.github.io/Pmetrics_rust/
-BugReports: https://github.com/LAPKB/Pmetrics_rust/issues
+URL: https://lapkb.github.io/Pmetrics
+BugReports: https://github.com/LAPKB/Pmetrics/issues
 SystemRequirements: Cargo (>= 1.82) (Rust's package manager), rustc
 Depends: 
     R (>= 4.2)
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 47770eb84..8efa331d4 100755
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -113,6 +113,9 @@ reference:
 - title: Utility
   desc: Utility functions
 - contents:
+  - bestdose
+  - PM_bestdose
+  - PM_bestdose_problem
   - PM_build
   - PM_help
   - PMcheck