diff --git a/src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py b/src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py
index caf033d..6ee3735 100644
--- a/src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py
+++ b/src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py
@@ -1,4 +1,5 @@
 from src.wrappers.OsipiBase import OsipiBase
+import warnings
 import numpy as np
 import matlab.engine
 
@@ -47,7 +48,7 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
         self.use_initial_guess = {"f" : True, "D" : True, "Dp" : True, "S0" : True}
 
         if eng is None:
-            print('initiating matlab; this may take some time. For repeated testing one could use the optional input eng as an already initiated matlab engine')
+            warnings.warn('Initiating MATLAB; this may take some time. For repeated testing one could use the optional input eng as an already initiated MATLAB engine', UserWarning, stacklevel=2)
             self.eng=matlab.engine.start_matlab()
             self.keep_alive=False
         else:
@@ -100,7 +101,7 @@ def clean(self):
                 try:
                     self.eng.quit()
                 except Exception as e:
-                    print(f"Warning: Failed to quit MATLAB engine cleanly: {e}")
+                    warnings.warn(f"Failed to quit MATLAB engine cleanly: {e}", UserWarning, stacklevel=2)
 
     def __del__(self):
         self.clean()
diff --git a/src/standardized/ETP_SRI_LinearFitting.py b/src/standardized/ETP_SRI_LinearFitting.py
index ec33ecf..46d97c4 100644
--- a/src/standardized/ETP_SRI_LinearFitting.py
+++ b/src/standardized/ETP_SRI_LinearFitting.py
@@ -47,7 +47,7 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
 
         super(ETP_SRI_LinearFitting, self).__init__(bvalues, thresholds, bounds, initial_guess)
         if bounds is not None:
-            print('warning, bounds from wrapper are not (yet) used in this algorithm')
+            warnings.warn('Bounds from wrapper are not (yet) used in this algorithm', UserWarning, stacklevel=2)
         self.use_bounds = {"f": False, "Dp": False, "D": False, "S0": False}
         self.use_initial_guess = {"f": False, "Dp": False, "D": False, "S0": False}
 
diff --git a/src/standardized/IAR_LU_modified_topopro.py b/src/standardized/IAR_LU_modified_topopro.py
index caad456..58e6323 100644
--- a/src/standardized/IAR_LU_modified_topopro.py
+++ b/src/standardized/IAR_LU_modified_topopro.py
@@ -1,3 +1,4 @@
+import warnings
 import numpy as np
 from dipy.core.gradients import gradient_table
 from src.wrappers.OsipiBase import OsipiBase
@@ -53,7 +54,7 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
 
         # Check the inputs
         if self.bounds["Dp"][0] == self.bounds["D"][1]:
-            print('warning, bounds for D* and D are equal, this will likely cause fitting errors. Setting D_upper to 99 percent of D_upper')
+            warnings.warn('Bounds for D* and D are equal, this will likely cause fitting errors. Setting D_upper to 99 percent of D_upper', UserWarning, stacklevel=2)
             self.bounds["D"][1] = self.bounds["D"][1]*0.99
         # Check the inputs
         
diff --git a/src/standardized/IVIM_NEToptim.py b/src/standardized/IVIM_NEToptim.py
index c3cbdbe..e7c7719 100644
--- a/src/standardized/IVIM_NEToptim.py
+++ b/src/standardized/IVIM_NEToptim.py
@@ -66,7 +66,7 @@ def initialize(self, bounds, initial_guess, fitS0, traindata, SNR, n):
                 SNR = (5, 100)
             self.training_data(self.bvalues,n=n,SNR=SNR)
         self.arg=Arg()
-        print('note that the bounds in the network are soft bounds and implemented by a sigmoid transform. In order for the network to be sensitive over the range, we extend the bounds ny 30%')
+        warnings.warn('Note that the bounds in the network are soft bounds and implemented by a sigmoid transform. In order for the network to be sensitive over the range, we extend the bounds by 30%', UserWarning, stacklevel=2)
         if bounds is not None:
             self.bounds = bounds
         else:
diff --git a/src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py b/src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py
index e35d07f..a4e68e7 100644
--- a/src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py
+++ b/src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py
@@ -1,5 +1,6 @@
 from src.wrappers.OsipiBase import OsipiBase
 from src.original.OGC_AmsterdamUMC.LSQ_fitting import flat_neg_log_prior, fit_bayesian, empirical_neg_log_prior, fit_segmented, fit_bayesian_array, fit_segmented_array
