diff --git a/configs/.small_network_tf.yml.swp b/configs/.small_network_tf.yml.swp deleted file mode 100644 index d2f1af4..0000000 Binary files a/configs/.small_network_tf.yml.swp and /dev/null differ diff --git a/configs/automated_tf.yml b/configs/automated_tf.yml new file mode 100644 index 0000000..ba587c6 --- /dev/null +++ b/configs/automated_tf.yml @@ -0,0 +1,62 @@ +FileSettings: + # FileName : "/home/henryi/sft/MaCh3Tutorial/Test.root" + # FileName: "/home/henryi/scratch/Delayed_Fit_Long_Chains_no_delay_adapt_0_4.root" + SkipFileLoading: True + FileName: "/home/henryi/scratch/thinned_partial_chain.root" + ChainName : "posteriors" + Verbose : False + MakeMLModel: True + RunLLHScan: False + RunMCMC: True + + +ParameterSettings: + ParameterNames : ["sin2th","delm2", "delta", "xsec", "sk", "nd"] + LabelName : "LogL" + IgnoredParameters : ["LogL_systematic_xsec_cov", "Log", + "LogL_systematic_nddet_cov"] + ParameterCuts : ["LogL<5000", "LogL>4900", "delm2_23>0"] + CircularParameter: ["delta_cp"] + + +MLSettings: + FitterPackage : "TensorFlow" + FitterName : "autotune" + TestSize : 0.2 + + MLOutputFile: "tf_model_normal_order_small.keras" + MLScalerOutputName: "scaler_tf_model_normal_order_small.pkl" + + PlotOutputName : "tf_normal_order.pdf" + + AddFromExternalModel: False + ExternalModel: "tf_model_normal_order.keras" + ExternalScaler: "scaler_tf_model_normal_order_small.pkl" + + FitterKwargs: + BuildSettings: + epochs: 10000 + batch_size: 4096 + # Tuning features + n_trials: 30 + + hyperband_iterations: 10 + tuning_dir: "~/scratch/t2k_full_tuning" + project_name: "t2k_linear_tune" + + n_layers: [5, 25, 2] + activation: ['tanh', 'relu', 'swish', 'selu', 'elu', 'leak_relu', 'softplus'] + neurons_per_layer: [24, 1024, 100] + learning_rate: [0.001, 0.0001, 0.00001] + regularization: [0.0001, 0.1, 1.0] + +MCMCSettings: + NSteps: 1000000 + NChains: 20 + UpdateStep: 100 + MaxUpdateSteps: 500000 + + MCMCOutput: "/home/henryi/scratch/mcmc_chain_GPU.h5" + +LikelihoodScanSettings: + NDivisions: 10000 \ No newline at end of file diff --git a/configs/big_network.yml b/configs/big_network.yml deleted file mode 100644 index 37ffa36..0000000 --- a/configs/big_network.yml +++ /dev/null @@ -1,103 +0,0 @@ -FileSettings: - # FileName : "/home/henryi/sft/MaCh3Tutorial/Test.root" - # FileName: "/home/henryi/sft/MaCh3Tutorial/Test.root" - # FileName: "/home/henryi/scratch/Delayed_Fit_Long_Chains_no_delay_adapt_0_4.root" - - FileName: "/home/henryi/scratch/all_par_chain_without_adaption.root" - ChainName : "posteriors" - - Verbose : False - MakeMLModel: True - RunLLHScan: False - RunMCMC: False - -ParameterSettings: - ParameterNames : ["sin2th","delm2", "delta", "xsec", "sk", "nd"] - LabelName : "LogL" - IgnoredParameters : ["LogL_systematic_xsec_cov", "Log", - "LogL_systematic_nddet_cov"] - ParameterCuts : ["LogL<12345678", "delm2_23>0"] - CircularParameter: ["delta_cp"] - - -MLSettings: - FitterPackage : "TensorFlow" - FitterName : "Sequential" - TestSize : 0.2 - MLOutputFile: "tf_model_normal_order.keras" - PlotOutputName : "tf_tutorial_chain_no.pdf" - - FitterKwargs: - BuildSettings: - optimizer: 'adam' - loss: 'mse' - metrics: ['mae', 'mse'] - - Layers: - - dense: - units: 64 - activation: 'relu' - - dense: - units: 64 - activation: 'relu' - - dense: - units: 64 - activation: 'relu' - - dense: - units: 64 - activation: 'relu' - - dense: - units: 64 - activation: 'relu' - - dense: - units: 64 - activation: 'relu' - - dense: - units: 64 - activation: 'relu' - - dense: - units: 64 - activation: 'relu' - - dense: - units: 64 - activation: 'relu' - - dense: - units: 64 - activation: 'relu' - - dense: - units: 64 - activation: 'relu' - - dense: - units: 64 - activation: 'relu' - - dense: - units: 64 - activation: 'relu' - - dense: - units: 32 - activation: 'relu' - - dense: - units: 16 - activation: 'relu' - - dense: - units: 8 - activation: 'relu' - - - dense: - units: 1 - activation: 'linear' - - FitSettings: - batch_size: 2048 - epochs: 512 - -MCMCSettings: - NSteps: 1000000 - NChains: 1 - UpdateStep: 1000 - BatchSize: 10 - - MCMCOutput: "/home/henryi/scratch/tutorial_chain_test.h5" - -LikelihoodScanSettings: - NDivisions: 10000 \ No newline at end of file diff --git a/configs/normalizing_flow.yml b/configs/normalizing_flow.yml new file mode 100644 index 0000000..d92ec63 --- /dev/null +++ b/configs/normalizing_flow.yml @@ -0,0 +1,58 @@ +FileSettings: + # FileName : "/home/henryi/sft/MaCh3Tutorial/Test.root" + # FileName: "/home/henryi/scratch/Delayed_Fit_Long_Chains_no_delay_adapt_0_4.root" + SkipFileLoading: True + FileName: "/home/henryi/scratch/thinned_partial_chain.root" + ChainName : "posteriors" + Verbose : False + MakeMLModel: True + RunLLHScan: False + RunMCMC: False + + +ParameterSettings: + ParameterNames : ["sin2th","delm2", "delta", "xsec", "sk", "nd"] + LabelName : "LogL" + IgnoredParameters : ["LogL_systematic_xsec_cov", "Log", + "LogL_systematic_nddet_cov"] + ParameterCuts : ["LogL<5000", "LogL>4900", "delm2_23>0"] + CircularParameter: ["delta_cp"] + + +MLSettings: + FitterPackage : "TensorFlow" + FitterName : "normalizing_flow" + TestSize : 0.2 + + MLOutputFile: "tf_model_normal_order_small.keras" + MLScalerOutputName: "scaler_tf_model_normal_order_small.pkl" + + PlotOutputName : "tf_normal_order.pdf" + + AddFromExternalModel: False + ExternalModel: "tf_model_normal_order.keras" + ExternalScaler: "scaler_tf_model_normal_order_small.pkl" + + FitterKwargs: + BuildSettings: + learning_rate: 0.01 + hidden_units: [2000, 2000, 2000, 2000, 2000, 2000] + loc: 0 + scale: 1 + + FitSettings: + batch_size: 4096 + epochs: 250 + validation_split: 0.1 + shuffle: True + steps_per_epoch: 1 +MCMCSettings: + NSteps: 1000000 + NChains: 20 + UpdateStep: 100 + MaxUpdateSteps: 500000 + + MCMCOutput: "/home/henryi/scratch/mcmc_chain_GPU.h5" + +LikelihoodScanSettings: + NDivisions: 10000 \ No newline at end of file diff --git a/configs/normalizing_flow_config.yml b/configs/normalizing_flow_config.yml deleted file mode 100644 index 5f3b1a4..0000000 --- a/configs/normalizing_flow_config.yml +++ /dev/null @@ -1,23 +0,0 @@ -FileSettings: - FileName : "/home/henryi/sft/MaCh3Tutorial/Test.root" - ChainName : "posteriors" - Verbose : True - MakeMLModel: True - -ParameterSettings: - ParameterNames : ["sin2th","delm2", "delta", "xsec"] - LabelName : "LogL" - IgnoredParameters : ["LogL_systematic_xsec_cov", "LogL_systematic_osc_cov"] - ParameterCuts : [LogL<12345678] - -MLSettings: - FitterPackage : "NormalizingFlow" - FitterName : "realnvp" - TestSize : 0.6 - - FitterKwargs: - n_iter: 1000 - layer_structure: [8, 64, 64, 1] - n_layers: 32 - - MLOutputName: "flow_model_full.pkl" \ No newline at end of file diff --git a/configs/plotting_config.yml b/configs/plotting_config.yml deleted file mode 100644 index 22bd673..0000000 --- a/configs/plotting_config.yml +++ /dev/null @@ -1,50 +0,0 @@ -FileSettings: - FileName : "/users/php20hti/php20hti/public/chains/all_par_chain_without_adaption.root" - # FileName : "/users/php20hti/t2ksft/regeneration_times_project/inputs/markov_chains/oscpar_only_datafit.root" - ChainName : "posteriors" - # ChainName : "osc_posteriors" - Verbose : False - # Model Settings - MakeDiagnostics: True - MakePosteriors: False - -ParameterSettings: - ParameterNames : ["xsec"] - # ParameterNames : ["d", "th"] - IgnoredParameters : ["LogL_systematic_xsec_cov", "LogL_systematic_nddet_cov", ] - ParameterCuts : ["step>10000"] - CircularParameters: [] - -PlottingSettings: - DiagnosticsSettings: - # Where are we plotting? - DiagnosticsOutputFile: "diagnostics_output.pdf" - # Make Trace/AC Plot? - MakeTraceAC: True - # Make Violin Plot? - MakeViolin: False - #Make ESS Plot? - MakeESS: False - # Make MCSE Plot - MakeMCSE: False - # Make suboptimality Plot - MakeSuboptimality: False - # Steps/calculation for subopt. - SuboptimalitySteps: 10000 - # Print summary stats - PrintSummary: True - - PosteriorSettings: - # Output PDF - PosteriorOutputFile: "posteriors.pdf" - # Do you want 1D CIs? - Make1DPosteriors: True - # Plotted CIs - CredibleIntervals: [0.6, 0.90, 0.95] - # Do you want 2D CIs? - Make2DPosteriors: False - # Do you want a triangle plot? - MakeTrianglePlot: False - # variables to put in the triangle plot - TrianglePlot: [] - \ No newline at end of file diff --git a/configs/small_network_tf.yml b/configs/small_network_tf.yml index 778846e..94c9344 100644 --- a/configs/small_network_tf.yml +++ b/configs/small_network_tf.yml @@ -7,7 +7,7 @@ FileSettings: Verbose : False MakeMLModel: True RunLLHScan: False - RunMCMC: True + RunMCMC: False ParameterSettings: @@ -15,13 +15,13 @@ ParameterSettings: LabelName : "LogL" IgnoredParameters : ["LogL_systematic_xsec_cov", "Log", "LogL_systematic_nddet_cov"] - ParameterCuts : ["LogL<12345678", "delm2_23>0"] + ParameterCuts : ["LogL<5000", "LogL>4900", "delm2_23>0"] CircularParameter: ["delta_cp"] MLSettings: FitterPackage : "TensorFlow" - FitterName : "Sequential" + FitterName : "sequential" TestSize : 0.2 MLOutputFile: "tf_model_normal_order_small.keras" @@ -29,100 +29,77 @@ MLSettings: PlotOutputName : "tf_normal_order.pdf" - AddFromExternalModel: True - ExternalModel: "tf_model_normal_order.keras" + AddFromExternalModel: False + ExternalModel: "tf_model_normal_order.keras" ExternalScaler: "scaler_tf_model_normal_order_small.pkl" FitterKwargs: BuildSettings: - optimizer: 'adam' loss: 'mse' metrics: ['mae', 'mse'] - - + learning_rate: 0.0001 Layers: + # Input worked out automatically + # 0 + # 1 - dense: - units: 256 - activation: 'relu' - - dense: - units: 256 - activation: 'relu' - - dense: - units: 256 - activation: 'relu' - - dense: - units: 256 - activation: 'relu' - - dense: - units: 256 - activation: 'relu' - - dense: - units: 256 - activation: 'relu' + units: 1024 + activation: 'PReLU' + # 2 - dense: - units: 256 - activation: 'relu' - - dense: - units: 256 - activation: 'relu' - - dense: - units: 256 - activation: 'relu' - - dense: - units: 256 - activation: 'relu' - - dense: - units: 256 - activation: 'relu' + units: 1024 + activation: 'PReLU' + kernel_regularizer: 0.01 + - dense: - units: 256 - activation: 'relu' + units: 1024 + activation: 'PReLU' + - batchnorm: + momentum: 0.99 + # 3 + - dense: - units: 256 - activation: 'relu' + units: 1024 + kernel_regularizer: 0.01 + + activation: 'PReLU' + - batchnorm: + momentum: 0.99 - dense: - units: 256 - activation: 'relu' + units: 1024 + activation: 'PReLU' + - batchnorm: + momentum: 0.99 - dense: - units: 256 - activation: 'relu' + units: 1024 + activation: 'PReLU' + - batchnorm: + momentum: 0.99 - dense: - units: 256 - activation: 'relu' + units: 1024 + activation: 'PReLU' + - batchnorm: + momentum: 0.99 - dense: - units: 256 - activation: 'relu' - - dense: - units: 256 - activation: 'relu' + units: 1024 + activation: 'PReLU' + - batchnorm: + momentum: 0.99 - dense: - units: 256 - activation: 'relu' - - dense: - units: 256 - activation: 'relu' - - - # - dense: - # units: 32 - # activation: 'relu' - # - dense: - # units: 16 - # activation: 'relu' - # - dense: - # units: 8 - # activation: 'relu' + units: 1024 + activation: 'PReLU' + # output - dense: units: 1 - activation: 'linear' + activation: 'relu' FitSettings: - batch_size: 10000 - epochs: 10 + batch_size: 4096 + epochs: 250 validation_split: 0.1 MCMCSettings: diff --git a/configs/tensorflow_config.yml b/configs/tensorflow_config.yml index 42046e7..c14520f 100644 --- a/configs/tensorflow_config.yml +++ b/configs/tensorflow_config.yml @@ -1,10 +1,13 @@ FileSettings: - FileName : "/home/henryi/sft/MaCh3Tutorial/Test.root" + # FileName : "/home/henryi/sft/MaCh3Tutorial/Test.root" + # FileName: "/home/henryi/scratch/Delayed_Fit_Long_Chains_no_delay_adapt_0_4.root" + SkipFileLoading: True + FileName: "/home/henryi/scratch/thinned_partial_chain.root" ChainName : "posteriors" - Verbose : True + Verbose : False MakeMLModel: True RunLLHScan: False - RunMCMC: True + RunMCMC: False ParameterSettings: @@ -12,88 +15,93 @@ ParameterSettings: LabelName : "LogL" IgnoredParameters : ["LogL_systematic_xsec_cov", "Log", "LogL_systematic_nddet_cov"] - ParameterCuts : ["LogL<12345678"] + ParameterCuts : ["LogL<5000", "LogL>4900", "delm2_23>0"] + CircularParameter: ["delta_cp"] MLSettings: FitterPackage : "TensorFlow" - FitterName : "Sequential" - TestSize : 0.9 - MLOutputFile: "tf_disable_dropout.keras" - PlotOutputName : "tf_disable_dropout.pdf" + FitterName : "sequential" + TestSize : 0.2 + + MLOutputFile: "tf_model_normal_order_small.keras" + MLScalerOutputName: "scaler_tf_model_normal_order_small.pkl" + + PlotOutputName : "tf_normal_order.pdf" + + AddFromExternalModel: False + ExternalModel: "tf_model_normal_order.keras" + ExternalScaler: "scaler_tf_model_normal_order_small.pkl" FitterKwargs: BuildSettings: - optimizer: 'adam' loss: 'mse' metrics: ['mae', 'mse'] + learning_rate: 0.0001 Layers: + - batchnorm: + momentum: 0.99 - dense: - units: 75 - activation: 'relu' - - dense: - units: 75 - activation: 'relu' - - dense: - units: 75 - activation: 'relu' - - dense: - units: 75 - activation: 'relu' - - dense: - units: 75 - activation: 'relu' + units: 1024 + activation: 'swish' + + # 2 - dense: - units: 75 - activation: 'relu' + units: 1024 + activation: 'swish' + kernel_regularizer: 0.01 + - dense: - units: 75 - activation: 'relu' + units: 1024 + activation: 'swish' + - dense: - units: 75 - activation: 'relu' + units: 1024 + kernel_regularizer: 0.01 + kernel_regularizer: True + activation: 'swish' - dense: - units: 75 - activation: 'relu' + units: 1024 + kernel_regularizer: True + activation: 'swish' + - batchnorm: + momentum: 0.99 - dense: - units: 75 - activation: 'relu' + units: 1024 + activation: 'swish' + kernel_regularizer: True + - dense: - units: 50 - activation: 'relu' + units: 1024 + activation: 'swish' - dense: - units: 25 - activation: 'relu' + units: 1024 + activation: 'swish' + - batchnorm: + momentum: 0.99 + - dense: - units: 8 - activation: 'relu' + units: 1024 + activation: 'elu' + + # output - dense: units: 1 activation: 'linear' FitSettings: - batch_size: 2048 - epochs: 8 + batch_size: 4096 + epochs: 250 + validation_split: 0.1 MCMCSettings: - NSteps: 100 - NWalkers: 400 + NSteps: 1000000 + NChains: 20 + UpdateStep: 100 + MaxUpdateSteps: 500000 + MCMCOutput: "/home/henryi/scratch/mcmc_chain_GPU.h5" LikelihoodScanSettings: - xsec_0: [1, 0.7, 1.3] # Norm 0 - xsec_1: [1, 0.5, 1.5] # Norm - xsec_2: [1, 0, 2] # Norm 2 - xsec_3: [1.21, 0.0, 1.4] #spline 0 - xsec_4: [1, 0, 3] # Spline 1 - xsec_5: [0, -10, 15] # Func 0x - xsec_6: [0, -10, 15] # Func 1 - xsec_7: [0, -10, 15] # Func 2 - xsec_8: [0, -10, 15] # Func 3 - sin2th_12: [0.307, 0.25, 0.35] - sin2th_23: [0.546, 0.4, 0.7] - sin2th_13: [0.0220, 0.0, 0.12] - delm2_12: [0.0000753, 0.00007, 0.00008] - delm2_23: [0, -0.0026, 0.0026] - delta_cp: [0, -3.14159, 3.14159] \ No newline at end of file + NDivisions: 10000 \ No newline at end of file diff --git a/src/MaCh3PythonUtils/__init__.py b/src/MaCh3PythonUtils/__init__.py index c7e0ba0..08b2af4 100644 --- a/src/MaCh3PythonUtils/__init__.py +++ b/src/MaCh3PythonUtils/__init__.py @@ -1,2 +1,2 @@ -__version__="1.2.1" +__version__="2.0.0" __author__="Henry Wallace" diff --git a/src/MaCh3PythonUtils/config_reader/config_reader.py b/src/MaCh3PythonUtils/config_reader/config_reader.py index b6f457b..63bcf94 100644 --- a/src/MaCh3PythonUtils/config_reader/config_reader.py +++ b/src/MaCh3PythonUtils/config_reader/config_reader.py @@ -32,12 +32,6 @@ class ConfigReader: "ChainName": "", # More printouts "Verbose": False, - # Make posteriors from chain? - "MakePosteriors": False, - # Run Diagnostics code? - "MakeDiagnostics": False, - # Make an ML model to replicate the chain likelihood model - "MakeMLModel": False, # Run an LLH Scan? "RunLLHScan": False, # Run MCMC? @@ -56,42 +50,6 @@ class ConfigReader: "IgnoredParameters":[] }, - # Settings for plotting tools - "PlottingSettings":{ - # Specific Settings for posterior plots - "PosteriorSettings": { - # Make 2D Posterior Plots? - "Make2DPosteriors": False, - # Make a triangle plot - "MakeTrianglePlot": False, - # Variables in the triangle plot - "TrianglePlot": [], - # 1D credible intervals - "MakeCredibleIntervals": False, - # Output file - "PosteriorOutputFile": "posteriors.pdf" - }, - - # Specific Settings for diagnostic plots - "DiagnosticsSettings": { - # Make violin plot? - "MakeViolin": False, - # Make trace + AC plot? - "MakeTraceAC": False, - # Make effective sample size plot? - "MakeESS": False, - # Make MCSE plot? - "MakeMCSE": False, - # Make suboptimality plot - "MakeSuboptimality": False, - # Step for calculation - "SuboptimalitySteps": 0, - # Output file - "DiagnosticsOutputFile": "diagnostics.pdf", - # Print summary statistic - "PrintSummary": False - } - }, # Specific Settings for ML Applications "MLSettings": { # Name of plots @@ -157,69 +115,6 @@ def make_file_handler(self)->None: self._file_handler.convert_ttree_to_array() - - def make_posterior_plots(self): - """Generates posterior plots - """ - if self._plot_interface is None: - self._plot_interface = PlottingInterface(self._file_handler) - - posterior_labels = [] - - if self.__chain_settings['PlottingSettings']['PosteriorSettings']['Make2DPosteriors']: - self._plot_interface.initialise_new_plotter(PosteriorPlotter2D(self._file_handler), 'posterior_2d') - posterior_labels.append('posterior_2d') - - if self.__chain_settings['PlottingSettings']['PosteriorSettings']['MakeTrianglePlot']: - self._plot_interface.initialise_new_plotter(TrianglePlotter(self._file_handler), 'posterior_triangle') - posterior_labels.append('posterior_triangle') - - # Which variables do we actually want 2D plots for? - self._plot_interface.set_variables_to_plot(self.__chain_settings['PlottingSettings']['PosteriorSettings']['TrianglePlot'], posterior_labels) - - if self.__chain_settings['PlottingSettings']['PosteriorSettings']['Make1DPosteriors']: - self._plot_interface.initialise_new_plotter(PosteriorPlotter1D(self._file_handler), 'posterior_1d') - posterior_labels.append('posterior_1d') - - - self._plot_interface.set_credible_intervals(self.__chain_settings['PlottingSettings']['PosteriorSettings']['CredibleIntervals']) - self._plot_interface.set_is_circular(self.__chain_settings['ParameterSettings']['CircularParameters']) - - self._plot_interface.make_plots(self.__chain_settings['PlottingSettings']['PosteriorSettings']['PosteriorOutputFile'], posterior_labels) - - def make_diagnostics_plots(self): - """Generates diagnostics plots - """ - if self._plot_interface is None: - self._plot_interface = PlottingInterface(self._file_handler) - - diagnostic_labels = [] - - if self.__chain_settings['PlottingSettings']['DiagnosticsSettings']['MakeViolin']: - self._plot_interface.initialise_new_plotter(ViolinPlotter(self._file_handler), 'violin_plotter') - diagnostic_labels.append('violin_plotter') - if self.__chain_settings['PlottingSettings']['DiagnosticsSettings']['MakeTraceAC']: - self._plot_interface.initialise_new_plotter(AutocorrelationTracePlotter(self._file_handler), 'trace_autocorr') - diagnostic_labels.append('trace_autocorr') - if self.__chain_settings['PlottingSettings']['DiagnosticsSettings']['MakeESS']: - self._plot_interface.initialise_new_plotter(EffectiveSampleSizePlotter(self._file_handler), 'ess_plot') - diagnostic_labels.append('ess_plot') - if self.__chain_settings['PlottingSettings']['DiagnosticsSettings']['MakeMCSE']: - self._plot_interface.initialise_new_plotter(MarkovChainStandardError(self._file_handler), 'msce_plot') - diagnostic_labels.append('msce_plot') - - self._plot_interface.make_plots(self.__chain_settings['PlottingSettings']['DiagnosticsSettings']['DiagnosticsOutputFile'], diagnostic_labels) - - # Final one, covariance plotter - if self.__chain_settings['PlottingSettings']['DiagnosticsSettings']['MakeSuboptimality']: - suboptimality_obj = CovarianceMatrixUtils(self._file_handler) - suboptimality_obj.calculate_suboptimality(self.__chain_settings['PlottingSettings']['DiagnosticsSettings']['SuboptimalitySteps'], self.__chain_settings['PlottingSettings']['DiagnosticsSettings']['SubOptimalityMin']) - suboptimality_obj.plot_suboptimality(f"suboptimality_{self.__chain_settings['PlottingSettings']['DiagnosticsSettings']['DiagnosticsOutputFile']}") - - # Finally let's make a quick simmary - if self.__chain_settings['PlottingSettings']['DiagnosticsSettings']['PrintSummary']: - self._plot_interface.print_summary(f"summary_{self.__chain_settings['PlottingSettings']['DiagnosticsSettings']['DiagnosticsOutputFile']}.txt") - def make_ml_interface(self)->None: """Generates ML interface objects """ @@ -245,36 +140,29 @@ def make_ml_interface(self)->None: self._interface.save_model(self.__chain_settings["MLSettings"]["MLOutputFile"]) self._interface.save_scaler(self.__chain_settings['MLSettings']['MLScalerOutputName']) - def __call__(self) -> None: - """Runs over all files from config - """ - - if self.__chain_settings["FileSettings"]["SkipFileLoading"]: - self.make_file_handler() - - if self.__chain_settings["FileSettings"]["MakePosteriors"]: - self.make_posterior_plots() - - if self.__chain_settings["FileSettings"]["MakeDiagnostics"]: - self.make_diagnostics_plots() - - if self.__chain_settings["FileSettings"]["MakeMLModel"]: - self.make_ml_interface() - - if self.__chain_settings["FileSettings"]["RunLLHScan"] and self._interface is not None: - self._interface.run_likelihood_scan(self.__chain_settings["LikelihoodScanSettings"]["NDivisions"]) - - if self.__chain_settings["FileSettings"]["RunMCMC"] and self._interface is not None: - print("WARNING: MCMC HAS ONLY BEEN TESTED WITH TENSORFLOW INTERFACES!") + def run_mcmc(self): + print("WARNING: MCMC HAS ONLY BEEN TESTED WITH TENSORFLOW INTERFACES!") + mcmc = MCMCMultGPU(self._interface, + self.__chain_settings["MCMCSettings"]["NChains"], + self.__chain_settings["ParameterSettings"]["CircularParameters"], + self.__chain_settings["MCMCSettings"]["UpdateStep"], + self.__chain_settings["MCMCSettings"]["MaxUpdateSteps"]) - mcmc = MCMCMultGPU(self._interface, - self.__chain_settings["MCMCSettings"]["NChains"], - self.__chain_settings["ParameterSettings"]["CircularParameters"], - self.__chain_settings["MCMCSettings"]["UpdateStep"], - self.__chain_settings["MCMCSettings"]["MaxUpdateSteps"]) + mcmc(self.__chain_settings["MCMCSettings"]["NSteps"], + self.__chain_settings["MCMCSettings"]["MCMCOutput"],) - mcmc(self.__chain_settings["MCMCSettings"]["NSteps"], - self.__chain_settings["MCMCSettings"]["MCMCOutput"],) + + def __call__(self) -> None: + """Runs over all files from config + """ + self.make_file_handler() + + self.make_ml_interface() + if self.__chain_settings["FileSettings"]["RunLLHScan"] and self._interface is not None: + self._interface.run_likelihood_scan(self.__chain_settings["LikelihoodScanSettings"]["NDivisions"]) + + if self.__chain_settings["FileSettings"]["RunMCMC"] and self._interface is not None: + self.run_mcmc() diff --git a/src/MaCh3PythonUtils/diagnostics/interface/__init__.py b/src/MaCh3PythonUtils/diagnostics/interface/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/src/MaCh3PythonUtils/diagnostics/interface/plotting_interface.py b/src/MaCh3PythonUtils/diagnostics/interface/plotting_interface.py deleted file mode 100644 index 929966b..0000000 --- a/src/MaCh3PythonUtils/diagnostics/interface/plotting_interface.py +++ /dev/null @@ -1,152 +0,0 @@ -''' -Interface class for making plots! -''' -from typing import List -from MaCh3PythonUtils.file_handling.chain_handler import ChainHandler -from MaCh3PythonUtils.diagnostics.mcmc_plots.plotter_base import _PlottingBaseClass -from MaCh3PythonUtils.diagnostics.mcmc_plots.posteriors.posterior_base_classes import _PosteriorPlottingBase -from matplotlib.backends.backend_pdf import PdfPages -import arviz as az - -class PlottingInterface: - def __init__(self, chain_handler: ChainHandler)->None: - """Interface for handling plotting objects - - :param chain_handler: ChainHandler instance - :type chain_handler: ChainHandler - """ - self._chain_handler = chain_handler - self._plotter_object_dict = {} # dict of objects from plotting tools - - self._chain_handler.make_arviz_tree() - - def initialise_new_plotter(self, new_plotter: _PlottingBaseClass , plot_label: str)->None: - """Initialises a new plotter object - - :param new_plotter: Instance of new plotting class - :type new_plotter: pt.plotter_base._PlottingBaseClass - :param plot_label: Label for calling this instance - :type plot_label: str - """ - self._plotter_object_dict[plot_label] = new_plotter - - def set_credible_intervals(self, credible_intervals: List[float])->None: - """Sets credible intervals for for posterior plots - - :param credible_intervals: List of credivle intervals - :type credible_intervals: List[float] - :raises ValueError: Checks if crredible intervals are a valid type - """ - - print(f"Setting credible intervals as {credible_intervals}") - - # bit of defensive programming - if not isinstance(credible_intervals, list): - raise ValueError(f"Cannot set credible intervals to {credible_intervals}") - - for plotter in list(self._plotter_object_dict.values()): - if not isinstance(plotter, pt.posterior_base_classes._PosteriorPlottingBase): - continue - - # set credible intervals - plotter.credible_intervals = credible_intervals - - def set_variables_to_plot(self, plot_variables: List[str], plot_labels: List[str]=[])->None: - """Set variables to plot - - :param plot_variables: Variables to plot - :type plot_variables: List[str] - :param plot_labels: Plots we want to only plot these variables for, defaults to [] - :type plot_labels: List[str], optional - """ - for plotter in plot_labels: - self._plotter_object_dict[plotter].plot_params = plot_variables - - - def set_is_multimodal(self, param_ids: List[int | str])->None: - ''' - [Summary] - Lets posteriors know which parameters are multimodal - :param param_ids: list of multi-modal parameter ids/names - :type param_ids: List[str] - - ''' - print(f"Setting {param_ids} to be multimodal") - - # Loop over our plotters - for plotter in self._plotter_object_dict.values(): - if not isinstance(plotter, pt._posterior_plotting_base): - continue - - plotter.set_pars_multimodal(param_ids) - - - def set_is_circular(self, param_ids: List[int | str])->None: - ''' - [Summary] - Lets posteriors know which parameters are multimodal - - :param param_ids: list of multi-modal parameter ids/names - :type param_ids: List[str] - ''' - print(f"Setting {param_ids} to be circular") - - # Loop over our plotters - for plotter in self._plotter_object_dict.values(): - if not isinstance(plotter, pt.posterior_base_classes._PosteriorPlottingBase): - continue - - plotter.set_pars_circular(param_ids) - - def add_text_to_plots(self, text: str, location: tuple=(0.05, 0.95)): - ''' - [Summary] - Adds text to some plots - :param text: Text to add to plot - :param location: Location of text on plot - ''' - for plotter in self._plotter_object_dict.values(): - plotter.add_text_to_figures(text, location) - - def make_plots(self, output_file_name: str, plot_labels: List[str] | str): - ''' - [Summary] - Outputs all plots from a list of labels to an output PDF - :param output_file_name: Output file pdf name - :type output_file_name: str: - :param plot_labels: names of plots in self._plotter_object_dict - :type plot_labels: List[str] - ''' - # Cast our labels to hist - if not isinstance(plot_labels, list): - plot_labels = list(plot_labels) - - with PdfPages(output_file_name) as pdf_file: - for plotter_id in plot_labels: - try: - plotter_obj = self._plotter_object_dict.get(plotter_id) - except KeyError: - print(f"Warning:Key not found {plotter_id}, skipping") - continue - - print(f"Generating Plots for {plotter_id}") - plotter_obj.generate_all_plots() - print(f"Printing to {output_file_name}") - plotter_obj.write_to_pdf(existing_pdf_fig=pdf_file) - - def print_summary(self, latex_output_name:str=None): - ''' - [Summary] - Print stats summary to terminal and output as a LaTeX table [text file] - :param latex_output_name: name of output file, if empty or None doesn't print to file - :type latex_output_name: str, optional - ''' - summary = az.summary(self._chain_handler.arviz_tree, kind='stats', hdi_prob=0.9) - if latex_output_name is None: - return - - print(summary) - - with open(latex_output_name, "w") as output_file: - output_file.write(summary.to_latex()) - \ No newline at end of file diff --git a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/__init__.py b/src/MaCh3PythonUtils/diagnostics/mcmc_plots/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/diagnostics/__init__.py b/src/MaCh3PythonUtils/diagnostics/mcmc_plots/diagnostics/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/diagnostics/autocorrelation_trace_plotter.py b/src/MaCh3PythonUtils/diagnostics/mcmc_plots/diagnostics/autocorrelation_trace_plotter.py deleted file mode 100644 index abc30fa..0000000 --- a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/diagnostics/autocorrelation_trace_plotter.py +++ /dev/null @@ -1,45 +0,0 @@ -''' -HI : Class to make autocorrelations and traces, puts them all onto a single plot -''' -import arviz as az -from matplotlib import pyplot as plt -from MaCh3PythonUtils.diagnostics.mcmc_plots.plotter_base import _PlottingBaseClass -from matplotlib.figure import Figure - -class AutocorrelationTracePlotter(_PlottingBaseClass): - def _generate_plot(self, parameter_name: str) -> Figure: - """Generates auto-correlation and trace plots for a single parameter - - :param parameter_name: Name of parameter to plot - :type parameter_name: str - :return: Figure containing auto-correlation and trace plots - :rtype: Figure - """ - fig, (trace_ax, autocorr_ax) = plt.subplots(nrows=2, sharex=False) - - # We want the numpy array containing our parameter - param_array = self._chain_handler.arviz_tree[parameter_name].to_numpy()[0] - - # Okay now we can plot our trace (might as well!) - trace_ax.plot(param_array, linewidth=0.05, color='purple') - - #next we need to grab our autocorrelations - # auto_corr = sm.tsa.acf(param_array, nlags=total_lags) - auto_corr = az.autocorr(param_array) - autocorr_ax.plot(auto_corr, color='purple') - # Now we do some tidying - trace_ax.set_ylabel(f"{parameter_name} variation", fontsize=11) - trace_ax.set_xlabel("Step", fontsize=11) - trace_ax.set_title(f"Trace for {parameter_name}", fontsize=15) - trace_ax.tick_params(labelsize=10) - # autocorr_ax.set_ylabel("Autocorrelation Function") - autocorr_ax.set_xlabel("Lag", fontsize=12) - autocorr_ax.set_ylabel("Autocorrelation", fontsize=12) - autocorr_ax.set_title(f"Autocorrelation for {parameter_name}", fontsize=15, verticalalignment='center_baseline') - autocorr_ax.tick_params(labelsize=10) - fig.suptitle(f"{parameter_name} diagnostics")#, fontsize=40) - # fig.subplots_adjust(hspace=0.01) - fig.subplots_adjust(wspace=0.0, hspace=0.3) - - plt.close() - return fig diff --git a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/diagnostics/covariance_matrix_utils.py b/src/MaCh3PythonUtils/diagnostics/mcmc_plots/diagnostics/covariance_matrix_utils.py deleted file mode 100644 index 25a4c6b..0000000 --- a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/diagnostics/covariance_matrix_utils.py +++ /dev/null @@ -1,120 +0,0 @@ -''' -Additional Diagnostics that can be used with MCMC but don't rely on plotting - -Suboptimality : https://www.jstor.org/stable/25651249?seq=3 -''' - -from MaCh3PythonUtils.file_handling.chain_handler import ChainHandler -import numpy as np -from scipy.linalg import sqrtm -from tqdm.auto import tqdm -from matplotlib import pyplot as plt -from matplotlib.backends.backend_pdf import PdfPages -import warnings -import mplhep as hep -from concurrent.futures import ThreadPoolExecutor, as_completed - - -class CovarianceMatrixUtils: - def __init__(self, chain_handler: ChainHandler)->None: - """Constructor for handling covariance matrix - - :param chain_handler: Instance of a chain handler objec - :type chain_handler: ChainHandler - """ - # Let's just ignore some warnings :grin: - warnings.filterwarnings("ignore", category=DeprecationWarning) #Some imports are a little older - warnings.filterwarnings("ignore", category=UserWarning) #Some imports are a little older - - chain_handler.get_ttree_as_arviz() # Ensure that we have things in the correct format - self._parameter_names = list(chain_handler.arviz_tree.keys()) - self._total_parameters = len(self._parameter_names) - self._ttree_data_frame = chain_handler.arviz_tree.to_dataframe() # All our handling will be done using thi s object - # Various useful class properties - self._suboptimality_array = [] # For filling with suboptimality - self._suboptimality_evaluation_points = [] #Where did we evalutate this? - self._full_covariance = self.calculate_covariance_matrix() - self._sqrt_full_covariance_inv = np.linalg.inv(sqrtm(self._full_covariance)) - plt.style.use(hep.style.ROOT) - - def calculate_covariance_matrix(self, min_step: int=0, max_step: int=-1)->np.ndarray: - """_summary_ - - :param min_step: Steps at which to calculate covariance matrix at, defaults to 0 - :type min_step: int, optional - :param max_step: Maximum step to calculate covariance matrix at, defaults to -1 indicating you want to evaluate the full chain - :type max_step: int, optional - :return: Covariance matrix at each Nth step - :rtype: np.ndarray - """ - if(max_step<=0): - sub_array = self._ttree_data_frame[min_step : ] - else: - sub_array = self._ttree_data_frame[min_step : max_step] - - - return sub_array.cov().to_numpy() - - def _calculate_matrix_suboptimality(self, step_number: int)->float: - """Calculate the suboptimality value for the mtarix for step N - - :param step_number: calculate the covariance matrix for the first step_number steps - :type step_number: int - :return: Suboptimality - :rtype: float - """ - new_covariance = self.calculate_covariance_matrix(max_step=step_number) - - sqrt_input_cov = sqrtm(new_covariance) - - # Get product of square roots - matrix_prod = sqrt_input_cov @ self._sqrt_full_covariance_inv - - #Get eigen values - eigenvalues, _ = np.linalg.eig(matrix_prod) - - return self._total_parameters * np.sum(eigenvalues**(-2))/((np.sum(eigenvalues)**(-1))**2) - - def calculate_suboptimality(self, step_skip : int = 1000, min_step=0)->None: - """ - Calculates the suboptimality for every step_skip steps - - :param step_skip: Number of steps to skip, defaults to 1000 - :type step_skip: int, optional - :param min_step: smallest number of starting steps, defaults to 0 - :type min_step: int, optional - """ - self._suboptimality_evaluation_points = np.arange(min_step, len(self._ttree_data_frame), step_skip) - # Make sure we have the last step as well! - self._suboptimality_evaluation_points = np.append(self._suboptimality_evaluation_points, len(self._ttree_data_frame)-1) - # make array to fill with suboptimality values - self._suboptimality_array = np.zeros(len(self._suboptimality_evaluation_points)) # Reset our suboptimality values - print(f"Calculating suboptimality for {len(self._suboptimality_evaluation_points)} points") - - # Lets speed this up a bit! - with ThreadPoolExecutor() as executor: - futures = {executor.submit(self._calculate_matrix_suboptimality, step): step - for step in self._suboptimality_evaluation_points} - - for future in tqdm(as_completed(futures), ascii="▖▘▝▗▚▞█", total=len(self._suboptimality_evaluation_points)): - i = np.where(self._suboptimality_evaluation_points==futures[future])[0] - self._suboptimality_array[i]=future.result() - - - def plot_suboptimality(self, output_file: str): - """Plots suboptimality value - - :param output_file: Output PDF file - :type output_file: str - """ - print(f"Saving to {output_file}") - with PdfPages(output_file) as pdf: - fig, axes = plt.subplots() - - axes.plot(self._suboptimality_evaluation_points, self._suboptimality_array) - axes.set_xlabel("Step Number") - axes.set_label("Suboptimality") - axes.set_title("Suboptimality Plot") - axes.set_yscale('log') - - pdf.savefig(fig) diff --git a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/diagnostics/simple_diag_plots.py b/src/MaCh3PythonUtils/diagnostics/mcmc_plots/diagnostics/simple_diag_plots.py deleted file mode 100644 index 4ada1c4..0000000 --- a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/diagnostics/simple_diag_plots.py +++ /dev/null @@ -1,57 +0,0 @@ -''' -HI : Several simple diagnostics -''' - -from typing import List -import arviz as az -from matplotlib import pyplot as plt -from MaCh3PythonUtils.diagnostics.mcmc_plots.plotter_base import _PlottingBaseClass -from matplotlib.figure import Figure - -class EffectiveSampleSizePlotter(_PlottingBaseClass): - def _generate_plot(self, parameter_name: str) -> Figure: - """Generates ESS plot for single parameter - - :param parameter_name: Name of parameter to plot - :type parameter_name: str - :return: ESS plot for parameter_name - :rtype: Figure - """ - fig, axes = plt.subplots() - - az.plot_ess(self._chain_handler.arviz_tree, var_names=parameter_name, - ax=axes, textsize=30, color='purple', drawstyle="steps-mid", linestyle="-") - plt.close() - return fig - -class MarkovChainStandardError(_PlottingBaseClass): - def _generate_plot(self, parameter_name: str) -> Figure: - """Generates MCSE plot for single parameter - - :param parameter_name: Name of parameter to plot - :type parameter_name: str - :return: ESS plot for parameter_name - :rtype: Figure - """ - fig, axes = plt.subplots() - az.plot_mcse(self._chain_handler.arviz_tree, var_names=parameter_name, ax=axes, textsize=10, color='purple') - plt.close() - return fig - - -class ViolinPlotter(_PlottingBaseClass): - def _generate_plot(self, parameter_name: str | List[str]) -> Figure: - """Generates Violin plot for single parameter - - :param parameter_name: Name of parameter to plot - :type parameter_name: str - :return: ESS plot for parameter_name - :rtype: Figure - """ - # total number of axes we need - fig, axes = plt.subplots() - az.plot_violin(self._chain_handler.arviz_tree, - var_names=parameter_name, ax=axes, textsize=10, - shade_kwargs={'color':'purple'}) - plt.close() - return fig diff --git a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/plotter_base.py b/src/MaCh3PythonUtils/diagnostics/mcmc_plots/plotter_base.py deleted file mode 100644 index 4af10de..0000000 --- a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/plotter_base.py +++ /dev/null @@ -1,226 +0,0 @@ -''' -HI : Mostly abstract base class, contains common methods for use by most other plotting classes -''' - -from MaCh3PythonUtils.file_handling.chain_handler import ChainHandler -import arviz as az -from matplotlib import pyplot as plt -from matplotlib.figure import Figure -from matplotlib.backends.backend_pdf import PdfPages -import numpy as np -import warnings -from abc import ABC, abstractmethod -from typing import List -import mplhep as hep -from tqdm.auto import tqdm -from concurrent.futures import ThreadPoolExecutor, as_completed -import numpy.typing as npt -from typing import Any -from collections.abc import Iterable - - -# Base class with common methods -class _PlottingBaseClass(ABC): - def __init__(self, chain_handler: ChainHandler)->None: - """Abstract plotting class for generating plots - - :param chain_handler: ChainHandler instance - :type chain_handler: ChainHandler - """ # I'm very lazy but this should speed up parallelisation since we won't have this warning - plt.rcParams.update({'figure.max_open_warning': 0}) - - self._chain_handler=chain_handler # Make sure we don't use crazy amount of memory - # Setup plotting styles - az.style.use(hep.style.ROOT) - plt.style.use(hep.style.ROOT) - - self._figure_list = np.empty([]) # List of all figures - - # Let's just ignore some warnings :grin: - self._all_plots_generated = False # Stops us making plots multiple times - - # Setting it like this replicates the old behaviour! - self._circular_params=np.ones(len(self._chain_handler.arviz_tree)).astype(bool) - warnings.filterwarnings("ignore", category=DeprecationWarning) #Some imports are a little older - warnings.filterwarnings("ignore", category=UserWarning) #Some imports are a little older - - az.utils.Dask().enable_dask(dask_kwargs={"dask": "parallelized", "output_dtypes": [float]}) # let it be a bit cleverer - az.utils.Numba().enable_numba() - - # Do we want figure text? - self._figure_text = None - self._text_location = (0.05, 0.85) - - # Default option - self._params_to_plot = list(self._chain_handler.arviz_tree.keys()) - - def __str__(self) -> str: - return "plotting_base_class" - - @property - def figure_list(self)->npt.NDArray[Any]: - """Get list of figures - - :return: List of figures generated - :rtype: npt.NDArray[Any] - """ - return self._figure_list - - @figure_list.setter - def figure_list(self, new_figure)->None: - raise NotImplementedError("Cannot set figure list using property!") - - @abstractmethod - def _generate_plot(self, parameter_name: str | List[str])->None: - """Abstract method for generating plot - - :param parameter_name: Name/Names of parameters wanted for plotting - :type parameter_name: str | List[str] - """ - # Abstract method to generate a single plot - pass - - def generate_all_plots(self)->None: - """Generates all plots requested by user - """ - # Generates plots for every parameter [can be overwritten] - if self._all_plots_generated: # No need to make all possible plots again! - return - - self._figure_list = np.empty(len(self._params_to_plot), dtype=Figure) - - # Parallelised loop - with ThreadPoolExecutor() as executor: - # Set of threadpools - futures = {executor.submit(self._generate_plot, param) : param for param in self._params_to_plot} - # Begin loop - for future in tqdm(as_completed(futures), ascii="▖▘▝▗▚▞█", total=len(self._params_to_plot)): - param_id = list(futures).index(future) # Unique plot ID - self._figure_list[param_id] = future.result() - plt.close() - - self._all_plots_generated=True - - # Lets us select a subset/add list of parameters we'd like to plot - @property - def plot_params(self)-> List[str] | List[List[str]]: - """Get all parameters that are going to be plotted - - :return: List of parameters/groups of parameteres used in each plot - :rtype: List[str] | List[List[str]] - """ - return self._params_to_plot - - @plot_params.setter - def plot_params(self, new_plot_parameter_list: List[str]|List[List[str]]): - ''' - Allows for the addition of new parameters - ''' - if len(new_plot_parameter_list)==0: - raise ValueError("Parameter list cannot have length 0") - - # Check our new parameters are in our list of keys - if isinstance(new_plot_parameter_list[0], str): - for parameter in new_plot_parameter_list: - self._parameter_not_found_error(parameter) - - - elif isinstance(new_plot_parameter_list[0][0], str): - for param_list in new_plot_parameter_list: - for param in param_list: - self._parameter_not_found_error(param) - - else: - raise ValueError("Plot params must be of type List[str] or List[List[str]]") - - self._params_to_plot = new_plot_parameter_list - - - - def _parameter_not_found_error(self, parameter_name: str): - if parameter_name not in list(self._chain_handler.arviz_tree.keys()): - raise ValueError(f"{parameter_name} not in list of parameters!") - - def _get_param_index_from_name(self, parameter_name: str)->int: - """Gets parameter index in list based on parameter name - - :param parameter_name: Parameter name - :type parameter_name: str - :return: Parameter index in self._params_to_plot - :rtype: int - """ - # Gets index of parameter in our arviz array - - self._parameter_not_found_error(parameter_name) - param_id = list(self._chain_handler.arviz_tree.keys()).index(parameter_name) - return param_id - - def set_pars_circular(self, par_id_list: List[str] | List[int])->None: - """Sets parameter to be cyclical about -pi, pi - - :param par_id_list: List of a parameter names/parameter ids - :type par_id_list: List[str] | List[int] - """ - # If it's an int this is easy - if not isinstance(par_id_list, list): - par_id_list = list(par_id_list) - - for par_id in par_id_list: - if isinstance(par_id, int): - true_index = par_id - else: - true_index = self._get_param_index_from_name(par_id) - - self._circular_params[true_index] = True - - def add_text_to_figures(self, text: str, text_location: tuple=(0.05, 0.95))->None: - """Add text to figures - - :param text: Text to be added to figure - :type text: str - :param text_location: Location of text on figures, defaults to (0.05, 0.95) - :type text_location: tuple, optional - """ - self._figure_text = text - self._text_location = text_location - - def write_to_pdf(self, output_pdf_name: str=None, existing_pdf_fig: PdfPages=None)->None: - """Dump all plots to PDF file - - :param output_pdf_name: Output PDF file, defaults to None - :type output_pdf_name: str, optional - :param existing_pdf_fig: Write to existing PDF instance, defaults to None - :type existing_pdf_fig: PdfPages, optional - :raises ValueError: No file used - """ ''' - Dump all our plots to PDF file must either set output name or existing_pdf_fig - - inputs : - output_pdf_name: [type=string] Output name for NEW pdf - existing_pdf_fig: [type=PDFPages] - returns : - pdf_fig: [type=PDFPages] pdf file reader [REMEMBER TO CLOSE THE PDF FILE AT THE END!] - ''' - if len(self._figure_list)==0: - return - - if output_pdf_name is not None and existing_pdf_fig is None: - pdf_file = PdfPages(output_pdf_name) - elif existing_pdf_fig is not None: - pdf_file = existing_pdf_fig - else: - raise ValueError("ERROR:Must set EITHER output_pdf_name OR existing_pdf_fig") - - # For some arrays we might want to make them 1D - if not isinstance(self._figure_list[0], Figure): - self._figure_list = [fig for sublist in self._figure_list for fig in sublist] - - - for fig in tqdm(self._figure_list, ascii=" ▖▘▝▗▚▞█"): - # Add text to all plots! - if self._figure_text is not None: - if len(fig.get_axes())==1: - fig.text(*self._text_location, self._figure_text, transform=fig.get_axes()[0].transAxes, fontsize=60, fontstyle='oblique') - - pdf_file.savefig(fig) - plt.close() diff --git a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/posteriors/__init__.py b/src/MaCh3PythonUtils/diagnostics/mcmc_plots/posteriors/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/posteriors/posterior_base_classes.py b/src/MaCh3PythonUtils/diagnostics/mcmc_plots/posteriors/posterior_base_classes.py deleted file mode 100644 index 5650216..0000000 --- a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/posteriors/posterior_base_classes.py +++ /dev/null @@ -1,68 +0,0 @@ -''' -HI : Contains a set of base classes with methods common to posterior plotting code -Separated into 1D and 2D-like objects to make life easier -''' -from MaCh3PythonUtils.file_handling.chain_handler import ChainHandler -from MaCh3PythonUtils.diagnostics.mcmc_plots.plotter_base import _PlottingBaseClass -from typing import List -import numpy as np -from tqdm.auto import tqdm -import numpy.typing as npt - - -# Base class for all posterior plotters -class _PosteriorPlottingBase(_PlottingBaseClass): - # Small extension of _plotting_base class for posterior specific stuff - def __init__(self, file_loader: ChainHandler)->None: - """Base class for all posterior plotters - - :param file_loader: ChainHandler instance - :type file_loader: ChainHandler - """ - # Setup additional features for posteriors - super().__init__(file_loader) - self._credible_intervals = np.array([0.6, 0.9, 0.95]) - self._parameter_multimodal = np.zeros(len(self._file_loader.arviz_tree)).astype(bool) - - def __str__(self) -> str: - return "posterior_plotting_base_class" - - @property - def credible_intervals(self)->List[float]: - """Get Credible Intervals - - :return: List of credible intervals - :rtype: List[float] - """ - return self._credible_intervals - - @credible_intervals.setter - def credible_intervals(self, new_creds: List[float])->None: - """Set new credible intervals - - :param new_creds: Sets list of credible intervals - :type new_creds: List[float] - """ - # Sets credible intervals from list - self._credible_intervals = np.array(new_creds) # Make sure it's 1D - self._credible_intervals.sort() # Flatten it - - def set_pars_multimodal(self, par_id_list: List[str] | List[int], is_multimodal: bool=True)->None: - """Let arviz know which parameters are multi-modal - - :param par_id_list: List of multimodal parameter names/ids - :type par_id_list: List[str] | List[int] - :param is_multimodal: Sets parameters to either be multimodal [True] or uni-modal [False], defaults to True - :type is_multimodal: bool, optional - """ - # If it's an int this is easy - if not isinstance(par_id_list, list): - par_id_list = list(par_id_list) - - for par_id in par_id_list: - if isinstance(par_id, int): - true_index = par_id - else: - true_index = self._get_param_index_from_name(par_id) - - self._parameter_multimodal[true_index] = is_multimodal diff --git a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/posteriors/posteriors_1d.py b/src/MaCh3PythonUtils/diagnostics/mcmc_plots/posteriors/posteriors_1d.py deleted file mode 100644 index 934db49..0000000 --- a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/posteriors/posteriors_1d.py +++ /dev/null @@ -1,87 +0,0 @@ -import arviz as az -from matplotlib import pyplot as plt -import numpy as np -from MaCh3PythonUtils.diagnostics.mcmc_plots.posteriors.posterior_base_classes import _PosteriorPlottingBase -from matplotlib.figure import Figure - -''' -Plotting class for 1D plots. Slightly overcomplicated by various special cases but works a treat! -''' -class PosteriorPlotter1D(_PosteriorPlottingBase): - def _generate_plot(self, parameter_name: str) -> Figure: - """Generate 1D posterior plot for a single parameter - - :param parameter_name: Name of parameter to plot - :type parameter_name: str - :raises ValueError: if parameter is not a string - :return: 1D Posterior plot - :rtype: Figure - """ ''' - Generates a single posterior plot for a parameter - - inputs : - parameter name : [type=str] Name of parameter - ''' - if not isinstance(parameter_name, str): - raise ValueError(f"Can only pass single parameters to posterior plotting class. Cannot plot {parameter_name}") - - # Checks if parameter is in our array+gets the index - param_index = self._get_param_index_from_name(parameter_name) - - fig, axes = plt.subplots(figsize=(30, 30)) - # Set bin number and range - n_bins = 50 - # Grab our density plot - # plt.rcParams['image.cmap'] = 'Purples' # To avoid any multi-threaded weirdness - line_colour_generator = iter(plt.cm.Purples(np.linspace(0.4, 0.8, len(self.credible_intervals)+1))) - line_colour = next(line_colour_generator) - hist_kwargs={'density' : True, 'bins' : n_bins, "alpha": 1.0, 'linewidth': None, 'edgecolor': line_colour, 'color': 'white'} - - _, bins, patches = axes.hist(self._file_loader.arviz_tree[parameter_name].to_numpy()[0], **hist_kwargs) - # Make lines for each CI - - cred_rev = self.credible_intervals[::-1] - for credible in cred_rev: - line_colour = next(line_colour_generator) - - # We want the bayesian credible interval, for now we set the maximum number of modes for multi-modal parameters to 2 - hdi = az.hdi(self._file_loader.arviz_tree, var_names=[parameter_name], hdi_prob=credible, - multimodal=self._parameter_multimodal[param_index], max_modes=20, circular=self._circular_params[param_index]) - - # Might be multimodal so we want all our credible intervals in a 1D array to make plotting easier! - credible_bounds = hdi[parameter_name].to_numpy() - - if isinstance(credible_bounds[0], float): - credible_bounds = np.array([credible_bounds]) #when we're not multi-modal we need to be a little careful - - # set up credible interval - plot_label = f"{100*credible}% credible interval " - for bound in credible_bounds: - # Set up plotting options - # Reduce the plotting array to JUST be between our boudnaries - mask = (bins>=bound[0]) & (bins<=bound[1]) - if bound[0]>bound[1]: - # We need a SLIGHTLY different treatment since we loop around for some parameters (delta_cp) - mask = (bins>=bound[0]) | (bins<=bound[1]) - - for patch_index in np.where(mask)[0]: - patches[patch_index-1].set_facecolor(line_colour) - patches[patch_index-1].set_edgecolor(None) - patches[patch_index-1].set_label(plot_label) - - # add legend - # Set Some labels - axes.set_xlabel(parameter_name, fontsize=50) - axes.tick_params(labelsize=40) - axes.set_title(f"Posterior Density for {parameter_name}", fontsize=60) - axes.set_ylim(ymin=0) - axes.set_ylabel("Posterior Density", fontsize=50) - - # Generate Unique set of labels and handles - plot_handles, plot_labels = axes.get_legend_handles_labels() - # Sometimes it loses track of the posterior histogram - unique_labs = [(h, l) for i, (h, l) in enumerate(zip(plot_handles, plot_labels)) if l not in plot_labels[:i]] - axes.legend(*zip(*unique_labs), loc=(0.6,0.85), fontsize=35, facecolor="white", edgecolor="black", frameon=True) - #Stop memory issues - plt.close() - return fig # Add to global figure list diff --git a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/posteriors/posteriors_2d.py b/src/MaCh3PythonUtils/diagnostics/mcmc_plots/posteriors/posteriors_2d.py deleted file mode 100644 index cb604d1..0000000 --- a/src/MaCh3PythonUtils/diagnostics/mcmc_plots/posteriors/posteriors_2d.py +++ /dev/null @@ -1,106 +0,0 @@ -''' -HI : Makes 2D posterior plots. Currently no way of putting a legend on the plot (thanks arviz...) -''' -import arviz as az -from typing import List, Any -from matplotlib import pyplot as plt -from itertools import combinations -from MaCh3PythonUtils.diagnostics.mcmc_plots.posteriors.posterior_base_classes import _PosteriorPlottingBase -from MaCh3PythonUtils.file_handling.chain_handler import ChainHandler -from matplotlib.figure import Figure -import numpy as np -import numpy.typing as npt - -class PosteriorPlotter2D(_PosteriorPlottingBase): - def _generate_plot(self, parameter_names: List[str]) -> npt.NDArray[Any]: - """Generates 2D posterior plot - - :param parameter_names: List of parameters to plot [plots all combinations of provided parameters] - :type parameter_names: List[str] - :return: List of figures containing plots - :rtype: npt.NDArray[Any] - """ - # Let's get pairs of names - name_pairs_list = list(combinations(parameter_names, 2)) - fig_list = np.empty(len(name_pairs_list), dtype=Figure) - # Now we loop over our pairs of names - for i, (par_1, par_2) in enumerate(name_pairs_list): - - fig, axes = plt.subplots(figsize=(30, 30)) - - par_1_numpy_arr = self._file_loader.arviz_tree[par_1].to_numpy()[0] - par_2_numpy_arr = self._file_loader.arviz_tree[par_2].to_numpy()[0] - - ciruclar = self._circular_params[self._get_param_index_from_name(par_1)] | self._circular_params[self._get_param_index_from_name(par_2)] - - az.plot_kde(par_1_numpy_arr, par_2_numpy_arr, - hdi_probs= self._credible_intervals, - contourf_kwargs={"cmap": "Purples"}, - ax=axes, is_circular=ciruclar, legend=True) - - axes.set_xlabel(par_1, fontsize=50) - axes.set_ylabel(par_2, fontsize=50) - axes.tick_params(labelsize=40) - - - cred_as_str=",".join(f"{100*i}%" for i in self._credible_intervals) - axes.set_title(f"{par_1} vs {par_2} : [{cred_as_str}] credible intervals", fontsize="60") - plt.close() # CLose the canvas - fig_list[i] = fig - return fig_list - - - -class TrianglePlotter(_PosteriorPlottingBase): - # Makes triangle plots - def __init__(self, file_loader: ChainHandler)->None: - ''' - Constructor - ''' - # Inherit the abstract base class - super().__init__(file_loader) - - def _generate_plot(self, parameter_names: List[str]) -> Figure: - ''' - Generates a single triangle plot and adds it to the figure list - inputs : - -> parameter_names : [List[str]] List of parameter names to put in the triangle - returns : - -> Figure - ''' - if not isinstance(parameter_names, list): - raise ValueError("Parameter names must be list when plotting triangle plots") - - # Check we are using a valid parameter set - for param in parameter_names: - self._parameter_not_found_error(param) # Check if parameters exist - - fig, axes = plt.subplots(nrows=len(parameter_names), ncols=len(parameter_names), figsize=(30, 30)) - - az.plot_pair(self._file_loader.arviz_tree, var_names=parameter_names, - marginals=True, ax=axes, colorbar=True, figsize=(30,30), - kind='kde', - textsize=30, - kde_kwargs={ - "hdi_probs": self._credible_intervals, # Plot HDI contours - "contourf_kwargs": {"cmap": "Purples"}, - 'legend':True - }, - marginal_kwargs={ - # 'fill_kwargs': {'alpha': 0.0}, - 'plot_kwargs': {"linewidth": 4.5, "color": "purple"}, - # "quantiles": credible_intervals, - "rotated": False - }, - point_estimate='mean', - ) - - cred_as_str=",".join(f"{100*i}%" for i in self._credible_intervals) - fig.suptitle(f"Triangle Plot for : {cred_as_str} credible intervals", fontsize="60") - - # axes[-1, -1].legend(axes, [f"{i}% Credible Interval" for i in self._credible_intervals], frameon=True, loc='right') - - fig.subplots_adjust(wspace=0.01, hspace=0.01) - - plt.close() - return fig \ No newline at end of file diff --git a/src/MaCh3PythonUtils/file_handling/chain_handler.py b/src/MaCh3PythonUtils/file_handling/chain_handler.py index 5ffbc3e..997b5ca 100644 --- a/src/MaCh3PythonUtils/file_handling/chain_handler.py +++ b/src/MaCh3PythonUtils/file_handling/chain_handler.py @@ -8,7 +8,6 @@ from concurrent.futures import ThreadPoolExecutor import gc import numpy as np -import arviz as az from numpy.typing import NDArray class ChainHandler: @@ -209,19 +208,7 @@ def ttree_array(self, new_array: Any=None)->None: :type new_array: Any ''' # Implemented in case someone tries to do something daft! - raise NotImplementedError("Cannot set converted TTree array to new type") - - def make_arviz_tree(self): - """Generate Arivz tree - - :raises RuntimeError: Not built TTree - """ - if self._ttree_array is None: - raise RuntimeError("Error have not converted ROOT TTree to pandas data frame yet!") - - print("Generating Arviz data struct [this may take some time!]") - self._arviz_tree = az.dict_to_dataset(self._ttree_array.to_dict(orient='list')) - + raise NotImplementedError("Cannot set converted TTree array to new type") @property def ndim(self)->int: @@ -243,14 +230,4 @@ def upper_bounds(self)->NDArray: if self._ttree_array is None: return np.empty(self.ndim) - return self._ttree_array.max(axis=0).to_numpy() - - - @property - def arviz_tree(self)->az.InferenceData: - """Gets arviz version of dataframe - - :return: Arviz version of tree - :rtype: arviz.InferenceData - """ - return self._arviz_tree \ No newline at end of file + return self._ttree_array.max(axis=0).to_numpy() \ No newline at end of file diff --git a/src/MaCh3PythonUtils/fitters/adaption_handler.py b/src/MaCh3PythonUtils/fitters/adaption_handler.py deleted file mode 100644 index a38426f..0000000 --- a/src/MaCh3PythonUtils/fitters/adaption_handler.py +++ /dev/null @@ -1,68 +0,0 @@ -import numpy as np -from numba import njit - -@njit -def update_covariance(mean, covariance, new_data, count): - """Update the mean and covariance matrix efficiently.""" - delta = new_data - mean - mean += delta / count - - delta2 = new_data - mean - covariance += np.outer(delta, delta2) * (count - 1) / count - - return mean, covariance - - -class CovarianceUpdater: - def __init__(self, n_dimensions, update_step=1000): - self.n = n_dimensions - self.mean = np.zeros(n_dimensions) - self.covariance = np.zeros((n_dimensions, n_dimensions)) - self.frozen_covariance = np.identity(n_dimensions) - - self.count = 0 - self.update_step=update_step - self._rng = np.random.default_rng() - - - def update(self, new_data): - """Update the mean and covariance with a new data point. - - Parameters: - new_data (np.ndarray): A 1D array representing the new data point (shape: (n_dimensions,)) - """ - # if new_data.shape[0] != self.n: - # raise ValueError(f"Expected data with {self.n} dimensions, got {new_data.shape[0]}.") - - self.count += 1 - - # Update mean and covariance using the Numba-optimized function - self.mean, self.covariance = update_covariance(self.mean, self.covariance, new_data, self.count) - - if self.count%self.update_step and self.count>100: - self.frozen_covariance = self.get_covariance().copy() - - def get_mean(self): - """Return the current mean.""" - return self.mean - - def get_covariance(self): - """Return the current covariance matrix.""" - return self.covariance / (self.count - 1) if self.count > 1 else np.zeros((self.n, self.n)) - - def get_frozen_covariance(self): - return (2.4**2)/(float(self.n)) * self.frozen_covariance + np.identity(self.n)*1e-6 - - def sample(self, n_samples: int): - """Sample from a multivariate Gaussian centered at 'center'. - - Parameters: - center (np.ndarray): A 1D array representing the center for the Gaussian distribution (shape: (n_dimensions,)) - num_samples (int): The number of samples to draw. - - Returns: - np.ndarray: Samples drawn from the multivariate Gaussian distribution. - """ - return self._rng.multivariate_normal(mean=np.zeros(self.n), size=n_samples, - cov=self.get_frozen_covariance(), tol=1e-6) - diff --git a/src/MaCh3PythonUtils/fitters/mcmc.py b/src/MaCh3PythonUtils/fitters/mcmc.py deleted file mode 100644 index 788ec91..0000000 --- a/src/MaCh3PythonUtils/fitters/mcmc.py +++ /dev/null @@ -1,114 +0,0 @@ -import emcee -import numpy as np -from numpy.typing import NDArray -from tqdm import tqdm - -import matplotlib.pyplot as plt -from matplotlib.backends.backend_pdf import PdfPages -import h5py -from typing import List - - -from MaCh3PythonUtils.machine_learning.file_ml_interface import FileMLInterface -from MaCh3PythonUtils.fitters.adaption_handler import CovarianceUpdater - -class MCMC: - def __init__(self, interface: FileMLInterface, circular_params: List[str]=[], update_step: int=10, start_matrix_throw: int=0): - print("MCMC let's go!") - - self._interface = interface - self._n_dim = interface.chain.ndim - 1 - self._n_chains = 1 - - # HACK, we fit to scaled rather than anything else because I AM LAZY (and slight efficiency saving since we invert this at the end) - self._upper_bounds = interface.chain.upper_bounds[:-1] - self._lower_bounds = interface.chain.lower_bounds[:-1] - self._sampler = None - - # Get initial state etc. - self._chain_state: NDArray = self._interface.training_data.iloc[[1]].to_numpy()[0] - - - self._circular_indices = [self._interface.chain.plot_branches.index(par) for par in circular_params] - self._matrix_scale = 2.4**2/float(self._n_dim) - self._start_matrix_throw = start_matrix_throw - self._update_step = update_step - - self._matrix_handler = CovarianceUpdater(self._n_dim, update_step) - - self._current_loglikelihood = self._calc_likelihood(self._chain_state) - self._current_step = 0 - - def _wrap_circular(self, state): - return (state + np.pi) % (2 * np.pi) - np.pi - - def _calc_likelihood(self, state: NDArray): - # state[self._circular_indices] = self._wrap_circular(state[self._circular_indices]) - - # if np.any(state>self._upper_bounds) or np.any(state -1*float(np.inf): - u = np.random.uniform(0, 1) - - if proposed_loglikelihood>self._current_loglikelihood: - acc_prob = 1 - else: - acc_prob = np.exp(proposed_loglikelihood-self._current_loglikelihood) - - if min(1, acc_prob)>u: - self._chain_state = proposed_state - self._current_loglikelihood = proposed_loglikelihood - - self._matrix_handler.update(self._chain_state[0]) - - - self._dataset[self._current_step]=self._chain_state - self._current_step += 1 - - - def save_mcmc_chain_to_pdf(self, filename: str, output_pdf: str): - # Open the HDF5 file - with h5py.File(filename, 'r') as f: - # Load the chain dataset - chain = f['chain'][:] - - _, n_params = chain.shape - - # Create a PdfPages object to save plots - with PdfPages(output_pdf) as pdf: - for i in tqdm(range(n_params)): - fig, ax = plt.subplots(figsize=(10, 6)) - - # Plot the chain for the i-th parameter - ax.plot(chain[:, i], lw=0.5) - ax.set_ylabel(self._interface.chain.plot_branches[i]) - ax.set_title(f"Parameter {self._interface.chain.plot_branches[i]} MCMC Chain") - ax.set_xlabel('Step') - - # Save the current figure to the PDF - pdf.savefig(fig) - plt.close(fig) # Close the figure to save memory - - print(f"MCMC chain plots saved to {output_pdf}") - - - def __call__(self, n_steps, output_file_name: str): - print(f"Running MCMC for {n_steps} steps") - - # Open the HDF5 file and create the dataset with the correct shape upfront - with h5py.File(output_file_name, 'w') as f: - self._dataset = f.create_dataset('chain', (n_steps, self._n_dim), chunks=True) - - for _ in tqdm(range(n_steps)): - self.propose_step() - - self.save_mcmc_chain_to_pdf(output_file_name, "traces.pdf") \ No newline at end of file diff --git a/src/MaCh3PythonUtils/fitters/multi_mcmc.py b/src/MaCh3PythonUtils/fitters/multi_mcmc.py deleted file mode 100644 index d5ae909..0000000 --- a/src/MaCh3PythonUtils/fitters/multi_mcmc.py +++ /dev/null @@ -1,116 +0,0 @@ -import emcee -import numpy as np -from numpy.typing import NDArray -from tqdm import tqdm -import matplotlib.pyplot as plt -from matplotlib.backends.backend_pdf import PdfPages -import h5py -from typing import List -from MaCh3PythonUtils.machine_learning.file_ml_interface import FileMLInterface -from MaCh3PythonUtils.fitters.adaption_handler import CovarianceUpdater - -class MCMC: - def __init__(self, interface: FileMLInterface, n_chains: int = 1024, circular_params: List[str] = [], update_step: int = 10, start_matrix_throw: int = 0): - print("MCMC let's go!") - - self._interface = interface - self._n_dim = interface.chain.ndim - 1 - self._n_chains = n_chains - - # HACK, we fit to scaled rather than anything else because I AM LAZY (and slight efficiency saving since we invert this at the end) - self._upper_bounds = interface.chain.upper_bounds[:-1] - self._lower_bounds = interface.chain.lower_bounds[:-1] - self._sampler = None - - # Get initial states for all chains - initial_state = self._interface.training_data.iloc[[1]].to_numpy()[0] - self._chain_states: NDArray = np.full((n_chains, self._n_dim),initial_state) # Create a state for each chain - - self._circular_indices = [self._interface.chain.plot_branches.index(par) for par in circular_params] - self._matrix_scale = 2.4**2 / float(self._n_dim) - self._start_matrix_throw = start_matrix_throw - self._update_step = update_step - - # Use the first chain to update the covariance matrix - self._matrix_handler = CovarianceUpdater(self._n_dim, update_step) - - # Calculate likelihoods for all chains - self._current_loglikelihoods = self._calc_likelihood(self._chain_states) - self._current_step = 0 - - def _wrap_circular(self, state): - return (state + np.pi) % (2 * np.pi) - np.pi - - def _calc_likelihood(self, states: NDArray): - # states[:, self._circular_indices] = self._wrap_circular(states[:, self._circular_indices]) - # Apply boundary checks - # Apply model prediction to each chain - return -1 * self._interface.model_predict(states) - - def propose_step(self): - # Propose steps for all chains - proposed_states = self._matrix_handler.sample(self._n_chains)+self._chain_states - - # Calculate likelihoods for all proposed states (batch processing) - proposed_loglikelihoods = self._calc_likelihood(proposed_states) - - # Metropolis-Hastings acceptance step (vectorized) - valid_loglikelihoods = proposed_loglikelihoods > -1 * float(np.inf) - - # Compute acceptance probabilities - log_diff = proposed_loglikelihoods - self._current_loglikelihoods - acc_probs = np.minimum(1, np.exp(np.clip(log_diff, -100, 0))) # clip to avoid overflow - - # Generate uniform random numbers for comparison - u = np.random.uniform(0, 1, size=self._n_chains) - - # Update chains where proposed step is accepted - accept = (acc_probs > u) & valid_loglikelihoods - self._chain_states[accept] = proposed_states[accept] - self._current_loglikelihoods[accept] = proposed_loglikelihoods[accept] - - # Update covariance matrix only based on the first chain - self._matrix_handler.update(self._chain_states[0]) - - # Save the current state of all chains - self._dataset[self._current_step, :] = self._chain_states - self._current_step += 1 - - def save_mcmc_chain_to_pdf(self, filename: str, output_pdf: str): - # Open the HDF5 file - with h5py.File(filename, 'r') as f: - # Load the chain dataset - chain = f['chain'][:] - - _, n_params = chain.shape[1:] - - # Create a PdfPages object to save plots - with PdfPages(output_pdf) as pdf: - for i in tqdm(range(n_params)): - fig, ax = plt.subplots(figsize=(10, 6)) - - # Plot the chain for the i-th parameter - for j in range(self._n_chains): - ax.plot(chain[:, j, i], lw=0.5, label=f'Chain {j+1}') - ax.set_ylabel(self._interface.chain.plot_branches[i]) - ax.set_title(f"Parameter {self._interface.chain.plot_branches[i]} MCMC Chain") - ax.set_xlabel('Step') - # ax.legend() - - # Save the current figure to the PDF - pdf.savefig(fig) - plt.close(fig) # Close the figure to save memory - - print(f"MCMC chain plots saved to {output_pdf}") - - def __call__(self, n_steps, output_file_name: str): - print(f"Running MCMC for {n_steps} steps with {self._n_chains} chains") - - # Open the HDF5 file and create the dataset with the correct shape upfront - with h5py.File(output_file_name, 'w') as f: - self._dataset = f.create_dataset('chain', (n_steps, self._n_chains, self._n_dim), chunks=True) - - for _ in tqdm(range(n_steps)): - self.propose_step() - - self.save_mcmc_chain_to_pdf(output_file_name, "traces.pdf") diff --git a/src/MaCh3PythonUtils/fitters/multi_mcmc_gpu.py b/src/MaCh3PythonUtils/fitters/multi_mcmc_gpu.py index bd0dbd0..f4df8c7 100644 --- a/src/MaCh3PythonUtils/fitters/multi_mcmc_gpu.py +++ b/src/MaCh3PythonUtils/fitters/multi_mcmc_gpu.py @@ -28,7 +28,10 @@ def __init__(self, interface: TfInterface, n_chains: int = 1024, circular_params # boundaries self._upper_bounds = tf.convert_to_tensor(self._interface.scale_data(self._interface.chain.upper_bounds[:-1].reshape(1,-1)), dtype=tf.float32) self._lower_bounds = tf.convert_to_tensor(self._interface.scale_data(self._interface.chain.lower_bounds[:-1].reshape(1,-1)), dtype=tf.float32) - + + max_val = self._interface.chain.upper_bounds[-1] + + self._max_val_tf = tf.constant(max_val) self._circular_indices = self._get_circular_indices(circular_params) print(self._circular_indices) @@ -107,7 +110,10 @@ def apply_circular_bounds(idx): proposed_loglikelihoods = self._calc_likelihood(proposed_states) # Metropolis-Hastings acceptance step - log_diff = proposed_loglikelihoods - self._current_loglikelihoods + log_diff = proposed_loglikelihoods - self._current_loglikelihoods # Because network scales itself + + # Need to undo the scaling done in the NN + log_diff = self._max_val_tf * log_diff acc_probs = tf.minimum(1.0, tf.exp(tf.clip_by_value(log_diff, -100, 0))) # Generate uniform random values for comparison @@ -164,10 +170,10 @@ def _flush_async(self, final_flush=False): self._dataset[end_idx-len(steps_to_write):end_idx, :] = steps_to_write - def save_mcmc_chain_to_pdf(self, filename: str, output_pdf: str): + def save_mcmc_chain_to_pdf(self, filename: str, output_pdf: str, thinning: int=10, burnin: int=10000): # Open the HDF5 file and read the chain with h5py.File(filename, 'r') as f: - chain = f['chain'][:] + chain = f['chain'][burnin::thinning] # Need it to reflect the actual parameters in our fit so let's combine everything! rescaled_chain = [self._interface.invert_scaling(chain[1000:,i]) for i in range(self._n_chains)] diff --git a/src/MaCh3PythonUtils/machine_learning/file_ml_interface.py b/src/MaCh3PythonUtils/machine_learning/file_ml_interface.py index c67e36c..ffc5c1e 100644 --- a/src/MaCh3PythonUtils/machine_learning/file_ml_interface.py +++ b/src/MaCh3PythonUtils/machine_learning/file_ml_interface.py @@ -15,7 +15,7 @@ from scipy.optimize import minimize, OptimizeResult from sklearn import metrics -from sklearn.preprocessing import StandardScaler +from sklearn.preprocessing import StandardScaler, MinMaxScaler from sklearn.decomposition import PCA import matplotlib.pyplot as plt from matplotlib.backends.backend_pdf import PdfPages @@ -57,7 +57,11 @@ def __init__(self, chain: ChainHandler, prediction_variable: str, fit_name: str) # Scaling components self._scaler = StandardScaler() - self._pca_matrix = PCA() + # self._pca_matrix = PCA(n_components=0.95) + + self._label_scaler = MinMaxScaler(feature_range=(0, 1)) + + def __separate_dataframe(self)->Tuple[pd.DataFrame, pd.DataFrame]: """Split data frame into feature + label objects @@ -82,14 +86,19 @@ def set_training_test_set(self, test_size: float): self._training_data, self._test_data, self._training_labels, self._test_labels = train_test_split(features, labels, test_size=test_size) # Fit scaling pre-processors. These get applied properly when scale_data is called - scaled_training= self._scaler.fit_transform(self._training_data) - self._pca_matrix.fit(scaled_training) + _= self._scaler.fit_transform(self._training_data) + self._label_scaler.fit_transform(self._training_labels) + + # self._pca_matrix.fit(scaled_training) def scale_data(self, input_data): # Applies transformations to data set scale_data = self._scaler.transform(input_data) # scale_data = self._pca_matrix.transform(scale_data) return scale_data + + def scale_labels(self, labels): + return self._label_scaler.transform(labels) def invert_scaling(self, input_data): # Inverts transform @@ -119,7 +128,7 @@ def training_data(self)->pd.DataFrame: :rtype: pd.DataFrame """ if self._training_data is None: - return self._chain.ttree_array + return self._chain.ttree_array.iloc[:,:-1] return self._training_data @@ -131,7 +140,7 @@ def test_data(self)->pd.DataFrame: :rtype: pd.DataFrame """ if self._test_data is None: - return self._chain.ttree_array + return self._chain.ttree_array.iloc[:,:-1] return self._test_data @@ -204,14 +213,16 @@ def test_model(self): print("Training Results!") train_prediction = self.model_predict(self._training_data) - train_as_numpy = self._training_labels.to_numpy().T[0] + print(self._training_data) + print(train_prediction) + train_as_numpy = self.scale_labels(self._training_labels).T[0] self.evaluate_model(train_prediction, train_as_numpy, "train_qq_plot.pdf") print("=====\n\n") print("Testing Results!") test_prediction = self.model_predict(self._test_data) - test_as_numpy = self._test_labels.to_numpy().T[0] + test_as_numpy = self.scale_labels(self._test_labels).T[0] self.evaluate_model(test_prediction, test_as_numpy, outfile=f"{self._fit_name}") print("=====\n\n") @@ -262,7 +273,6 @@ def run_likelihood_scan(self, n_divisions: int = 500): plt.ylabel("-2*loglikelihood") pdf.savefig() plt.close() - @@ -276,6 +286,8 @@ def evaluate_model(self, predicted_values: Iterable, true_values: Iterable, outf :param outfile: File to output plots to, defaults to "" :type outfile: str, optional """ + + print(predicted_values) print(f"Mean Absolute Error : {metrics.mean_absolute_error(predicted_values,true_values)}") diff --git a/src/MaCh3PythonUtils/machine_learning/ml_factory.py b/src/MaCh3PythonUtils/machine_learning/ml_factory.py index 8867e1e..39a0d61 100644 --- a/src/MaCh3PythonUtils/machine_learning/ml_factory.py +++ b/src/MaCh3PythonUtils/machine_learning/ml_factory.py @@ -2,13 +2,19 @@ ML Factory implementation, effectively a selector for making models """ -from MaCh3PythonUtils.machine_learning.scikit_interface import SciKitInterface -from MaCh3PythonUtils.machine_learning.tf_interface import TfInterface -from MaCh3PythonUtils.machine_learning.normalizing_flow_interface import NormalisingFlowInterface +from MaCh3PythonUtils.machine_learning.scikit.scikit_interface import SciKitInterface +from MaCh3PythonUtils.machine_learning.tensorflow.tf_autotune_interface import TfAutotuneInterface +from MaCh3PythonUtils.machine_learning.tensorflow.tf_sequential_model import TfSequentialModel +from MaCh3PythonUtils.machine_learning.tensorflow.tf_residual_model import TfResidualModel +from MaCh3PythonUtils.machine_learning.tensorflow.tf_normalizing_flow_model import TfNormalizingFlowModel + +from MaCh3PythonUtils.machine_learning.tensorflow.tf_manual_interface import TfManualLayeredInterface +from MaCh3PythonUtils.machine_learning.tensorflow.tf_interface import TfInterface + +from MaCh3PythonUtils.file_handling.chain_handler import ChainHandler + import sklearn.ensemble as ske import tensorflow.keras as tfk -from MaCh3PythonUtils.file_handling.chain_handler import ChainHandler -import MaCh3PythonUtils.machine_learning.algorithms.normalizing_flow_structures as nfs class MLFactory: @@ -21,7 +27,10 @@ class MLFactory: "histboost" : ske.HistGradientBoostingRegressor }, "tensorflow": { - "sequential" : tfk.Sequential + "sequential" : TfSequentialModel, + "residual": TfResidualModel, + "normalizing_flow": TfNormalizingFlowModel, + "autotune": TfAutotuneInterface }, } @@ -77,32 +86,35 @@ def __make_scikit_model(self, algorithm: str, **kwargs)->SciKitInterface: interface = SciKitInterface(self._chain, self._prediction_variable, self._plot_name) interface.add_model(self.__setup_package_factory(package="scikit", algorithm=algorithm, **kwargs)) - return interface + return interface + - def __make_tensorflow_model(self, algorithm: str, **kwargs)->TfInterface: - """Generates TensorFlow model interface + def __make_tensorflow_layered_model(self, interface: TfManualLayeredInterface, layers: dict)->TfManualLayeredInterface: + for layer in layers: + layer_id = list(layer.keys())[0] + interface.add_layer(layer_id, layer[layer_id]) - :param algorithm: TensorFlow algorithm [NOT layers] - :type algorithm: str - :return: TfInterface wrapper around model - :rtype: _type_ - """ - interface = TfInterface(self._chain, self._prediction_variable, self._plot_name) + return interface + + def __make_tensorflow_model(self, algorithm: str, **kwargs)->TfInterface: + model_func = self.__IMPLEMENTED_ALGORITHMS["tensorflow"].get(algorithm.lower(), None) - interface.add_model(self.__setup_package_factory(package="tensorflow", algorithm=algorithm)) + if model_func is None: + raise Exception(f"Cannot find {algorithm}") - for layer in kwargs["Layers"]: - layer_id = list(layer.keys())[0] - - interface.add_layer(layer_id, layer[layer_id]) - - interface.build_model(kwargs["BuildSettings"]) + model: TfInterface = model_func(self._chain, self._prediction_variable, self._plot_name) + + # Ugh + if algorithm=="sequential" or algorithm=="residual": + print("HERE") + model = self.__make_tensorflow_layered_model(model, kwargs["Layers"]) + model.set_training_settings(kwargs.get("FitSettings")) + + + model.build_model(**kwargs["BuildSettings"]) - interface.set_training_settings(kwargs["FitSettings"]) + return model - return interface - - def make_interface(self, interface_type: str, algorithm: str, **kwargs): interface_type = interface_type.lower() match(interface_type): diff --git a/src/MaCh3PythonUtils/machine_learning/scikit_interface.py b/src/MaCh3PythonUtils/machine_learning/scikit/scikit_interface.py similarity index 100% rename from src/MaCh3PythonUtils/machine_learning/scikit_interface.py rename to src/MaCh3PythonUtils/machine_learning/scikit/scikit_interface.py diff --git a/src/MaCh3PythonUtils/diagnostics/__init__.py b/src/MaCh3PythonUtils/machine_learning/tensorflow/__init__.py similarity index 100% rename from src/MaCh3PythonUtils/diagnostics/__init__.py rename to src/MaCh3PythonUtils/machine_learning/tensorflow/__init__.py diff --git a/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_autotune_interface.py b/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_autotune_interface.py new file mode 100644 index 0000000..0bd65cf --- /dev/null +++ b/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_autotune_interface.py @@ -0,0 +1,104 @@ +from MaCh3PythonUtils.machine_learning.tf_interface import TfInterface +import tensorflow as tf +import pandas as pd +import numpy as np +import tensorflow.keras as tfk +import keras_tuner as kt + + +class TfAutotuneInterface(TfInterface): + # just make sure things don't break! + # gpus = tf.config.experimental.list_physical_devices('GPU') + # for gpu in gpus: + # tf.config.experimental.set_memory_growth(gpu, True) + + _epochs = 1000 + _val_split = 0.2 + + def build_model(self, **kwargs): + # Set up layers + n_layers = kwargs.get("n_layers", [5, 20, 1]) + + + activation_functions = kwargs.get("layer_activation", ["tanh", "relu"]) + + neurons_per_layer = kwargs.get("neurons_per_layer", [16, 2048, 100]) + learning_rate = kwargs.get("learning_rate", [1, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6]) + + regularization_rate = kwargs.get("regularization", [0.001, 0.01, 0.1, 1]) + + # Actually make the model + def model_builder(hp): + model = tfk.Sequential() + + hp_model_layers = hp.Int('n_layers', min_value = n_layers[0], + max_value=n_layers[1], + step=n_layers[2]) + + for i in range(hp_model_layers): + # Add layer + hp_layer_units = hp.Int(f'units_{i}', min_value = neurons_per_layer[0], + max_value=neurons_per_layer[1], + step=neurons_per_layer[2]) + hp_layer_activation = hp.Choice(f'activation_{i}', values=activation_functions) + + hp_layer_regularization =hp.Choice(f'regularization_{i}', regularization_rate) + + model.add(tfk.layers.Dense(units = hp_layer_units, + activation = hp_layer_activation, + kernel_regularizer = tfk.regularizers.L2(hp_layer_regularization))) + + # Add batch normalization + model.add(tfk.layers.BatchNormalization()) + + # Add output + model.add(tfk.layers.Dense(units=1, activation="linear")) + + hp_learning_rate = hp.Choice('learning_rate', learning_rate) + + model.compile(optimizer=tfk.optimizers.Adam(learning_rate=hp_learning_rate), + loss=kwargs.get('loss', 'mse'), + metrics=kwargs.get('metrics', ['mse', 'mae']) + ) + + return model + + + self._epochs = kwargs.get('epochs', 1000) + self._val_split = kwargs.get('validation_split', 0.2) + self._batch_size = kwargs.get('batch_size', 2048) + + # Hyperband parameters + hyperband_iterations = kwargs.get("hyperband_iterations", 100) + model_directory = kwargs.get("tuning_dir", "tuning-model") + project_name = kwargs.get("project_name", "tuning-project") + + self._model_tuner= kt.Hyperband(model_builder, + objective='val_loss', + max_epochs=self._epochs, + hyperband_iterations=hyperband_iterations, + directory=model_directory, + # distribution_strategy=tf.distribute.MirroredStrategy(), + project_name=project_name, + overwrite=False,) + + def train_model(self): + + scaled_data = self.scale_data(self._training_data) + scaled_labels = self.scale_labels(self._training_labels) + + stop_early = tfk.callbacks.EarlyStopping(monitor='val_loss', patience=5) + lr_schedule = tfk.callbacks.ReduceLROnPlateau(monitor="loss", patience=5, factor=0.5, min_lr=1e-6, verbose=1) + + self._model_tuner.search(scaled_data, scaled_labels, + epochs= self._epochs, + validation_split=self._val_split, + batch_size=self._batch_size, + callbacks=[stop_early]) + + best_hps=self._model_tuner.get_best_hyperparameters()[0] + print("Finished auto-tuning") + + self._model = self._model_tuner.hypermodel.build(best_hps) + + self._model.fit(scaled_data, scaled_labels, epochs=self._epochs, batch_size=self._batch_size , validation_split=0.2, callbacks=[stop_early, lr_schedule]) \ No newline at end of file diff --git a/src/MaCh3PythonUtils/machine_learning/tf_interface.py b/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_interface.py similarity index 70% rename from src/MaCh3PythonUtils/machine_learning/tf_interface.py rename to src/MaCh3PythonUtils/machine_learning/tensorflow/tf_interface.py index 3b60036..8ef6a9a 100644 --- a/src/MaCh3PythonUtils/machine_learning/tf_interface.py +++ b/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_interface.py @@ -4,6 +4,9 @@ import pandas as pd import numpy as np import tensorflow.keras as tfk +import tensorflow_probability as tfp + + class TfInterface(FileMLInterface): __TF_LAYER_IMPLEMENTATIONS = { @@ -12,15 +15,9 @@ class TfInterface(FileMLInterface): "batchnorm": tfk.layers.BatchNormalization } - __TF_REGULARIZERS = { - "l2" : tf.keras.regularizers.L2 - } - _layers = [] _training_settings = {} - - def add_layer(self, layer_id: str, layer_args: dict): """Add new layer to TF model @@ -33,35 +30,16 @@ def add_layer(self, layer_id: str, layer_args: dict): if layer_id not in self.__TF_LAYER_IMPLEMENTATIONS.keys(): raise ValueError(f"{layer_id} not implemented yet!") - if "kernel_regularizer" in layer_args.keys(): + if layer_args.get("kernel_regularizer"): # Hacky, swaps string value of regularliser for proper one - reg = layer_args["kernel_regularizer"] - reg_name = list(reg.keys())[0] - layer_args["kernel_regularizer"] = self.__TF_REGULARIZERS[reg_name.lower()](reg[reg_name]) + layer_args["kernel_regularizer"] = tfk.regularizers.L2(0.2) self._layers.append(self.__TF_LAYER_IMPLEMENTATIONS[layer_id.lower()](**layer_args)) - def build_model(self, model_args: dict): - """Build and compile TF model - - :param model_args: Model arguments as dictionary - :type model_args: dict - :raises ValueError: Model not set up yet - """ - if self._model is None or not self._layers: - raise ValueError("No model can be built! Please setup model and layers") - - for layer in self._layers: - self._model.add(layer) - - self._model.build() - optimizer = tfk.optimizers.AdamW(learning_rate=model_args.get("learning_rate", 1e-3), - weight_decay=1e-4, clipnorm=1.0) - - self._model.compile(**model_args, optimizer=optimizer) - - + def build_model(self, _: dict): + return None + def set_training_settings(self, kwargs): """Set training settings, needs to be done early for...reasons @@ -78,7 +56,8 @@ def train_model(self): lr_schedule = tfk.callbacks.ReduceLROnPlateau(monitor="loss", patience=5, factor=0.5, min_lr=1e-6, verbose=1) self._model.fit(scaled_data, self._training_labels, **self._training_settings, callbacks=[lr_schedule]) - + print(f"Using loss function: {self._model.loss}") + def save_model(self, output_file: str): """Save model to file @@ -95,8 +74,8 @@ def load_model(self, input_file: str): :param input_file: Name offile to load model from :type input_file: str - """ - + """ + print(f"Loading model from {input_file}") self._model = tf.keras.models.load_model(input_file) @@ -114,10 +93,11 @@ def model_predict(self, test_data: pd.DataFrame): if self._model is None: return np.zeros(len(test_data)) - return self._model.predict(scaled_data, verbose=False).T[0] - + pred = self._model.predict(scaled_data, verbose=False) + print(f"PREDICTION: {pred}") + # return self._model.predict(scaled_data, verbose=False).T[0] + return pred def model_predict_no_scale(self, test_data): # Same as above but specifically for TF, optimised to avoid if statement... - return self._model(test_data, training=False) - \ No newline at end of file + return self._model(test_data, training=False) \ No newline at end of file diff --git a/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_manual_interface.py b/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_manual_interface.py new file mode 100644 index 0000000..4f05408 --- /dev/null +++ b/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_manual_interface.py @@ -0,0 +1,64 @@ +from MaCh3PythonUtils.machine_learning.tf_interface import TfInterface +import tensorflow as tf +import pandas as pd +import numpy as np +import tensorflow.keras as tfk +import keras_tuner as kt + +class TfManualInterface(TfInterface): + _training_settings = {} + + def set_training_settings(self, kwargs): + """Set training settings, needs to be done early for...reasons + + :param kwargs: training kwargs + :type kwargs: kwargs + """ + self._training_settings = kwargs + + def train_model(self): + """train model + """ + scaled_data = self.scale_data(self._training_data) + scaled_labels = self.scale_labels(self._training_labels) + + lr_schedule = tfk.callbacks.ReduceLROnPlateau(monitor="val_loss", patience=10, factor=0.5, min_lr=1e-8, verbose=1) + stop_early = tfk.callbacks.EarlyStopping(monitor='val_loss', patience=20) + + self._model.fit(scaled_data, scaled_labels, **self._training_settings, callbacks=[lr_schedule, stop_early]) + + +class TfManualLayeredInterface(TfManualInterface): + __TF_LAYER_IMPLEMENTATIONS = { + "dense": tfk.layers.Dense, + "dropout": tfk.layers.Dropout, + "batchnorm": tfk.layers.BatchNormalization + } + + _layers = [] + + def add_layer(self, layer_id: str, layer_args: dict): + """Add new layer to TF model + + :param layer_id: Layer type [dense/dropout] + :type layer_id: str + :param layer_args: kwargs for layer + :type layer_args: dict + :raises ValueError: Layer type not implemented in __TF_LAYER_IMPLEMENTATIONS yet + """ + if layer_id not in self.__TF_LAYER_IMPLEMENTATIONS.keys(): + raise ValueError(f"{layer_id} not implemented yet!") + + if layer_args.get("kernel_regularizer", False): + # Hacky, swaps string value of regularliser for proper one + layer_args["kernel_regularizer"] = tfk.regularizers.L2(layer_args["kernel_regularizer"]) + + self._layers.append(self.__TF_LAYER_IMPLEMENTATIONS[layer_id.lower()](**layer_args)) + + + + + + + + diff --git a/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_normalizing_flow_model.py b/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_normalizing_flow_model.py new file mode 100644 index 0000000..e40577b --- /dev/null +++ b/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_normalizing_flow_model.py @@ -0,0 +1,94 @@ +from MaCh3PythonUtils.machine_learning.tf_manual_interface import TfManualInterface + +import tensorflow_probability as tfp +import tensorflow.keras as tfk +import tensorflow as tf +from typing import List +from matplotlib.backends.backend_pdf import PdfPages +import matplotlib.pyplot as plt +from tqdm import tqdm +import seaborn as sns + +tfd = tfp.distributions +tfb = tfp.bijectors + + +# class FlowLayer(tfk.layers.Layer): +# def __init__(self, dist): +# super().__init__() +# self._layer_dist = dist.log_prob + +# def call(self, x): +# return self._layer_dist(x) + + +class NormalizingFlow: + def __init__(self, hidden_units: List[int], n_bijectors: int, input_dim: int): + self._hidden_units = hidden_units + self._n_bijectors = n_bijectors + self._input_dim=input_dim + + # Create an autoregressive bijector + def _create_autoregressive_bijector(self): + return tfb.MaskedAutoregressiveFlow( + shift_and_log_scale_fn=tfb.AutoregressiveNetwork( + params=2, # Shift and log scale + hidden_units=self._hidden_units, + input_shape=(None,self._input_dim), + ) + ) + + def _create_normalizing_flow(self): + bijectors = [self._create_autoregressive_bijector() for _ in range(self._n_bijectors)] + # Add a final bijector for numerical stability (optional) + bijectors.append(tfb.Permute(permutation=list(range(self._input_dim))[::-1])) + return tfb.Chain(bijectors) + + def __call__(self): + flow_bijector = self._create_normalizing_flow() + base_distribution = tfd.MultivariateNormalDiag(loc=tf.zeros(self._input_dim), scale_diag=tf.ones(self._input_dim)) + # Grab the distribution + return tfd.TransformedDistribution(distribution=base_distribution, bijector=flow_bijector) + + + + +class TfNormalizingFlowModel(TfManualInterface): + def build_model(self, model_args): + input_dim = self.chain.ndim-1 + self._model = NormalizingFlow(model_args.get("hidden_units", [100]), model_args.get("n_bijectors", 1), input_dim)() + self._optimizer = tfk.optimizers.Adam(model_args.get("learning_rate", 1e-3)) + + def nll_loss(self, features): + return -tf.reduce_mean(self._model.log_prob(features)) + + + def train_model(self): + epochs = self._training_settings["epochs"] + batch_size = self._training_settings["batch_size"] + + scaled_data = tf.data.Dataset.from_tensor_slices(self.scale_data(self._training_data)).batch(batch_size) + + for epoch in tqdm(range(epochs)): + epoch_loss = 0 + for batch in scaled_data: + with tf.GradientTape() as tape: + loss = self.nll_loss(batch) + grads = tape.gradient(loss, self._model.trainable_variables) + self._optimizer.apply_gradients(zip(grads, self._model.trainable_variables)) + epoch_loss += loss.numpy() + print(f"Epoch {epoch+1}, Loss: {epoch_loss / len(scaled_data)}") + + def test_model(self): + print("Testing") + surrogate_samples = self._model.sample(10000).numpy() + + with PdfPages("likelihood_free_inference.pdf") as pdf: + for i in range(self._chain.ndim-1): + sns.histplot(surrogate_samples[:, i], label="Surrogate", color="blue", fill=False) + sns.histplot(self._chain.ttree_array[:, i], label="Test Data", color="blue", fill=False) + plt.legend() + plt.xlabel(f"{self._chain.plot_branches[i]}") + plt.ylabel("Density") + pdf.savefig() + plt.close() diff --git a/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_residual_model.py b/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_residual_model.py new file mode 100644 index 0000000..2c235c9 --- /dev/null +++ b/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_residual_model.py @@ -0,0 +1,31 @@ +from MaCh3PythonUtils.machine_learning.tensorflow.tf_manual_interface import TfManualLayeredInterface +import tensorflow.keras as tfk + + +class TfResidualModel(TfManualLayeredInterface): + def build_model(self, model_args: dict): + input_shape = self.training_data.shape[1:] # Assuming shape is (batch_size, features) + network_input = tfk.layers.Input(shape=input_shape) + + # Initial layer + x = self._layers[0](network_input) + + # Apply each layer in self._layers, adding residual connections every other layer + for i, layer in enumerate(self._layers[1:], start=1): + # Apply layer + x_new = layer(x) + + # Add residual connection every other layer if shapes match + if i % 2 == 0 and x.shape == x_new.shape: + x = tfk.layers.add([x, x_new]) + else: + x = x_new # Update x with current layer output + + # Define and compile the model + self._model = tfk.Model(inputs=network_input, outputs=x) + optimizer = tfk.optimizers.AdamW(learning_rate=model_args.get("learning_rate", 1e-5), + weight_decay=1e-4, clipnorm=1.0) + + _ = model_args.pop("learning_rate", None) + + self._model.compile(optimizer=optimizer, **model_args) \ No newline at end of file diff --git a/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_sequential_model.py b/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_sequential_model.py new file mode 100644 index 0000000..4a88408 --- /dev/null +++ b/src/MaCh3PythonUtils/machine_learning/tensorflow/tf_sequential_model.py @@ -0,0 +1,29 @@ +from MaCh3PythonUtils.machine_learning.tensorflow.tf_manual_interface import TfManualLayeredInterface + +import tensorflow.keras as tfk + +class TfSequentialModel(TfManualLayeredInterface): + + def build_model(self, model_args: dict): + """Build and compile TF model + + :param model_args: Model arguments as dictionary + :type model_args: dict + :raises ValueError: Model not set up yet + """ + self._model = tfk.Sequential() + + if self._model is None or not self._layers: + raise ValueError("No model can be built! Please setup model and layers") + + for layer in self._layers: + self._model.add(layer) + + self._model.build() + optimizer = tfk.optimizers.AdamW(learning_rate=model_args.get("learning_rate", 1e-5), + weight_decay=1e-4, clipnorm=1.0) + + model_args.pop("learning_rate", None) + + + self._model.compile(**model_args, optimizer=optimizer) \ No newline at end of file