diff --git a/Cargo.toml b/Cargo.toml
index 671c2d81..bee8c913 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -43,6 +43,8 @@ nalgebra = "0.32.5"
 faer = "0.18.2"
 faer-ext = {version = "0.1.0", features = ["nalgebra", "ndarray"]}
 serde_with = "3.8.1"
+anyhow = "1.0.83"
+num-traits = "0.2.19"
 
 
 [profile.release]
diff --git a/benches/compare.rs b/benches/compare.rs
index 48d4bcd2..854dae09 100644
--- a/benches/compare.rs
+++ b/benches/compare.rs
@@ -1,4 +1,4 @@
-use pmcore::routines::data::{parse_pmetrics::read_pmetrics, DataTrait};
+use pmcore::prelude::data::{parse_pmetrics::read_pmetrics, DataTrait};
 use pmcore::simulator::analytical::one_compartment_with_absorption;
 use pmcore::{prelude::*, simulator::Equation};
 use std::path::Path;
diff --git a/examples/new_simulator.rs b/examples/new_simulator.rs
index e969574d..2b87a9ed 100644
--- a/examples/new_simulator.rs
+++ b/examples/new_simulator.rs
@@ -1,5 +1,5 @@
+use pmcore::prelude::data::{parse_pmetrics::read_pmetrics, DataTrait};
 use pmcore::prelude::*;
-use pmcore::routines::data::{parse_pmetrics::read_pmetrics, DataTrait};
 use pmcore::simulator::analytical::one_compartment_with_absorption;
 use pmcore::simulator::Equation;
 use std::path::Path;