+import warnings
 import numpy as np
 
 class OGC_AmsterdamUMC_Bayesian_biexp(OsipiBase):
@@ -62,10 +63,10 @@ def initialize(self, bounds=None, initial_guess=None, fitS0=True, prior_in=None,
         self.thresholds = thresholds
 
         if prior_in is None:
-            print('using a flat prior between bounds')
+            warnings.warn('Using a flat prior between bounds', UserWarning, stacklevel=2)
             self.neg_log_prior=flat_neg_log_prior([self.bounds["D"][0],self.bounds["D"][1]],[self.bounds["f"][0],self.bounds["f"][1]],[self.bounds["Dp"][0],self.bounds["Dp"][1]],[self.bounds["S0"][0],self.bounds["S0"][1]])
         else:
-            print('warning, bounds are not used, as a prior is used instead')
+            warnings.warn('Bounds are not used, as a prior is used instead', UserWarning, stacklevel=2)
             if len(prior_in) == 4:
                 self.neg_log_prior = empirical_neg_log_prior(prior_in[0], prior_in[1], prior_in[2],prior_in[3])
             else:
diff --git a/src/standardized/OGC_AmsterdamUMC_biexp_segmented.py b/src/standardized/OGC_AmsterdamUMC_biexp_segmented.py
index 5cfe777..c3ac076 100644
--- a/src/standardized/OGC_AmsterdamUMC_biexp_segmented.py
+++ b/src/standardized/OGC_AmsterdamUMC_biexp_segmented.py
@@ -1,5 +1,6 @@
 from src.wrappers.OsipiBase import OsipiBase
 from src.original.OGC_AmsterdamUMC.LSQ_fitting import fit_segmented, fit_segmented_array
+import warnings
 import numpy as np
 
 class OGC_AmsterdamUMC_biexp_segmented(OsipiBase):
@@ -53,7 +54,7 @@ def initialize(self, thresholds):
 
         if self.thresholds is None:
             self.thresholds = 150
-            print('warning, no thresholds were defined, so default threshold of 150 was used')
+            warnings.warn('No thresholds were defined, so default threshold of 150 was used', UserWarning, stacklevel=2)
         else:
             self.thresholds = thresholds
 
diff --git a/src/standardized/OJ_GU_bayesMATLAB.py b/src/standardized/OJ_GU_bayesMATLAB.py
index 889e843..72ccaa8 100644
--- a/src/standardized/OJ_GU_bayesMATLAB.py
+++ b/src/standardized/OJ_GU_bayesMATLAB.py
@@ -1,4 +1,5 @@
 from src.wrappers.OsipiBase import OsipiBase
+import warnings
 import numpy as np
 import matlab.engine
 
@@ -49,7 +50,7 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
         self.use_bounds = {"f" : True, "D" : True, "Dp" : True, "S0" : True}
 
         if eng is None:
-            print('initiating matlab; this may take some time. For repeated testing one could use the optional input eng as an already initiated matlab engine')
+            warnings.warn('Initiating MATLAB; this may take some time. For repeated testing one could use the optional input eng as an already initiated MATLAB engine', UserWarning, stacklevel=2)
             self.eng=matlab.engine.start_matlab()
             self.keep_alive=False
         else:
@@ -101,7 +102,7 @@ def clean(self):
                 try:
                     self.eng.quit()
                 except Exception as e:
-                    print(f"Warning: Failed to quit MATLAB engine cleanly: {e}")
+                    warnings.warn(f"Failed to quit MATLAB engine cleanly: {e}", UserWarning, stacklevel=2)
 
     def __del__(self):
         self.clean()
diff --git a/src/standardized/OJ_GU_segMATLAB.py b/src/standardized/OJ_GU_segMATLAB.py
index 37b2399..e2d6812 100644
--- a/src/standardized/OJ_GU_segMATLAB.py
+++ b/src/standardized/OJ_GU_segMATLAB.py
@@ -1,4 +1,5 @@
 from src.wrappers.OsipiBase import OsipiBase
+import warnings
 import numpy as np
 import matlab.engine
 
@@ -47,7 +48,7 @@ def __init__(self, bvalues=None, thresholds=200, bounds=None, initial_guess=None
         self.use_bounds = {"f" : True, "D" : True, "Dp" : True, "S0" : True}
         self.use_initial_guess = {"f" : False, "D" : False, "Dp" : False, "S0" : False}
         if eng is None:
-            print('initiating matlab; this may take some time. For repeated testing one could use the optional input eng as an already initiated matlab engine')
+            warnings.warn('Initiating MATLAB; this may take some time. For repeated testing one could use the optional input eng as an already initiated MATLAB engine', UserWarning, stacklevel=2)
             self.eng=matlab.engine.start_matlab()
             self.keep_alive=False
         else:
@@ -95,7 +96,7 @@ def clean(self):
                 try:
                     self.eng.quit()
                 except Exception as e:
-                    print(f"Warning: Failed to quit MATLAB engine cleanly: {e}")
+                    warnings.warn(f"Failed to quit MATLAB engine cleanly: {e}", UserWarning, stacklevel=2)
 
     def __del__(self):
         self.clean()
diff --git a/src/standardized/PvH_KB_NKI_IVIMfit.py b/src/standardized/PvH_KB_NKI_IVIMfit.py
index 20bf8f8..febeeaa 100644
--- a/src/standardized/PvH_KB_NKI_IVIMfit.py
+++ b/src/standardized/PvH_KB_NKI_IVIMfit.py
@@ -47,7 +47,7 @@ def __init__(self, bvalues=None, thresholds=None,bounds=None,initial_guess=None)
         super(PvH_KB_NKI_IVIMfit, self).__init__(bvalues=bvalues, thresholds=thresholds,bounds=bounds,initial_guess=initial_guess)
         self.NKI_algorithm = generate_IVIMmaps_standalone
         if bounds is not None:
-            print('warning, bounds from wrapper are not (yet) used in this algorithm')
+            warnings.warn('Bounds from wrapper are not (yet) used in this algorithm', UserWarning, stacklevel=2)
         self.use_bounds = {"f" : False, "D" : False, "Dp" : False, "S0" : False}
         self.use_initial_guess = {"f" : False, "D" : False, "Dp" : False, "S0" : False}
 
diff --git a/src/standardized/TCML_TechnionIIT_SLS.py b/src/standardized/TCML_TechnionIIT_SLS.py
index 1a07c44..991b2a5 100644
--- a/src/standardized/TCML_TechnionIIT_SLS.py
+++ b/src/standardized/TCML_TechnionIIT_SLS.py
@@ -55,7 +55,7 @@ def initialize(self, bounds, fitS0,thresholds):
         self.fitS0=fitS0
         if thresholds is None:
             self.thresholds = 150
-            print('warning, no thresholds were defined, so default bounds are used of  150')
+            warnings.warn('No thresholds were defined, so default threshold of 150 is used', UserWarning, stacklevel=2)
         else:
             self.thresholds = thresholds
 
diff --git a/src/standardized/TCML_TechnionIIT_lsq_sls_BOBYQA.py b/src/standardized/TCML_TechnionIIT_lsq_sls_BOBYQA.py
index a80bb2e..7c9c81e 100644
--- a/src/standardized/TCML_TechnionIIT_lsq_sls_BOBYQA.py
+++ b/src/standardized/TCML_TechnionIIT_lsq_sls_BOBYQA.py
@@ -1,5 +1,6 @@
 from src.wrappers.OsipiBase import OsipiBase
 from super_ivim_dc.source.Classsic_ivim_fit import IVIM_fit_sls_BOBYQA
+import warnings
 import numpy as np
 
 class TCML_TechnionIIT_lsq_sls_BOBYQA(OsipiBase):
@@ -50,7 +51,7 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
     def initialize(self, bounds, fitS0, thresholds):
         if self.thresholds is None:
             self.thresholds = 200
-            print('warning, no thresholds were defined, so default bounds are used of  200')
+            warnings.warn('No thresholds were defined, so default threshold of 200 is used', UserWarning, stacklevel=2)
         else:
             self.thresholds = thresholds
         self.fitS0=fitS0
diff --git a/src/standardized/TCML_TechnionIIT_lsq_sls_lm.py b/src/standardized/TCML_TechnionIIT_lsq_sls_lm.py
index e72a16a..065fb92 100644
--- a/src/standardized/TCML_TechnionIIT_lsq_sls_lm.py
+++ b/src/standardized/TCML_TechnionIIT_lsq_sls_lm.py
@@ -54,7 +54,7 @@ def initialize(self, bounds, fitS0,thresholds):
 
         if thresholds is None:
             self.thresholds = 150
-            print('warning, no thresholds were defined, so default bounds are used of  150')
+            warnings.warn('No thresholds were defined, so default threshold of 150 is used', UserWarning, stacklevel=2)
         else:
             self.thresholds = thresholds
         self.fitS0=fitS0
diff --git a/src/standardized/TCML_TechnionIIT_lsq_sls_trf.py b/src/standardized/TCML_TechnionIIT_lsq_sls_trf.py
index bd37628..eeecd66 100644
--- a/src/standardized/TCML_TechnionIIT_lsq_sls_trf.py
+++ b/src/standardized/TCML_TechnionIIT_lsq_sls_trf.py
@@ -1,5 +1,6 @@
 from src.wrappers.OsipiBase import OsipiBase
 from super_ivim_dc.source.Classsic_ivim_fit import IVIM_fit_sls_trf
+import warnings
 import numpy as np
 
 class TCML_TechnionIIT_lsq_sls_trf(OsipiBase):
@@ -56,7 +57,7 @@ def initialize(self, bounds, fitS0, thresholds):
 
         if thresholds is None:
             self.thresholds = 200
-            print('warning, no thresholds were defined, so default bounds are used of  200')
+            warnings.warn('No thresholds were defined, so default threshold of 200 is used', UserWarning, stacklevel=2)
         else:
             self.thresholds = thresholds
         self.fitS0=fitS0
diff --git a/src/standardized/TF_reference_IVIMfit.py b/src/standardized/TF_reference_IVIMfit.py
index eba1630..eee58f3 100644
--- a/src/standardized/TF_reference_IVIMfit.py
+++ b/src/standardized/TF_reference_IVIMfit.py
@@ -1,5 +1,6 @@
 from src.wrappers.OsipiBase import OsipiBase
 from src.original.TF_reference.segmented_IVIMfit import segmented_IVIM_fit
+import warnings
 import numpy as np
 
 class TF_reference_IVIMfit(OsipiBase):
@@ -51,7 +52,7 @@ def initialize(self, bounds, thresholds):
         self.use_bounds = {"f" : True, "D" : True, "Dp" : True, "S0" : True}
         if self.thresholds is None:
             self.thresholds = 200
-            print('warning, no thresholds were defined, so default threshold are used of 200')
+            warnings.warn('No thresholds were defined, so default threshold of 200 is used', UserWarning, stacklevel=2)
         else:
             self.thresholds = thresholds
 
@@ -73,7 +74,7 @@ def ivim_fit(self, signals, bvalues=None):
 
         if np.any(signals < 0):
             signals = np.clip(signals,0.01, None)
-            print('warning, negative values in signal: values clipped to 0.01')
+            warnings.warn('Negative values in signal: values clipped to 0.01', UserWarning, stacklevel=2)
 
 
         fit_results = self.TF_reference_algorithm(bvalues,signals,b_cutoff=self.thresholds, bounds=bounds)
diff --git a/src/wrappers/OsipiBase.py b/src/wrappers/OsipiBase.py
index 1e5e166..3555dde 100644
--- a/src/wrappers/OsipiBase.py
+++ b/src/wrappers/OsipiBase.py
@@ -3,6 +3,7 @@
 from scipy.stats import norm
 import pathlib
 import sys
+import warnings
 from tqdm import tqdm
 from joblib import Parallel, delayed
 
@@ -149,12 +150,12 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             else:
                 # Generic defaults (original behavior)
                 if self.bounds is None:
-                    print('warning, no bounds were defined, so default bounds are used of [0, 0, 0.005, 0.7],[0.005, 1.0, 0.2, 1.3]')
+                    warnings.warn('No bounds were defined, so default bounds are used of [0, 0, 0.005, 0.7],[0.005, 1.0, 0.2, 1.3]', UserWarning, stacklevel=2)
                     self.bounds = {"S0" : [0.7, 1.3], "f" : [0, 1.0], "Dp" : [0.005, 0.2], "D" : [0, 0.005]} # These are defined as [lower, upper]
                     self.forced_default_bounds = True
 
                 if self.initial_guess is None:
-                    print('warning, no initial guesses were defined, so default initial guesses are used of  [0.001, 0.001, 0.01, 1]')
+                    warnings.warn('No initial guesses were defined, so default initial guesses are used of [0.001, 0.001, 0.01, 1]', UserWarning, stacklevel=2)
                     self.initial_guess = {"S0" : 1, "f" : 0.1, "Dp" : 0.01, "D" : 0.001}
                     self.forced_default_initial_guess = True