diff --git a/examples/parse_pmetrics.rs b/examples/parse_pmetrics.rs
index 110866bd..d1a41e75 100644
--- a/examples/parse_pmetrics.rs
+++ b/examples/parse_pmetrics.rs
@@ -1,4 +1,4 @@
-use pmcore::routines::data::*;
+use pmcore::prelude::data::*;
 
 fn main() {
     let path = std::path::Path::new("examples/data/bimodal_ke.csv");
diff --git a/src/algorithms/npag.rs b/src/algorithms/npag.rs
index a50c26fc..bf9d85ef 100644
--- a/src/algorithms/npag.rs
+++ b/src/algorithms/npag.rs
@@ -6,6 +6,7 @@ use crate::{
         output::{CycleLog, NPCycle, NPResult},
         qr,
         settings::Settings,
+        simulator::get_population_predictions,
     },
     routines::expansion::adaptative_grid::adaptative_grid,
     simulator::{likelihood::PopulationPredictions, Equation},
@@ -16,7 +17,7 @@ use ndarray::{Array, Array1, Array2, Axis};
 use ndarray_stats::{DeviationExt, QuantileExt};
 use tokio::sync::mpsc::UnboundedSender;
 
-use super::{data::Subject, get_population_predictions};
+use super::data::Subject;
 
 const THETA_E: f64 = 1e-4; // Convergence criteria
 const THETA_G: f64 = 1e-4; // Objective function convergence criteria
diff --git a/src/algorithms/npod.rs b/src/algorithms/npod.rs
index 772b1b3d..3b50136c 100644
--- a/src/algorithms/npod.rs
+++ b/src/algorithms/npod.rs
@@ -8,6 +8,7 @@ use crate::{
         output::{CycleLog, NPCycle, NPResult},
         qr,
         settings::Settings,
+        simulator::get_population_predictions,
     },
     simulator::{likelihood::PopulationPredictions, Equation},
     tui::ui::Comm,
@@ -17,7 +18,7 @@ use ndarray::{Array, Array1, Array2, Axis};
 use ndarray_stats::{DeviationExt, QuantileExt};
 use tokio::sync::mpsc::UnboundedSender;
 
-use super::{data::Subject, get_population_predictions};
+use super::data::Subject;
 
 const THETA_D: f64 = 1e-4;
 const THETA_F: f64 = 1e-2;
diff --git a/src/algorithms/postprob.rs b/src/algorithms/postprob.rs
index 680c11bd..4081d818 100644
--- a/src/algorithms/postprob.rs
+++ b/src/algorithms/postprob.rs
@@ -5,6 +5,7 @@ use crate::{
         ipm_faer::burke,
         output::NPResult,
         settings::Settings,
+        simulator::get_population_predictions,
     },
     simulator::Equation,
     tui::ui::Comm,
@@ -13,7 +14,7 @@ use crate::{
 use ndarray::{Array1, Array2};
 use tokio::sync::mpsc::UnboundedSender;
 
-use super::{data::Subject, get_population_predictions};
+use super::data::Subject;
 
 /// Posterior probability algorithm
 /// Reweights the prior probabilities to the observed data and error model
diff --git a/src/routines/data/mod.rs b/src/data/mod.rs
similarity index 100%
rename from src/routines/data/mod.rs
rename to src/data/mod.rs
diff --git a/src/routines/data/parse_pmetrics.rs b/src/data/parse_pmetrics.rs
similarity index 99%
rename from src/routines/data/parse_pmetrics.rs
rename to src/data/parse_pmetrics.rs
index cf7982e5..fd66868d 100644
--- a/src/routines/data/parse_pmetrics.rs
+++ b/src/data/parse_pmetrics.rs
@@ -1,4 +1,4 @@
-use crate::routines::data::*;
+use crate::prelude::data::*;
 use serde::de::{MapAccess, Visitor};
 use serde::{de, Deserialize, Deserializer, Serialize};
 use std::collections::HashMap;
diff --git a/src/lib.rs b/src/lib.rs
index 4f849abb..c928d239 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -52,10 +52,12 @@ pub mod algorithms;
 
 pub mod simulator;
 
+/// New data format
+pub mod data;
+
 /// Routines for the crate
 pub mod routines {
-    /// New data format
-    pub mod data;
+
     /// Routines for initializing the grid
     pub mod initialization;
     pub mod optimization {
@@ -84,9 +86,6 @@ pub mod routines {
         pub mod qr;
         pub mod sigma;
     }
-    pub mod simulation {
-        pub mod predict;
-    }
 }
 
 /// Entry points for external use of the framework.
@@ -112,12 +111,9 @@ pub mod prelude {
     pub use crate::logger;
     pub use crate::prelude::evaluation::{sigma, *};
     pub use crate::routines::condensation;
-    pub use crate::routines::data::CovariateTrait;
     pub use crate::routines::expansion::*;
     pub use crate::routines::initialization::*;
     pub use crate::routines::optimization;
-    pub use crate::routines::simulation::predict::*;
-    pub use crate::routines::simulation::*;
     pub use crate::routines::*;
     pub use crate::tui::ui::*;
     #[macro_export]
diff --git a/src/routines/optimization/d_optimizer.rs b/src/routines/optimization/d_optimizer.rs
index 481d3e54..c5d03c43 100644
--- a/src/routines/optimization/d_optimizer.rs
+++ b/src/routines/optimization/d_optimizer.rs
@@ -5,10 +5,8 @@ use argmin::{
 use ndarray::{Array1, Axis};
 
 use crate::{
-    routines::{
-        data::Subject, evaluation::sigma::ErrorModel,
-        simulation::predict::get_population_predictions,
-    },
+    prelude::{data::Subject, simulator::get_population_predictions},
+    routines::evaluation::sigma::ErrorModel,
     simulator::Equation,
 };
 
diff --git a/src/routines/simulation/predict.rs b/src/routines/simulation/predict.rs
deleted file mode 100644
index bacac2d0..00000000
--- a/src/routines/simulation/predict.rs
+++ /dev/null
@@ -1,182 +0,0 @@
-use crate::routines::data::Subject;
-
-// use crate::routines::evaluation::sigma::ErrorModel;
-use crate::simulator::likelihood::PopulationPredictions;
-// use crate::simulator::likelihood::Prediction;
-use crate::simulator::Equation;
-// use dashmap::mapref::entry::Entry;
-// use dashmap::DashMap;
-// use lazy_static::lazy_static;
-use ndarray::parallel::prelude::*;
-use ndarray::prelude::*;
-// use ndarray::Array1;
-use ndarray::{Array2, Axis};
-// use std::collections::HashMap;
-// use std::error;
-// use std::hash::{Hash, Hasher};
-
-/// Number of support points to cache for each scenario
-// const CACHE_SIZE: usize = 1000;
-
-// #[derive(Debug, Clone)]
-// pub struct Model {
-//     params: HashMap<String, f64>,
-//     _scenario: Scenario,
-//     _infusions: Vec<Infusion>,
-//     _cov: Option<HashMap<String, CovLine>>,
-// }
-// impl Model {
-//     pub fn get_param(&self, str: &str) -> f64 {
-//         *self.params.get(str).unwrap()
-//     }
-// }
-
-/// return the predicted values for the given scenario and parameters
-/// where the second element of the tuple is the predicted values
-/// one per observation time in scenario and in the same order
-/// it is not relevant the outeq of the specific event.
-// pub trait Predict<'a> {
-//     type Model: 'a + Clone;
-//     type State;
-//     fn initial_system(&self, params: &Vec<f64>, scenario: Scenario) -> (Self::Model, Scenario);
-//     fn initial_state(&self) -> Self::State;
-//     fn add_covs(&self, system: &mut Self::Model, cov: Option<HashMap<String, CovLine>>);
-//     fn add_infusion(&self, system: &mut Self::Model, infusion: Infusion);
-//     fn add_dose(&self, state: &mut Self::State, dose: f64, compartment: usize);
-//     fn get_output(&self, time: f64, state: &Self::State, system: &Self::Model, outeq: usize)
-//         -> f64;
-//     fn state_step(&self, state: &mut Self::State, system: &Self::Model, time: f64, next_time: f64);
-// }
-
-// #[derive(Clone, Debug, PartialEq)]
-// struct CacheKey {
-//     i: usize,
-//     support_point: Vec<f64>,
-// }
-
-// impl Eq for CacheKey {}
-
-// impl Hash for CacheKey {
-//     fn hash<H: Hasher>(&self, state: &mut H) {
-//         self.i.hash(state);
-//         for value in &self.support_point {
-//             value.to_bits().hash(state);
-//         }
-//     }
-// }
-
-// lazy_static! {
-//     static ref YPRED_CACHE: DashMap<CacheKey, Array1<f64>> =
-//         DashMap::with_capacity(CACHE_SIZE); // Adjust cache size as needed
-// }
-
-// fn get_ypred<S: Predict<'static> + Sync + Clone>(
-//     sim_eng: &Engine<S>,
-//     scenario: Scenario,
-//     support_point: Vec<f64>,
-//     i: usize,
-//     cache: bool,
-// ) -> Array1<f64> {
-//     let key = CacheKey {
-//         i,
-//         support_point: support_point.clone(),
-//     };
-//     if cache {
-//         match YPRED_CACHE.entry(key) {
-//             Entry::Occupied(entry) => entry.get().clone(), // Clone the cached value
-//             Entry::Vacant(entry) => {
-//                 let new_value = Array::from(sim_eng.pred(scenario, support_point));
-//                 entry.insert(new_value.clone());
-//                 new_value
-//             }
-//         }
-//     } else {
-//         Array::from(sim_eng.pred(scenario, support_point))
-//     }
-// }
-
-/// Simulate observations for multiple scenarios and support points.
-///
-/// This function performs simulation of observations for multiple scenarios and support points
-/// using the provided simulation engine `sim_eng`. It returns a 2D Array where each element
-/// represents the simulated observations for a specific scenario and support point.
-///
-/// # Arguments
-///
-/// * `sim_eng` - A reference to the simulation engine implementing the `Predict` trait.
-///
-/// * `scenarios` - A reference to a `Vec<Scenario>` containing information about different scenarios.
-///
-/// * `support_points` - A 2D Array `(Array2<f64>)` representing the support points. Each row
-///                     corresponds to a different support point scenario.
-///
-/// * `cache` - A boolean flag indicating whether to cache predicted values during simulation.
-///
-/// # Returns
-///
-/// A 2D Array `(Array2<Array1<f64>>)` where each element is an `Array1<f64>` representing the
-/// simulated observations for a specific scenario and support point.
-///
-/// # Example
-///
-///
-/// In this example, `observations` will contain the simulated observations for multiple scenarios
-/// and support points.
-///
-/// Note: This function allows for optional caching of predicted values, which can improve
-/// performance when simulating observations for multiple scenarios.
-///
-pub fn get_population_predictions<'a>(
-    equation: &'a Equation,
-    subjects: &Vec<Subject>,
-    support_points: &Array2<f64>,
-    _cache: bool,
-) -> PopulationPredictions {
-    let mut pred = Array2::default((subjects.len(), support_points.nrows()).f());
-    pred.axis_iter_mut(Axis(0))
-        .into_par_iter()
-        .enumerate()
-        .for_each(|(i, mut row)| {
-            row.axis_iter_mut(Axis(0))
-                .into_par_iter()
-                .enumerate()
-                .for_each(|(j, mut element)| {
-                    let subject = subjects.get(i).unwrap();
-                    let ypred =
-                        equation.simulate_subject(subject, support_points.row(j).to_vec().as_ref());
-                    element.fill(ypred);
-                });
-        });
-    pred.into()
-}
-
-// fn simple_sim<S>(sim_eng: &Engine<S>, scenario: Scenario, support_point: &Array1<f64>) -> Vec<f64>
-// where
-//     S: Predict<'static> + Sync + Clone,
-// {
-//     sim_eng.pred(scenario, support_point.to_vec())
-// }
-
-// fn post_predictions<S>(
-//     _equation: &Equation,
-//     _post: Array2<f64>,
-//     _subjects: &Vec<Subject>,
-// ) -> Result<Array1<Vec<f64>>, Box<dyn error::Error>> {
-//     unimplemented!();
-//     // if post.nrows() != scenarios.len() {
-//     //     return Err("Error calculating the posterior predictions, size mismatch.".into());
-//     // }
-//     // let mut predictions: Array1<Vec<f64>> = Array1::default(post.nrows());
-
-//     // predictions
-//     //     .axis_iter_mut(Axis(0))
-//     //     .into_par_iter()
-//     //     .enumerate()
-//     //     .for_each(|(i, mut pred)| {
-//     //         let scenario = scenarios.get(i).unwrap();
-//     //         let support_point = post.row(i).to_owned();
-//     //         pred.fill(simple_sim(sim_engine, scenario.clone(), &support_point))
-//     //     });
-
-//     // Ok(predictions)
-// }
diff --git a/src/simulator/analytical/mod.rs b/src/simulator/analytical/mod.rs
index 49eed400..2f163d65 100644
--- a/src/simulator/analytical/mod.rs
+++ b/src/simulator/analytical/mod.rs
@@ -1,4 +1,4 @@
-use crate::{routines::data::Covariates, simulator::*};
+use crate::{prelude::data::Covariates, simulator::*};
 
 // let eq = |x: &V, p: &V, t:T, rateiv: V, _cov: &Covariates|
 
diff --git a/src/simulator/likelihood.rs b/src/simulator/likelihood.rs
index 94c8a907..dbf2e6cf 100644
--- a/src/simulator/likelihood.rs
+++ b/src/simulator/likelihood.rs
@@ -1,6 +1,6 @@
-use crate::routines::{
-    data::Observation,
-    evaluation::sigma::{ErrorModel, ErrorType},
+use crate::{
+    prelude::data::Observation,
+    routines::evaluation::sigma::{ErrorModel, ErrorType},
 };
 
 use ndarray::{Array1, Array2, Axis};
diff --git a/src/simulator/mod.rs b/src/simulator/mod.rs
index a16cedf3..d73cfd2c 100644
--- a/src/simulator/mod.rs
+++ b/src/simulator/mod.rs
@@ -2,17 +2,22 @@ pub mod analytical;
 pub mod cache;
 pub mod likelihood;
 pub mod ode;
+use ndarray::{parallel::prelude::*, Axis};
 use std::collections::HashMap;
 
-use self::likelihood::SubjectPredictions;
+use self::likelihood::{PopulationPredictions, SubjectPredictions};
 use crate::{
-    prelude::data::Event,
-    routines::data::{Covariates, Infusion, OccasionTrait, Subject, SubjectTrait},
+    prelude::data::{Covariates, Event, Infusion, OccasionTrait, Subject, SubjectTrait},
     simulator::likelihood::ToPrediction,
 };
-// use diffsol::vector::Vector;
 
+// use dashmap::mapref::entry::Entry;
+// use dashmap::DashMap;
+// use lazy_static::lazy_static;
+use ndarray::prelude::*;
+// use ndarray::Array1;
 use cache::*;
+use ndarray::Array2;
 
 pub type T = f64;
 // pub type V = faer::Col<T>;
@@ -214,6 +219,30 @@ impl Equation {
     }
 }
 
+pub fn get_population_predictions<'a>(
+    equation: &'a Equation,
+    subjects: &Vec<Subject>,
+    support_points: &Array2<f64>,
+    _cache: bool,
+) -> PopulationPredictions {
+    let mut pred = Array2::default((subjects.len(), support_points.nrows()).f());
+    pred.axis_iter_mut(Axis(0))
+        .into_par_iter()
+        .enumerate()
+        .for_each(|(i, mut row)| {
+            row.axis_iter_mut(Axis(0))
+                .into_par_iter()
+                .enumerate()
+                .for_each(|(j, mut element)| {
+                    let subject = subjects.get(i).unwrap();
+                    let ypred =
+                        equation.simulate_subject(subject, support_points.row(j).to_vec().as_ref());
+                    element.fill(ypred);
+                });
+        });
+    pred.into()
+}
+
 trait FromVec {
     fn from_vec(vec: Vec<f64>) -> Self;
 }
diff --git a/src/simulator/ode/diffsol_traits.rs b/src/simulator/ode/diffsol_traits.rs
index 596d8ce1..e04aa511 100644
--- a/src/simulator/ode/diffsol_traits.rs
+++ b/src/simulator/ode/diffsol_traits.rs
@@ -1,11 +1,13 @@
-use crate::routines::data::{Covariates, Infusion};
+use crate::prelude::data::{Covariates, Infusion};
+use anyhow::Result;
 use diffsol::{
     matrix::Matrix,
     ode_solver::problem::OdeSolverProblem,
     op::{unit::UnitCallable, Op},
     vector::Vector,
-    OdeBuilder, OdeEquations, Result, Zero,
+    OdeBuilder, OdeEquations,
 };
+use num_traits::Zero;
 use std::rc::Rc;
 
 use super::closure::PMClosure;
diff --git a/src/simulator/ode/mod.rs b/src/simulator/ode/mod.rs
index 69acfb4b..000665cf 100644
--- a/src/simulator/ode/mod.rs
+++ b/src/simulator/ode/mod.rs
@@ -1,7 +1,7 @@
 pub mod closure;
 pub mod diffsol_traits;
 use crate::{
-    routines::data::{Covariates, Infusion},
+    prelude::data::{Covariates, Infusion},
     simulator::{ode::diffsol_traits::BuildPmOde, DiffEq, M, T, V},
 };