From bf254d729953f8954035f178c8541614ed75921e Mon Sep 17 00:00:00 2001 From: ashriva16 Date: Mon, 12 May 2025 16:27:47 -0400 Subject: [PATCH 1/8] temporary working andie over nomad branch --- nohup.out | 10 + pdm.lock | 12 +- pyproject.toml | 1 + scripts/andie_mock.py | 331 ++++++++++++++++++ src/dial_dataclass/dial_dataclass.py | 6 +- src/dial_service/backends/__init__.py | 1 + src/dial_service/backends/andie_backend.py | 329 +++++++++++++++++ src/dial_service/core.py | 5 + .../service_specific_dataclasses.py | 5 +- 9 files changed, 695 insertions(+), 5 deletions(-) create mode 100644 nohup.out create mode 100644 scripts/andie_mock.py create mode 100644 src/dial_service/backends/andie_backend.py diff --git a/nohup.out b/nohup.out new file mode 100644 index 0000000..1be0a31 --- /dev/null +++ b/nohup.out @@ -0,0 +1,10 @@ +train -------------------- +initial chisq 135694.5(47) constraints=9.94437 +# steps: 1, # draws: 40 +final chisq 0.0(47) constraints=13.839 +BK : 75(24) +J : 0.5(13) +M0 : 17.0(87) +TN : 0.005(38)e3 +predict -------------------- +Predicting diff --git a/pdm.lock b/pdm.lock index 0d536f7..f1cc828 100644 --- a/pdm.lock +++ b/pdm.lock @@ -5,7 +5,7 @@ groups = ["default", "dev", "lint", "test"] strategy = ["inherit_metadata"] lock_version = "4.5.0" -content_hash = "sha256:1cca3fe4ae7dd4c5871e031ddd5babaf11b073f2db7f050d34c5794acec1e85c" +content_hash = "sha256:ffa992ef9dce1e6c6e11033dac52e49143fc3a63aba9ef644839c3713eee3164" [[metadata.targets]] requires_python = ">=3.10" @@ -115,6 +115,16 @@ files = [ {file = "bson-0.5.10.tar.gz", hash = "sha256:d6511b2ab051139a9123c184de1a04227262173ad593429d21e443d6462d6590"}, ] +[[package]] +name = "bumps" +version = "0.9.3" +summary = "Data fitting with bayesian uncertainty analysis" +groups = ["dev"] +files = [ + {file = "bumps-0.9.3-py3-none-any.whl", hash = "sha256:269f5e86c50fde2b80015b8c9aff9b06eb045be675bdef065d2897670c5b56e5"}, + {file = "bumps-0.9.3.tar.gz", hash = "sha256:3295298f7fe1b2392bb2ff88c7a0ae69d77a7698af580a81570052ca05ba530f"}, +] + [[package]] name = "certifi" version = "2024.8.30" diff --git a/pyproject.toml b/pyproject.toml index ffa1483..c9f0c99 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -135,4 +135,5 @@ build-backend = "setuptools.build_meta" [dependency-groups] dev = [ "statsmodels>=0.14.4", + "bumps>=0.9.3", ] diff --git a/scripts/andie_mock.py b/scripts/andie_mock.py new file mode 100644 index 0000000..066e701 --- /dev/null +++ b/scripts/andie_mock.py @@ -0,0 +1,331 @@ +import argparse +import json +import logging +import os +import sys +from pathlib import Path + +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np + +from intersect_sdk import ( + INTERSECT_JSON_VALUE, + HierarchyConfig, + IntersectClient, + IntersectClientCallback, + IntersectClientConfig, + IntersectDirectMessageParams, + default_intersect_lifecycle_loop, +) + +from dial_dataclass import ( + DialInputPredictions, + DialInputSingleOtherStrategy, + DialWorkflowCreationParamsClient, + DialWorkflowDatasetUpdate, +) + +mpl.use('agg') + +logging.basicConfig(level=logging.INFO) +logger = logging.getLogger(__name__) + +def peak_val_at_T(T: np.ndarray) -> np.ndarray: + result = 40.0 * 1.0/(1.0+np.exp((T-50)/5)) + 35.0 + logger.debug(result) + return result + + + +class IntersectCallbackError(Exception): + def __init__(self, operation, payload): + message = f"Intersect callback error during operation '{operation}'. Payload: {payload}" + super().__init__(message) + +class IntersectCallbackEnd(Exception): + def __init__(self): + message = "Stopping Intersect Calls" + super().__init__(message) + +class ActiveLearningOrchestrator: + + def __init__(self, service_destination: str): + #Maximum Measurments in Temperature: + Temperature_Loops = 30 + + T_start = 5.0 #Kelvin + T_stop = 300.0 #Kelvin + T_step = .5 #Kelvin + T_grid = np.linspace(T_start, T_stop, int((T_stop-T_start)/T_step)+1).reshape(-1,1) + + self.T_step = T_step + above_TN = 0 + + + self.above_TN = above_TN + self.last_idx = 0 + + self.bounds = np.array([[T_start, T_stop]]) + self.num_dims = len(self.bounds) + + self.x_raw = np.array([[ T_start]]) + self.x_test = np.array(T_grid) + self.y_raw = peak_val_at_T(self.x_raw) + + self.meshgrid_size = len(T_grid) + # self.grid_points = [ + # np.linspace(dim_bounds[0], dim_bounds[1], self.meshgrid_size) + # for dim_bounds in self.bounds + # ] + # self.meshgrids = np.meshgrid(*self.grid_points, indexing='ij') + # self.x_grid = np.hstack([mg.reshape(-1, 1) for mg in self.meshgrids]) + + # Active learning variables + self.dataset_x = self.x_raw.reshape(-1, 1).tolist() + print(self.dataset_x) + self.dataset_y = self.y_raw.reshape(-1).tolist() + self.test_points = self.x_test.reshape(-1, 1).tolist() + + self.kernel = None + self.kernel_args = None + self.backend = 'andie' + self.backend_args = None + self.strategy = None + self.strategy_args = None + self.niter = 0 + self.max_iter = Temperature_Loops + self.at_grids = False + self.variance_grid = None + self.mean_grid = None + self.variance_test = None + self.mean_test = None + self.x_next = None + + # Intersect variables + self.workflow_id = None + self.service_destination = service_destination + self.initialize_workflow_message = IntersectClientCallback(messages_to_send=[ + IntersectDirectMessageParams( + destination=self.service_destination, + operation='dial.initialize_workflow', + payload=DialWorkflowCreationParamsClient( + dataset_x=self.dataset_x, + dataset_y=self.dataset_y, + bounds=self.bounds, + kernel=self.kernel, + backend=self.backend, + preprocess_standardize=False, + y_is_good=True, + seed=20)) + ]) + + def callback(self, operation: str, **kwargs) -> IntersectClientCallback: + print("send", operation) + + next_payload = None + self.at_grids = kwargs.get('at_grids', True) + + if operation == 'dial.get_surrogate_values': + + _points_to_predict = np.array(self.test_points).reshape(-1, self.num_dims) + + next_payload = DialInputPredictions( + workflow_id=self.workflow_id, + points_to_predict=_points_to_predict, + kernel_args=self.kernel_args, + extra_args={'above_TN': self.above_TN, + 'last_idx': self.last_idx}, + ) + + elif operation == 'dial.get_next_point': + next_payload = DialInputSingleOtherStrategy( + workflow_id=self.workflow_id, + strategy=self.strategy, + strategy_args=self.strategy_args, + bounds=self.bounds.tolist(), + kernel_args=self.kernel_args, + extra_args={'above_TN': self.above_TN, + 'last_idx': self.last_idx, + 'T_step': self.T_step}, + ) + + elif operation == 'dial.update_workflow_with_data': + next_payload = DialWorkflowDatasetUpdate( + workflow_id=self.workflow_id, + next_x=self.dataset_x[-1], + next_y=self.dataset_y[-1], + ) + + else: + err_msg = f'Unknown operation received: {operation}' + raise Exception(err_msg) # noqa: TRY002 (INTERSECT interaction mechanism) + + return IntersectClientCallback(messages_to_send=[ + IntersectDirectMessageParams( + destination=self.service_destination, + operation=operation, + payload=next_payload) + ]) + + def __call__(self, _source: str, operation: str, _has_error: bool, + payload: INTERSECT_JSON_VALUE) -> IntersectClientCallback: + print(operation) + if _has_error: + print('============ERROR==============', file=sys.stderr) + print(operation, file=sys.stderr) + print(payload, file=sys.stderr) + raise IntersectCallbackError(operation, payload) + + if operation == 'dial.initialize_workflow': + self.workflow_id = payload + print("\n","--"*20,"\n") + return self.callback('dial.get_surrogate_values') + # TODO: repplace this with update_workflow_with_data that can also train the model + + # ----------------- Active learning loop ----------------- + if operation == 'dial.get_surrogate_values': + self.handle_surrogate_values(payload) + + if self.at_grids: + print(f"Step {self.niter}") + return self.callback('dial.get_surrogate_values', at_grids=False) + else: + return self.callback('dial.get_next_point') + + elif operation == 'dial.get_next_point': + self.handle_next_points(payload) + return self.callback('dial.update_workflow_with_data') + + elif operation == 'dial.update_workflow_with_data': + self.niter += 1 + return self.callback('dial.get_surrogate_values') + + else: + raise IntersectCallbackError(operation, payload) + + def handle_surrogate_values(self, payload): + + self.variance_grid = np.array(payload[1]).reshape(\ + (self.meshgrid_size,) * self.num_dims) + self.mean_grid = np.array(payload[0]).reshape(\ + (self.meshgrid_size,) * self.num_dims) + + # end of active learning loop after max_iter + if self.niter > self.max_iter: + raise IntersectCallbackEnd() + + def handle_next_points(self, payload): + + print(payload) + self.x_next = [payload[0]] + self.above_TN = payload[1] + self.last_idx = payload[2] + + print(f'Running simulation at ({self.x_next}): ', end='', flush=True) + + y = peak_val_at_T(*self.x_next) + print(f'{y:.3f}') + print(f'Adding ({self.x_next}, {y}) to dataset') + print(f'Current dataset: {self.dataset_x}') + print(f'Current dataset: {self.dataset_y}') + + self.dataset_x.append(self.x_next) + self.dataset_y.append(y) + + print(f'new dataset: {self.dataset_x}') + print(f'new dataset: {self.dataset_y}') + + optpos = np.argmax(self.dataset_y) + y_opt = self.dataset_y[optpos] + optimal_coords = self.dataset_x[optpos] + coord_str = ', '.join([f'{coord:.2f}' for coord in optimal_coords]) + print(f'Optimal simulated datapoint at ({coord_str}), y={y_opt:.3f}\n') + + def graph(self): + + plt.clf() + + fig, axs = plt.subplots(2, 1, figsize=(10, 8), sharex=True) + + # First subplot: Mean and variance with training data + axs[0].plot(self.x_grid, self.mean_grid, label='Mean Prediction') + axs[0].fill_between( + self.x_grid[:, 0], + self.mean_grid + 2 * self.variance_grid, + self.mean_grid - 2 * self.variance_grid, + alpha=0.5, + label='Confidence Interval' + ) + axs[0].scatter( + np.array(self.dataset_x)[:-1, 0], + np.array(self.dataset_y)[:-1], + color='black', + marker='o', + label='Training Data' + ) + if self.x_next is not None: + axs[0].axvline( + x=self.x_next[0], + color='red', + linestyle='--') + axs[0].set_ylabel('Response, y') + axs[0].legend() + axs[0].grid(True) + + # Second subplot: Acquisition function + if self.strategy_args is not None: + if self.mean_grid is not None and self.variance_grid is not None: + exploit = self.strategy_args.get('exploit', 0.0) + explore = self.strategy_args.get('explore', 1.0) + acquisition_values = exploit * self.mean_grid + explore * np.sqrt(self.variance_grid) + else: + acquisition_values = np.zeros_like(self.x_grid) + + axs[1].plot(self.x_grid, acquisition_values) + if self.x_next is not None: + axs[1].axvline( + x=self.x_next[0], + color='red', + linestyle='--', + label='Next Point (x_next)' + ) + axs[1].set_xlabel('Features, x') + axs[1].set_ylabel('Acquisition Value') + axs[1].legend() + axs[1].grid(True) + + plt.tight_layout() + plt.savefig('graph.png') + plt.close(fig) + + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description='Automated client') + parser.add_argument( + '--config', + type=Path, + default=os.environ.get('DIAL_CONFIG_FILE', Path(__file__).parents[1] / 'remote-conf.json'), + ) + args = parser.parse_args() + + try: + with Path(args.config).open('rb') as f: + from_config_file = json.load(f) + except (json.decoder.JSONDecodeError, OSError) as e: + logger.critical('unable to load config file: %s', str(e)) + sys.exit(1) + + active_learning = ActiveLearningOrchestrator( + service_destination=HierarchyConfig( + **from_config_file['intersect-hierarchy'] + ).hierarchy_string('.') + ) + + config = IntersectClientConfig( + initial_message_event_config=active_learning.initialize_workflow_message, + **from_config_file['intersect'], + ) + + client = IntersectClient(config=config, user_callback=active_learning) + default_intersect_lifecycle_loop(client) diff --git a/src/dial_dataclass/dial_dataclass.py b/src/dial_dataclass/dial_dataclass.py index bd50ef7..86a5038 100644 --- a/src/dial_dataclass/dial_dataclass.py +++ b/src/dial_dataclass/dial_dataclass.py @@ -6,7 +6,7 @@ PositiveIntType = Annotated[int, Field(ge=0)] -_POSSIBLE_BACKENDS = ('sklearn', 'gpax') +_POSSIBLE_BACKENDS = ('sklearn', 'gpax', 'andie') class _DialWorkflowCreationParams(BaseModel): @@ -37,7 +37,7 @@ class _DialWorkflowCreationParams(BaseModel): description='If true, treat higher y values as better (e.g. y represents yield or profit). If false, opposite (e.g. y represents error or waste)' ), ] - kernel: Literal['rbf', 'matern'] + kernel: Literal['rbf', 'matern', None] bounds: list[ Annotated[ Annotated[list[float], Field(min_length=2, max_length=2)], @@ -115,7 +115,7 @@ class DialInputSingleConfidenceBound(BaseModel): class DialInputSingleOtherStrategy(BaseModel): workflow_id: ValidatedObjectId - strategy: Literal['random', 'uncertainty', 'expected_improvement', 'upper_confidence_bound'] + strategy: Literal['random', 'uncertainty', 'expected_improvement', 'upper_confidence_bound', None] strategy_args: dict[str, Union[float, int, bool]] | None = Field(default=None) kernel_args: dict[str, Union[float, int, bool, str, list[float], tuple]] | None = Field(default=None) backend_args: dict[str, Union[float, int, bool, str, list[float], tuple]] | None = Field(default=None) diff --git a/src/dial_service/backends/__init__.py b/src/dial_service/backends/__init__.py index 04f9459..02889ed 100644 --- a/src/dial_service/backends/__init__.py +++ b/src/dial_service/backends/__init__.py @@ -8,6 +8,7 @@ _BACKENDS = { 'gpax': ('dial_service.backends.gpax_backend', 'GpaxBackend'), 'sklearn': ('dial_service.backends.sklearn_backend', 'SklearnBackend'), + 'andie': ('dial_service.backends.andie_backend', 'AndieBackend'), } _MODEL = TypeVar('_MODEL') diff --git a/src/dial_service/backends/andie_backend.py b/src/dial_service/backends/andie_backend.py new file mode 100644 index 0000000..cca93c1 --- /dev/null +++ b/src/dial_service/backends/andie_backend.py @@ -0,0 +1,329 @@ +"""NOTE: This file should not be imported in application code except dynamically via the get_backend_module function in __init__.py .""" + +import numpy as np +import scipy as sp +from scipy.optimize import root, OptimizeResult + +from bumps.names import Curve, FitProblem, fit +from bumps.formatnum import format_uncertainty +from bumps.bounds import BoundedNormal + + +from . import AbstractBackend + + +class AndieBackend( + AbstractBackend[OptimizeResult, None, tuple[np.ndarray, np.ndarray]] +): + +### set correct model type in first argument + + ### --- These are hyperparameters --- ### + ### --- Define the priors on the Isothermal parameters --- ### + + #Transition Temperature + TN_guess = 60.0 + TN_std = 20 + TN_limits = (0.01, 200.0) + + #Background + BK_guess = 70.0 + BK_std = 20.0 + BK_limits = (0.0, 150.0) + + #Second Order Scale + M0_guess = 18.0 + M0_std = 3.0 + M0_limits = (0.01, 35.0) + + #Total Anglular Momentum + J_guess = 0.6 + J_std = 0.5 + J_limits = (0.01, 14.0) + + + #Define the prior distributions + TN_dist = BoundedNormal(TN_guess, TN_std, limits=TN_limits) + M0_dist = BoundedNormal(M0_guess, M0_std, limits=M0_limits) + J_dist = BoundedNormal(J_guess, J_std, limits=J_limits) + BK_dist = BoundedNormal(BK_guess, BK_std, limits=BK_limits) + + ### --- End of hyperparameters --- ### + + ### --- The following functions can be put into a model.py file --- ### + ### --- Brillouin_J(), Wiess(), Mag_solve(), and second_order_I_vs_T() --- ### + + # Define the Brillouin Function + @staticmethod + def _Brillouin_J(x, J): + # x is a matrix + # J is a row-vector + + f1 = (2*J+1)/(2*J) #factor 1 + f2 = 1/(2*J) #factor 2 + B_J = f1/(np.tanh(f1*x)) - f2/(np.tanh(f2*x)) + return B_J + + # Define the Wiess equation (re-arranged so that it always = 0) + @staticmethod + def _Wiess(t, m, J): + # t is a collum vector + # J is a row vector + # m is a matrix + + f1 = 3.0*J/(J +1.0) # factor 1 + f2 = 1.0/t #factor 2 + + m = np.ones_like(t)*m + x = f1*m*f2 + + + W = m - AndieBackend._Brillouin_J(x, J) + return W.reshape(-1) + + #Define a Magnetization Solver for givin t and J: + @staticmethod + def _Mag_solve(t, J): + # t is a collum vector + # J is a row vector + + def fun(x): + x= x.reshape((t.shape[0], J.shape[1])) + return AndieBackend._Wiess(t, x, J) + + X0 = np.ones((t.shape[0], J.shape[1]), dtype=np.float64) # matrix of every t,J combination + + solution = root(fun, X0, method ='hybr') + + m = np.array(solution.x).reshape((t.shape[0], J.shape[1])) + return m + + #Define the Intensity as function of Temperature for Second order phase transitions + @staticmethod + def _second_order_I_vs_T(Temperatures, TN, M0, J, BK): + # TN, M0, J, BK are row vectors + TN = np.array(TN).reshape(1,-1) + M0 = np.array(M0).reshape(1,-1) + J = np.array(J).reshape(1,-1) + BK = np.array(BK).reshape(1,-1) + # Temperatures is a collumn vector + Temperatures = Temperatures.reshape(-1,1) + + Inverse_TN = 1/TN + + t = Temperatures @ Inverse_TN #Create a matrix of every (temperature, TN) combination + + #Find where all the temperatures below the TN. + t_low_truth = t <1.0 + t_low_truth = t_low_truth.reshape(t.shape) + t_low = t*t_low_truth + + #Because "Mag_solve" is expensive, only pass the non-zero rows of the t_low matrix. + t_low_cut = t_low[t_low.sum(axis=1) !=0] + + #Calculate the reduced magnetization for below Tn\ + m_low = AndieBackend._Mag_solve(t_low_cut, J) #returns a matrix where each collumn of m is associated with a collum of J + + #The mag is zero above the TN, construct that matrix + m_high = np.zeros((t.shape[0]-t_low_cut.shape[0], t.shape[1])) + #Combine the m matrixies for below and above the TN + m = np.concatenate((m_low,m_high), axis=0) + + #"t_low_cut" could have a raggid bottom edge, so only keep the m values below the TN + m = m*t_low_truth + + #Calcualte the Magnetization (not reduced) + M = m * M0 #row by row mulitplication of m by M0 + + ##Calculate the intensities + #Square of Magnetization + I = M**2 + + #Add the background + I = I + BK + + return I + + @staticmethod + def _Second_order_model(X_grid, Y_grid): + M = Curve(AndieBackend._second_order_I_vs_T, X_grid, Y_grid, + dy = np.sqrt(Y_grid/100), + TN = AndieBackend.TN_guess, + M0 = AndieBackend.M0_guess, + J = AndieBackend.J_guess, + BK = AndieBackend.BK_guess) + + M.TN.bounds = AndieBackend.TN_dist + M.M0.bounds = AndieBackend.M0_dist + M.J.bounds = AndieBackend.J_dist + M.BK.bounds = AndieBackend.BK_dist + + if len(X_grid) < 5 : + problem = FitProblem(M, partial = True) + else: + problem = FitProblem(M, partial = False) + + return problem + + @staticmethod + def train_model(data): + print("train", "-"*20) + + # X_grid = [temperature for temperature in list(np.array(data.X_train)[0,:])] + X_grid = np.array(data.X_train) + Y_grid = np.array(data.Y_train) + # X_grid = data.measuredTemperatures + # Y_grid = data.measuredTemperatureDependentIntensities + problem = AndieBackend._Second_order_model(X_grid, Y_grid) + + method = 'dream' + + print("initial chisq", problem.chisq_str()) + result = fit(problem, method=method, xtol=1e-6, ftol=1e-8, + samples = 10, + burn = 100, + pop = 10, + init = 'eps', + thin = 1, + alpha = 0.1, + outliers = 'none', + trim = False, + steps = 0) + print("final chisq", problem.chisq_str()) + for k, v, dv in zip(problem.labels(), result.x, result.dx): + print(k, ":", format_uncertainty(v, dv)) + + return result + + @staticmethod + def predict(posterior_results, data): + print("predict", "-"*20) + + #Extract the Posterior of the parameters + draw = posterior_results.state.draw() + + #Posterior results are in alphabetical order + Posterior_TN = draw.points[:,3] + Posterior_M0 = draw.points[:,2] + Posterior_J = draw.points[:,1] + Posterior_BK = draw.points[:,0] + + # #Calcualte the mean of the parameters + # TN_mean = np.mean(Posterior_TN).reshape(1,-1) + # M0_mean = np.mean(Posterior_M0).reshape(1,-1) + # J_mean = np.mean(Posterior_J).reshape(1,-1) + # BK_mean = np.mean(Posterior_BK).reshape(1,-1) + + # #Calculate the curve from the mean of the parameters + # Thermal_Mean_Posterior_parameter_curve = AndieBackend._second_order_I_vs_T(data.X_predict, + # TN_mean, M0_mean, J_mean, BK_mean) + + + + #Calcualte all the predictive curves by interating through the posterior samples + print('Predicting') + # Thermal_CI_curves = np.empty((data.temperatureGrid.shape[0],0)) + Thermal_CI_curves = np.empty((np.array(data.x_predict).shape[0],0)) + thin_posterior_TN_dist = Posterior_TN[::50] + thin_posterior_M0_dist = Posterior_M0[::50] + thin_posterior_J_dist = Posterior_J[::50] + thin_posterior_BK_dist = Posterior_BK[::50] + for i, _ in enumerate(thin_posterior_TN_dist): + I_interem = AndieBackend._second_order_I_vs_T(np.array(data.x_predict), + TN=thin_posterior_TN_dist[i].reshape(1,-1), + M0=thin_posterior_M0_dist[i].reshape(1,-1), + J=thin_posterior_J_dist[i].reshape(1,-1), + BK=thin_posterior_BK_dist[i].reshape(1,-1)) + Thermal_CI_curves = np.concatenate((Thermal_CI_curves, I_interem), axis=1) + + # q = np.array([0.95, 0.05]) + # I_ci = np.quantile(Thermal_CI_curves, q, axis=1)#.reshape((-1,2)) + # Thermal_lower_curve = I_ci[1,:].reshape(-1,1) + # Thermal_upper_curve = I_ci[0,:].reshape(-1,1) + + Thermal_mean_of_posterior_curves = np.mean(Thermal_CI_curves, axis=1).reshape(-1,1) + Thermal_variance = np.var(Thermal_CI_curves, axis=1).reshape(-1,1) + # dictionary = {'Thermal Mean Posterior Parameter Curve': Thermal_Mean_Posterior_parameter_curve, + # 'Thermal Mean of Posterior Curves': Thermal_mean_of_posterior_curves, + # 'Thermal Variance': Thermal_variance, + # 'Thermal CI Curves': Thermal_CI_curves, + # 'Thermal lower Curve': Thermal_lower_curve, + # 'Thermal upper Curve': Thermal_upper_curve} + return Thermal_mean_of_posterior_curves, Thermal_variance + + @staticmethod + def sample(module, model, data): + print("sample", "-"*20) + off_flag = False + + # Thermal_next_sample = np.max(data.measuredTemperatureIdx) + 1 + # Thermal_next_sample = np.max(data.X_train[-1][1]) + 1 + Thermal_next_sample = int(data.extra_args['last_idx']) + 1 + + above_TN = data.extra_args['above_TN'] + # above_TN = data.above_TN + T_start = data.bounds[0][0] + T_stop = data.bounds[0][1] + + T_step = data.extra_args['T_step'] + T2 = np.linspace(T_start, T_stop, int((T_stop-T_start)/T_step)+1).reshape(-1,1) + + TN_best = model.state.best()[0][3] + TN_std = model.dx[3] + TN_upper = TN_best + TN_std + + data.x_predict = T2 + Thermal_mean_of_posterior_curves, Thermal_variance = module.predict(model, data) + + # def find_next_temperature(T2, T_step, above_TN, Thermal_next_sample, Thermal_posterior_results, Thermal_predictions): + + + + next_sample_acquired = False + while next_sample_acquired == False: + print("**"*20, Thermal_next_sample) + uncertainty = Thermal_variance[Thermal_next_sample] + Bravery_factor = 8.5 + if T2[Thermal_next_sample] > TN_upper: + #If the next temp is bigger than the upper confidence bound TN, take a big step + #Or if the isothermal inference did not find a peak, take a big step + above_TN += 1 + print('We are in the end game now' ) + step = np.round(5*above_TN**2.5) #Degrees to increase the temperature. + if Thermal_next_sample + int(step/T_step) < T2.shape[0]: + #If the index after the step would still be within the search space, choose that temp. + Thermal_next_sample = Thermal_next_sample + int(step/T_step) + + else: + print('You have reached your destination.') + off_flag = True + + next_sample_acquired = True + #Optional further condition to take more data points within confidence interval of the predicted TN + # elif T2[Thermal_next_sample] >= TN_lower: + # #If the next temp is bigger than the lower confidence bound of TN, take a small step + # x_step = 0.5 #Degrees to increase the temperature + # Thermal_next_sample = Thermal_next_sample + int(x_step/T_step) + # next_sample_acquired = True + # print('Theres nothing for it, you got to take the next step' ) + elif uncertainty < Bravery_factor*Thermal_mean_of_posterior_curves[Thermal_next_sample]/100: #With the instrument error scaling factor + Thermal_next_sample = Thermal_next_sample+1 + else: + print('Acquired Next Temp on Main path') + next_sample_acquired = True + + + if off_flag: + Thermal_next_sample = T2.shape[0] - 1 + + + ## these updates may not be necessary (they are on the server side) + ## the updates on the client side are from the outputs of this function + ## the client updates are handled in handle_next_points() in the clinet file + data.extra_args['above_TN'] = above_TN + data.extra_args['last_idx'] = Thermal_next_sample + + return [T2[Thermal_next_sample], above_TN, Thermal_next_sample] + + + diff --git a/src/dial_service/core.py b/src/dial_service/core.py index ecee962..60428f5 100644 --- a/src/dial_service/core.py +++ b/src/dial_service/core.py @@ -76,4 +76,9 @@ def get_surrogate_values(data: ServersideInputPrediction) -> list[list[float]]: means, stddevs = module.predict(model, data) means = data.inverse_transform(means) transformed_stddevs = data.inverse_transform(stddevs, is_stddev=True) + + means = means.flatten() + transformed_stddevs = transformed_stddevs.flatten() + stddevs = stddevs.flatten() + return [means.tolist(), transformed_stddevs.tolist(), stddevs.tolist()] diff --git a/src/dial_service/service_specific_dataclasses.py b/src/dial_service/service_specific_dataclasses.py index 4470014..33fc9f9 100644 --- a/src/dial_service/service_specific_dataclasses.py +++ b/src/dial_service/service_specific_dataclasses.py @@ -14,8 +14,11 @@ def _get_permitted_backends() -> tuple[str, ...]: """ import importlib.util + # available_backends = [ + # bkend for bkend in _POSSIBLE_BACKENDS if importlib.util.find_spec(bkend) is not None + # ] available_backends = [ - bkend for bkend in _POSSIBLE_BACKENDS if importlib.util.find_spec(bkend) is not None + bkend for bkend in _POSSIBLE_BACKENDS ] if not available_backends: # TODO - provide explicit installation instructions in this message From 59f1c954d3a4c2b75f509b78133d526ff7861d80 Mon Sep 17 00:00:00 2001 From: Paul Laiu Date: Thu, 15 May 2025 13:41:48 -0400 Subject: [PATCH 2/8] attempts to fix the final point issue in ANDiE --- src/dial_service/backends/andie_backend.py | 88 +++++++++++----------- 1 file changed, 46 insertions(+), 42 deletions(-) diff --git a/src/dial_service/backends/andie_backend.py b/src/dial_service/backends/andie_backend.py index cca93c1..82d99b7 100644 --- a/src/dial_service/backends/andie_backend.py +++ b/src/dial_service/backends/andie_backend.py @@ -259,6 +259,10 @@ def sample(module, model, data): # Thermal_next_sample = np.max(data.measuredTemperatureIdx) + 1 # Thermal_next_sample = np.max(data.X_train[-1][1]) + 1 Thermal_next_sample = int(data.extra_args['last_idx']) + 1 + if Thermal_next_sample >= T2.shape[0]: + print('You have reached your destination.') + off_flag = True + above_TN = data.extra_args['above_TN'] # above_TN = data.above_TN @@ -278,50 +282,50 @@ def sample(module, model, data): # def find_next_temperature(T2, T_step, above_TN, Thermal_next_sample, Thermal_posterior_results, Thermal_predictions): + if (off_flag == False): + next_sample_acquired = False + while next_sample_acquired == False: + print("**"*20, Thermal_next_sample) + uncertainty = Thermal_variance[Thermal_next_sample] + Bravery_factor = 8.5 + if T2[Thermal_next_sample] > TN_upper: + #If the next temp is bigger than the upper confidence bound TN, take a big step + #Or if the isothermal inference did not find a peak, take a big step + above_TN += 1 + print('We are in the end game now' ) + step = np.round(5*above_TN**2.5) #Degrees to increase the temperature. + if Thermal_next_sample + int(step/T_step) < T2.shape[0]: + #If the index after the step would still be within the search space, choose that temp. + Thermal_next_sample = Thermal_next_sample + int(step/T_step) + + else: + print('You have reached your destination.') + off_flag = True + + next_sample_acquired = True + #Optional further condition to take more data points within confidence interval of the predicted TN + # elif T2[Thermal_next_sample] >= TN_lower: + # #If the next temp is bigger than the lower confidence bound of TN, take a small step + # x_step = 0.5 #Degrees to increase the temperature + # Thermal_next_sample = Thermal_next_sample + int(x_step/T_step) + # next_sample_acquired = True + # print('Theres nothing for it, you got to take the next step' ) + elif uncertainty < Bravery_factor*Thermal_mean_of_posterior_curves[Thermal_next_sample]/100: #With the instrument error scaling factor + Thermal_next_sample = Thermal_next_sample+1 + else: + print('Acquired Next Temp on Main path') + next_sample_acquired = True - next_sample_acquired = False - while next_sample_acquired == False: - print("**"*20, Thermal_next_sample) - uncertainty = Thermal_variance[Thermal_next_sample] - Bravery_factor = 8.5 - if T2[Thermal_next_sample] > TN_upper: - #If the next temp is bigger than the upper confidence bound TN, take a big step - #Or if the isothermal inference did not find a peak, take a big step - above_TN += 1 - print('We are in the end game now' ) - step = np.round(5*above_TN**2.5) #Degrees to increase the temperature. - if Thermal_next_sample + int(step/T_step) < T2.shape[0]: - #If the index after the step would still be within the search space, choose that temp. - Thermal_next_sample = Thermal_next_sample + int(step/T_step) - else: - print('You have reached your destination.') - off_flag = True - - next_sample_acquired = True - #Optional further condition to take more data points within confidence interval of the predicted TN - # elif T2[Thermal_next_sample] >= TN_lower: - # #If the next temp is bigger than the lower confidence bound of TN, take a small step - # x_step = 0.5 #Degrees to increase the temperature - # Thermal_next_sample = Thermal_next_sample + int(x_step/T_step) - # next_sample_acquired = True - # print('Theres nothing for it, you got to take the next step' ) - elif uncertainty < Bravery_factor*Thermal_mean_of_posterior_curves[Thermal_next_sample]/100: #With the instrument error scaling factor - Thermal_next_sample = Thermal_next_sample+1 - else: - print('Acquired Next Temp on Main path') - next_sample_acquired = True - - - if off_flag: - Thermal_next_sample = T2.shape[0] - 1 - - - ## these updates may not be necessary (they are on the server side) - ## the updates on the client side are from the outputs of this function - ## the client updates are handled in handle_next_points() in the clinet file - data.extra_args['above_TN'] = above_TN - data.extra_args['last_idx'] = Thermal_next_sample + if off_flag: + Thermal_next_sample = T2.shape[0] - 1 + + + ## these updates may not be necessary (they are on the server side) + ## the updates on the client side are from the outputs of this function + ## the client updates are handled in handle_next_points() in the clinet file + data.extra_args['above_TN'] = above_TN + data.extra_args['last_idx'] = Thermal_next_sample return [T2[Thermal_next_sample], above_TN, Thermal_next_sample] From 096fe1310a610b2512e39f007108372bcbc6f592 Mon Sep 17 00:00:00 2001 From: Paul Laiu Date: Thu, 15 May 2025 16:30:23 -0400 Subject: [PATCH 3/8] updated andie backend and client --- scripts/andie_mock.py | 6 ++++++ src/dial_service/backends/andie_backend.py | 18 ++++++++++-------- src/dial_service/core.py | 4 ---- 3 files changed, 16 insertions(+), 12 deletions(-) diff --git a/scripts/andie_mock.py b/scripts/andie_mock.py index 066e701..404b2ed 100644 --- a/scripts/andie_mock.py +++ b/scripts/andie_mock.py @@ -242,6 +242,12 @@ def handle_next_points(self, payload): coord_str = ', '.join([f'{coord:.2f}' for coord in optimal_coords]) print(f'Optimal simulated datapoint at ({coord_str}), y={y_opt:.3f}\n') + + plt.figure() + x_plot = np.linspace(0,310,300) + plt.plot(x_plot,peak_val_at_T(x_plot)) + plt.plot(self.dataset_x, self.dataset_y, 'o') + plt.savefig('andie.png') def graph(self): plt.clf() diff --git a/src/dial_service/backends/andie_backend.py b/src/dial_service/backends/andie_backend.py index 82d99b7..848fd8c 100644 --- a/src/dial_service/backends/andie_backend.py +++ b/src/dial_service/backends/andie_backend.py @@ -180,7 +180,7 @@ def train_model(data): print("initial chisq", problem.chisq_str()) result = fit(problem, method=method, xtol=1e-6, ftol=1e-8, - samples = 10, + samples = 100000, burn = 100, pop = 10, init = 'eps', @@ -249,7 +249,7 @@ def predict(posterior_results, data): # 'Thermal CI Curves': Thermal_CI_curves, # 'Thermal lower Curve': Thermal_lower_curve, # 'Thermal upper Curve': Thermal_upper_curve} - return Thermal_mean_of_posterior_curves, Thermal_variance + return Thermal_mean_of_posterior_curves.flatten(), Thermal_variance.flatten() @staticmethod def sample(module, model, data): @@ -259,10 +259,7 @@ def sample(module, model, data): # Thermal_next_sample = np.max(data.measuredTemperatureIdx) + 1 # Thermal_next_sample = np.max(data.X_train[-1][1]) + 1 Thermal_next_sample = int(data.extra_args['last_idx']) + 1 - if Thermal_next_sample >= T2.shape[0]: - print('You have reached your destination.') - off_flag = True - + above_TN = data.extra_args['above_TN'] # above_TN = data.above_TN @@ -281,6 +278,10 @@ def sample(module, model, data): # def find_next_temperature(T2, T_step, above_TN, Thermal_next_sample, Thermal_posterior_results, Thermal_predictions): + if Thermal_next_sample >= T2.shape[0]: + print('You have reached your destination.') + off_flag = True + if (off_flag == False): next_sample_acquired = False @@ -317,8 +318,9 @@ def sample(module, model, data): next_sample_acquired = True - if off_flag: - Thermal_next_sample = T2.shape[0] - 1 + if off_flag: + Thermal_next_sample = T2.shape[0] - 1 + ## these updates may not be necessary (they are on the server side) diff --git a/src/dial_service/core.py b/src/dial_service/core.py index 60428f5..8f04f0e 100644 --- a/src/dial_service/core.py +++ b/src/dial_service/core.py @@ -77,8 +77,4 @@ def get_surrogate_values(data: ServersideInputPrediction) -> list[list[float]]: means = data.inverse_transform(means) transformed_stddevs = data.inverse_transform(stddevs, is_stddev=True) - means = means.flatten() - transformed_stddevs = transformed_stddevs.flatten() - stddevs = stddevs.flatten() - return [means.tolist(), transformed_stddevs.tolist(), stddevs.tolist()] From 904aebfc469a5caf46f2a53b7c037dd636f7f0a1 Mon Sep 17 00:00:00 2001 From: Paul Laiu Date: Fri, 16 May 2025 11:05:59 -0400 Subject: [PATCH 4/8] simplified ANDiE backend workflow; fixed I/O shap issue; fixed end point issue; added client stopping criterion when the last grid point is sampled --- scripts/andie_mock.py | 106 +++++++-------------- src/dial_service/backends/andie_backend.py | 14 ++- 2 files changed, 44 insertions(+), 76 deletions(-) diff --git a/scripts/andie_mock.py b/scripts/andie_mock.py index 404b2ed..8a9a029 100644 --- a/scripts/andie_mock.py +++ b/scripts/andie_mock.py @@ -69,6 +69,9 @@ def __init__(self, service_destination: str): self.bounds = np.array([[T_start, T_stop]]) self.num_dims = len(self.bounds) + + self.plot_results = True + self.x_raw = np.array([[ T_start]]) self.x_test = np.array(T_grid) self.y_raw = peak_val_at_T(self.x_raw) @@ -124,7 +127,7 @@ def callback(self, operation: str, **kwargs) -> IntersectClientCallback: print("send", operation) next_payload = None - self.at_grids = kwargs.get('at_grids', True) + # self.at_grids = kwargs.get('at_grids', True) if operation == 'dial.get_surrogate_values': @@ -180,18 +183,22 @@ def __call__(self, _source: str, operation: str, _has_error: bool, if operation == 'dial.initialize_workflow': self.workflow_id = payload print("\n","--"*20,"\n") - return self.callback('dial.get_surrogate_values') + # return self.callback('dial.get_surrogate_values') + self.niter += 1 + + return self.callback('dial.get_next_point') + # TODO: repplace this with update_workflow_with_data that can also train the model # ----------------- Active learning loop ----------------- - if operation == 'dial.get_surrogate_values': - self.handle_surrogate_values(payload) + # if operation == 'dial.get_surrogate_values': + # self.handle_surrogate_values(payload) - if self.at_grids: - print(f"Step {self.niter}") - return self.callback('dial.get_surrogate_values', at_grids=False) - else: - return self.callback('dial.get_next_point') + # if self.at_grids: + # print(f"Step {self.niter}") + # return self.callback('dial.get_surrogate_values', at_grids=False) + # else: + # return self.callback('dial.get_next_point') elif operation == 'dial.get_next_point': self.handle_next_points(payload) @@ -199,7 +206,7 @@ def __call__(self, _source: str, operation: str, _has_error: bool, elif operation == 'dial.update_workflow_with_data': self.niter += 1 - return self.callback('dial.get_surrogate_values') + return self.callback('dial.get_next_point') else: raise IntersectCallbackError(operation, payload) @@ -243,68 +250,23 @@ def handle_next_points(self, payload): print(f'Optimal simulated datapoint at ({coord_str}), y={y_opt:.3f}\n') - plt.figure() - x_plot = np.linspace(0,310,300) - plt.plot(x_plot,peak_val_at_T(x_plot)) - plt.plot(self.dataset_x, self.dataset_y, 'o') - plt.savefig('andie.png') - def graph(self): - - plt.clf() - - fig, axs = plt.subplots(2, 1, figsize=(10, 8), sharex=True) - - # First subplot: Mean and variance with training data - axs[0].plot(self.x_grid, self.mean_grid, label='Mean Prediction') - axs[0].fill_between( - self.x_grid[:, 0], - self.mean_grid + 2 * self.variance_grid, - self.mean_grid - 2 * self.variance_grid, - alpha=0.5, - label='Confidence Interval' - ) - axs[0].scatter( - np.array(self.dataset_x)[:-1, 0], - np.array(self.dataset_y)[:-1], - color='black', - marker='o', - label='Training Data' - ) - if self.x_next is not None: - axs[0].axvline( - x=self.x_next[0], - color='red', - linestyle='--') - axs[0].set_ylabel('Response, y') - axs[0].legend() - axs[0].grid(True) - - # Second subplot: Acquisition function - if self.strategy_args is not None: - if self.mean_grid is not None and self.variance_grid is not None: - exploit = self.strategy_args.get('exploit', 0.0) - explore = self.strategy_args.get('explore', 1.0) - acquisition_values = exploit * self.mean_grid + explore * np.sqrt(self.variance_grid) - else: - acquisition_values = np.zeros_like(self.x_grid) - - axs[1].plot(self.x_grid, acquisition_values) - if self.x_next is not None: - axs[1].axvline( - x=self.x_next[0], - color='red', - linestyle='--', - label='Next Point (x_next)' - ) - axs[1].set_xlabel('Features, x') - axs[1].set_ylabel('Acquisition Value') - axs[1].legend() - axs[1].grid(True) - - plt.tight_layout() - plt.savefig('graph.png') - plt.close(fig) - + if self.plot_results: + plt.figure() + x_plot = np.linspace(0,310,300) + plt.plot(x_plot,peak_val_at_T(x_plot)) + plt.plot(self.dataset_x, self.dataset_y, 'o') + plt.savefig('andie.png') + + # end client when maximum iteration reached + if self.niter >= self.max_iter: + print("Maximum iteration reached") + raise IntersectCallbackEnd() + # end client when last point on the grid reached + if self.last_idx == len(self.x_test) - 1: + print("Last grid point reached") + raise IntersectCallbackEnd() + + if __name__ == '__main__': parser = argparse.ArgumentParser(description='Automated client') diff --git a/src/dial_service/backends/andie_backend.py b/src/dial_service/backends/andie_backend.py index 848fd8c..a179074 100644 --- a/src/dial_service/backends/andie_backend.py +++ b/src/dial_service/backends/andie_backend.py @@ -171,11 +171,14 @@ def train_model(data): # X_grid = [temperature for temperature in list(np.array(data.X_train)[0,:])] X_grid = np.array(data.X_train) - Y_grid = np.array(data.Y_train) + Y_grid = np.array(data.Y_train).reshape(-1,1) # X_grid = data.measuredTemperatures # Y_grid = data.measuredTemperatureDependentIntensities problem = AndieBackend._Second_order_model(X_grid, Y_grid) + print('Sampled temperature,', X_grid.flatten()) + print('Sampled intensity,', Y_grid.flatten()) + method = 'dream' print("initial chisq", problem.chisq_str()) @@ -221,8 +224,6 @@ def predict(posterior_results, data): #Calcualte all the predictive curves by interating through the posterior samples - print('Predicting') - # Thermal_CI_curves = np.empty((data.temperatureGrid.shape[0],0)) Thermal_CI_curves = np.empty((np.array(data.x_predict).shape[0],0)) thin_posterior_TN_dist = Posterior_TN[::50] thin_posterior_M0_dist = Posterior_M0[::50] @@ -286,7 +287,6 @@ def sample(module, model, data): if (off_flag == False): next_sample_acquired = False while next_sample_acquired == False: - print("**"*20, Thermal_next_sample) uncertainty = Thermal_variance[Thermal_next_sample] Bravery_factor = 8.5 if T2[Thermal_next_sample] > TN_upper: @@ -313,14 +313,20 @@ def sample(module, model, data): # print('Theres nothing for it, you got to take the next step' ) elif uncertainty < Bravery_factor*Thermal_mean_of_posterior_curves[Thermal_next_sample]/100: #With the instrument error scaling factor Thermal_next_sample = Thermal_next_sample+1 + if Thermal_next_sample >= T2.shape[0]: + print('You have reached your destination.') + off_flag = True + next_sample_acquired = True else: print('Acquired Next Temp on Main path') next_sample_acquired = True + if off_flag: Thermal_next_sample = T2.shape[0] - 1 + print("*"*20, 'Next sample:', T2[Thermal_next_sample]) ## these updates may not be necessary (they are on the server side) From b0ff3ab38925102bb84402f6127f375710f63d46 Mon Sep 17 00:00:00 2001 From: Paul Laiu Date: Fri, 23 May 2025 10:31:10 -0400 Subject: [PATCH 5/8] updated prior values and the bravery factor in andie_backend.py for temperature range 100-500K; updated temperature grid and mock intensity-temperature profile in andie_mock.py for temperature range 100-500K --- scripts/andie_mock.py | 12 +++++-- src/dial_service/backends/andie_backend.py | 39 +++++++++++++++++----- 2 files changed, 40 insertions(+), 11 deletions(-) diff --git a/scripts/andie_mock.py b/scripts/andie_mock.py index 8a9a029..d09649a 100644 --- a/scripts/andie_mock.py +++ b/scripts/andie_mock.py @@ -32,7 +32,9 @@ logger = logging.getLogger(__name__) def peak_val_at_T(T: np.ndarray) -> np.ndarray: - result = 40.0 * 1.0/(1.0+np.exp((T-50)/5)) + 35.0 + # result = 40.0 * 1.0/(1.0+np.exp((T-50)/5)) + 35.0 + result = 70.0 * 1.0/(1.0+np.exp((T-200)/7)) + 30.0 + logger.debug(result) return result @@ -54,9 +56,13 @@ def __init__(self, service_destination: str): #Maximum Measurments in Temperature: Temperature_Loops = 30 - T_start = 5.0 #Kelvin - T_stop = 300.0 #Kelvin + # T_start = 5.0 #Kelvin + # T_stop = 300.0 #Kelvin + # T_step = .5 #Kelvin + T_start = 100.0 #Kelvin + T_stop = 500.0 #Kelvin T_step = .5 #Kelvin + T_grid = np.linspace(T_start, T_stop, int((T_stop-T_start)/T_step)+1).reshape(-1,1) self.T_step = T_step diff --git a/src/dial_service/backends/andie_backend.py b/src/dial_service/backends/andie_backend.py index a179074..001b20a 100644 --- a/src/dial_service/backends/andie_backend.py +++ b/src/dial_service/backends/andie_backend.py @@ -21,14 +21,37 @@ class AndieBackend( ### --- These are hyperparameters --- ### ### --- Define the priors on the Isothermal parameters --- ### + ### --- OLD hyperparameters for temperature range 0.5-100 K --- ### + # #Transition Temperature + # TN_guess = 60.0 + # TN_std = 20 + # TN_limits = (0.01, 200.0) + + # #Background + # BK_guess = 70.0 + # BK_std = 20.0 + # BK_limits = (0.0, 150.0) + + # #Second Order Scale + # M0_guess = 18.0 + # M0_std = 3.0 + # M0_limits = (0.01, 35.0) + + # #Total Anglular Momentum + # J_guess = 0.6 + # J_std = 0.5 + # J_limits = (0.01, 14.0) + + ### --- NEW hyperparameters for temperature range 100-500 K --- ### + #Transition Temperature - TN_guess = 60.0 - TN_std = 20 - TN_limits = (0.01, 200.0) + TN_guess = 250.0 + TN_std = 80 + TN_limits = (90.0, 500.0) #Background - BK_guess = 70.0 - BK_std = 20.0 + BK_guess = 75.0 + BK_std = 50.0 BK_limits = (0.0, 150.0) #Second Order Scale @@ -176,8 +199,8 @@ def train_model(data): # Y_grid = data.measuredTemperatureDependentIntensities problem = AndieBackend._Second_order_model(X_grid, Y_grid) - print('Sampled temperature,', X_grid.flatten()) - print('Sampled intensity,', Y_grid.flatten()) + print('Sampled temperature:', X_grid.flatten()) + print('Sampled intensity:', Y_grid.flatten()) method = 'dream' @@ -288,7 +311,7 @@ def sample(module, model, data): next_sample_acquired = False while next_sample_acquired == False: uncertainty = Thermal_variance[Thermal_next_sample] - Bravery_factor = 8.5 + Bravery_factor = 4.0 if T2[Thermal_next_sample] > TN_upper: #If the next temp is bigger than the upper confidence bound TN, take a big step #Or if the isothermal inference did not find a peak, take a big step From 454744eb156a6b370fe65ea6cd24f9db1dceb823 Mon Sep 17 00:00:00 2001 From: Lance-Drane Date: Fri, 23 May 2025 13:33:35 -0400 Subject: [PATCH 6/8] move bumps and statsmodels dependencies from dev deps to project deps Signed-off-by: Lance-Drane --- pdm.lock | 24 ++++++++++++------------ pyproject.toml | 8 ++------ 2 files changed, 14 insertions(+), 18 deletions(-) diff --git a/pdm.lock b/pdm.lock index 7c25fa8..a822eaf 100644 --- a/pdm.lock +++ b/pdm.lock @@ -5,7 +5,7 @@ groups = ["default", "dev", "lint", "test"] strategy = ["inherit_metadata"] lock_version = "4.5.0" -content_hash = "sha256:53c785e022328218799c1a8af52a53be7159165d21459197218fc73dfa39244d" +content_hash = "sha256:17b29d9385eda10783ba9efdcab725cae8fcd7fbfdda110a83a87e8d4600532f" [[metadata.targets]] requires_python = ">=3.10" @@ -106,7 +106,7 @@ files = [ name = "bumps" version = "0.9.3" summary = "Data fitting with bayesian uncertainty analysis" -groups = ["dev"] +groups = ["default"] files = [ {file = "bumps-0.9.3-py3-none-any.whl", hash = "sha256:269f5e86c50fde2b80015b8c9aff9b06eb045be675bdef065d2897670c5b56e5"}, {file = "bumps-0.9.3.tar.gz", hash = "sha256:3295298f7fe1b2392bb2ff88c7a0ae69d77a7698af580a81570052ca05ba530f"}, @@ -1025,7 +1025,7 @@ name = "numpy" version = "1.26.4" requires_python = ">=3.9" summary = "Fundamental package for array computing in Python" -groups = ["default", "dev"] +groups = ["default"] files = [ {file = "numpy-1.26.4-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:9ff0f4f29c51e2803569d7a51c2304de5554655a60c5d776e35b4a41413830d0"}, {file = "numpy-1.26.4-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:2e4ee3380d6de9c9ec04745830fd9e2eccb3e6cf790d39d7b98ffd19b0dd754a"}, @@ -1088,7 +1088,7 @@ name = "packaging" version = "24.1" requires_python = ">=3.8" summary = "Core utilities for Python packages" -groups = ["default", "dev", "test"] +groups = ["default", "test"] files = [ {file = "packaging-24.1-py3-none-any.whl", hash = "sha256:5b8f2217dbdbd2f7f384c41c628544e6d52f2d0f53c6d0c3ea61aa5d1d7ff124"}, {file = "packaging-24.1.tar.gz", hash = "sha256:026ed72c8ed3fcce5bf8950572258698927fd1dbda10a5e981cdf0ac37f4f002"}, @@ -1108,7 +1108,7 @@ name = "pandas" version = "2.2.3" requires_python = ">=3.9" summary = "Powerful data structures for data analysis, time series, and statistics" -groups = ["dev"] +groups = ["default"] dependencies = [ "numpy>=1.22.4; python_version < \"3.11\"", "numpy>=1.23.2; python_version == \"3.11\"", @@ -1160,7 +1160,7 @@ name = "patsy" version = "1.0.1" requires_python = ">=3.6" summary = "A Python package for describing statistical models and for building design matrices." -groups = ["dev"] +groups = ["default"] dependencies = [ "numpy>=1.4", ] @@ -1518,7 +1518,7 @@ name = "python-dateutil" version = "2.9.0.post0" requires_python = "!=3.0.*,!=3.1.*,!=3.2.*,>=2.7" summary = "Extensions to the standard Python datetime module" -groups = ["default", "dev"] +groups = ["default"] dependencies = [ "six>=1.5", ] @@ -1531,7 +1531,7 @@ files = [ name = "pytz" version = "2025.2" summary = "World timezone definitions, modern and historical" -groups = ["dev"] +groups = ["default"] files = [ {file = "pytz-2025.2-py2.py3-none-any.whl", hash = "sha256:5ddf76296dd8c44c26eb8f4b6f35488f3ccbf6fbbd7adee0b7262d43f0ec2f00"}, {file = "pytz-2025.2.tar.gz", hash = "sha256:360b9e3dbb49a209c21ad61809c7fb453643e048b38924c765813546746e81c3"}, @@ -1778,7 +1778,7 @@ name = "scipy" version = "1.14.1" requires_python = ">=3.10" summary = "Fundamental algorithms for scientific computing in Python" -groups = ["default", "dev"] +groups = ["default"] dependencies = [ "numpy<2.3,>=1.23.5", ] @@ -1823,7 +1823,7 @@ name = "six" version = "1.16.0" requires_python = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*" summary = "Python 2 and 3 compatibility utilities" -groups = ["default", "dev"] +groups = ["default"] files = [ {file = "six-1.16.0-py2.py3-none-any.whl", hash = "sha256:8abb2f1d86890a2dfb989f9a77cfcfd3e47c2a354b01111771326f8aa26e0254"}, {file = "six-1.16.0.tar.gz", hash = "sha256:1e61c37477a1626458e36f7b1d82aa5c9b094fa4802892072e49de9c60c4c926"}, @@ -1834,7 +1834,7 @@ name = "statsmodels" version = "0.14.4" requires_python = ">=3.9" summary = "Statistical computations and models for Python" -groups = ["dev"] +groups = ["default"] dependencies = [ "numpy<3,>=1.22.3", "packaging>=21.3", @@ -1945,7 +1945,7 @@ name = "tzdata" version = "2025.2" requires_python = ">=2" summary = "Provider of IANA time zone data" -groups = ["dev"] +groups = ["default"] files = [ {file = "tzdata-2025.2-py2.py3-none-any.whl", hash = "sha256:1a403fada01ff9221ca8044d701868fa132215d84beb92242d9acd2147f667a8"}, {file = "tzdata-2025.2.tar.gz", hash = "sha256:b60a638fcc0daffadf82fe0f57e53d06bdec2f36c4df66280ae79bce6bd6f2b9"}, diff --git a/pyproject.toml b/pyproject.toml index e1fb61b..f9c6d82 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,6 +17,8 @@ dependencies = [ "scipy>=1.12.0,<2.0.0", "gpax>=0.1.8", # TODO consider making an optional dependency group "pymongo>=4.12.1", # TODO - this is only needed for dial_service, dial_dataclass can use a simple fixture for ObjectID representation which allows it to skip this dependency + "bumps>=0.9.3", + "statsmodels>=0.14.4", ] [tool.pdm.dev-dependencies] @@ -130,9 +132,3 @@ plugins = ["pydantic.mypy"] [build-system] requires = ["setuptools", "setuptools-scm"] build-backend = "setuptools.build_meta" - -[dependency-groups] -dev = [ - "statsmodels>=0.14.4", - "bumps>=0.9.3", -] From cb3642f1d592ca1876514d3c75a033a1c1d936cb Mon Sep 17 00:00:00 2001 From: Paul Laiu Date: Tue, 2 Dec 2025 12:36:12 -0500 Subject: [PATCH 7/8] initial attempt to generalize the model fitting framework in andie --- .../backends/andie_backend_general.py | 324 ++++++++++++++++++ .../utilities/parametric_models.py | 108 ++++++ 2 files changed, 432 insertions(+) create mode 100644 src/dial_service/backends/andie_backend_general.py create mode 100644 src/dial_service/utilities/parametric_models.py diff --git a/src/dial_service/backends/andie_backend_general.py b/src/dial_service/backends/andie_backend_general.py new file mode 100644 index 0000000..86fed82 --- /dev/null +++ b/src/dial_service/backends/andie_backend_general.py @@ -0,0 +1,324 @@ +"""NOTE: This file should not be imported in application code except dynamically via the get_backend_module function in __init__.py .""" + +import numpy as np +import scipy as sp +from scipy.optimize import OptimizeResult + +from bumps.names import Curve, FitProblem, fit +from bumps.formatnum import format_uncertainty +from bumps.bounds import BoundedNormal + +from ..utilities import parametric_models +from . import AbstractBackend + + +_MODELS_ANDIE = { + 'ANDIE_second_order_I_vs_T': parametric_models.ANDIE_second_order_I_vs_T +} + +class AndieBackend( + AbstractBackend[OptimizeResult, None, tuple[np.ndarray, np.ndarray]] +): + + @staticmethod + def get_parameter_names(data): + + return [parname for parname in data.extra_args['priors'].keys()] + + @staticmethod + def get_prior_distributions(data): + ### --- Define the priors on the Isothermal parameters --- ### + + ### --- OLD hyperparameters for temperature range 0.5-100 K --- ### + # #Transition Temperature + # TN_guess = 60.0 + # TN_std = 20 + # TN_limits = (0.01, 200.0) + + # #Background + # BK_guess = 70.0 + # BK_std = 20.0 + # BK_limits = (0.0, 150.0) + + # #Second Order Scale + # M0_guess = 18.0 + # M0_std = 3.0 + # M0_limits = (0.01, 35.0) + + # #Total Anglular Momentum + # J_guess = 0.6 + # J_std = 0.5 + # J_limits = (0.01, 14.0) + + ### --- NEW hyperparameters for temperature range 100-500 K --- ### + + + # #Transition Temperature + # TN_guess = 250.0 + # TN_std = 80 + # TN_limits = (90.0, 500.0) + + # #Background + # BK_guess = 75.0 + # BK_std = 50.0 + # BK_limits = (0.0, 150.0) + + # #Second Order Scale + # M0_guess = 18.0 + # M0_std = 3.0 + # M0_limits = (0.01, 35.0) + + # #Total Anglular Momentum + # J_guess = 0.6 + # J_std = 0.5 + # J_limits = (0.01, 14.0) + + + # #Define the prior distributions + # TN_dist = BoundedNormal(TN_guess, TN_std, limits=TN_limits) + # M0_dist = BoundedNormal(M0_guess, M0_std, limits=M0_limits) + # J_dist = BoundedNormal(J_guess, J_std, limits=J_limits) + # BK_dist = BoundedNormal(BK_guess, BK_std, limits=BK_limits) + + # ### --- End of hyperparameters --- ### + + # guesses = [TN_guess, M0_guess, J_guess, BK_guess] + # priors = [TN_dist, M0_dist, J_dist, BK_dist] + + prior_dict = data.extra_args['priors'] + + guesses = {} + priors = {} + for par_name, prior in prior_dict.items(): + prior_dist = BoundedNormal(prior['guess'], prior['std'], limits=prior['limits']) + guesses[par_name] = prior['guess'] + priors[par_name] = prior_dist + return guesses, priors + + @staticmethod + def get_parametric_model(data): + model_name = data.model_name.lower() + + if model_name not in _MODELS_ANDIE: + raise ValueError(f"Unknown startegy {model_name}") + + return _MODELS_ANDIE[model_name] + + @staticmethod + def _Second_order_model(X_grid, Y_grid, model_to_be_optimized, guesses, priors): + + M = Curve(model_to_be_optimized, X_grid, Y_grid, + dy = np.sqrt(Y_grid/100), + pars = guesses) + + for par_name, prior_dist in priors.items(): + M.__getattribute__(par_name).dist = prior_dist + + if len(X_grid) < 5 : + problem = FitProblem(M, partial = True) + else: + problem = FitProblem(M, partial = False) + + return problem + + @staticmethod + def train_model(data): + print("train", "-"*20) + + # X_grid = [temperature for temperature in list(np.array(data.X_train)[0,:])] + X_grid = np.array(data.X_train) + Y_grid = np.array(data.Y_train).reshape(-1,1) + # X_grid = data.measuredTemperatures + # Y_grid = data.measuredTemperatureDependentIntensities + + Andie_model = AndieBackend.get_parametric_model(data) + Andie_guesses, Andie_priors = AndieBackend.get_prior_distributions(data) + problem = AndieBackend._Second_order_model(X_grid, Y_grid, Andie_model, Andie_guesses, Andie_priors) + + print('Sampled temperature:', X_grid.flatten()) + print('Sampled intensity:', Y_grid.flatten()) + + method = 'dream' + + print("initial chisq", problem.chisq_str()) + result = fit(problem, method=method, xtol=1e-6, ftol=1e-8, + samples = 100000, + burn = 100, + pop = 10, + init = 'eps', + thin = 1, + alpha = 0.1, + outliers = 'none', + trim = False, + steps = 0) + print("final chisq", problem.chisq_str()) + for k, v, dv in zip(problem.labels(), result.x, result.dx): + print(k, ":", format_uncertainty(v, dv)) + + return result + + @staticmethod + def predict(posterior_results, data): + print("predict", "-"*20) + + #Extract the Posterior of the parameters + draw = posterior_results.state.draw() + + parameter_names = AndieBackend.get_parameter_names(data) + + parameter_names.sort() + + #Posterior results are in alphabetical order + Posterior = {} + for i, par_name in enumerate(parameter_names): + Posterior[par_name] = draw.points[:,i] + + # + # Posterior_TN = draw.points[:,3] + # Posterior_M0 = draw.points[:,2] + # Posterior_J = draw.points[:,1] + # Posterior_BK = draw.points[:,0] + + # #Calcualte the mean of the parameters + # TN_mean = np.mean(Posterior_TN).reshape(1,-1) + # M0_mean = np.mean(Posterior_M0).reshape(1,-1) + # J_mean = np.mean(Posterior_J).reshape(1,-1) + # BK_mean = np.mean(Posterior_BK).reshape(1,-1) + + # #Calculate the curve from the mean of the parameters + # Thermal_Mean_Posterior_parameter_curve = AndieBackend._second_order_I_vs_T(data.X_predict, + # TN_mean, M0_mean, J_mean, BK_mean) + + + + #Calcualte all the predictive curves by interating through the posterior samples + Thermal_CI_curves = np.empty((np.array(data.x_predict).shape[0],0)) + + thin_posterior = {} + for par_name in parameter_names: + thin_posterior[par_name] = Posterior[par_name][::50] + + # thin_posterior_TN_dist = Posterior_TN[::50] + # thin_posterior_M0_dist = Posterior_M0[::50] + # thin_posterior_J_dist = Posterior_J[::50] + # thin_posterior_BK_dist = Posterior_BK[::50] + + Andie_model = AndieBackend.get_parametric_model(data) + + for i, _ in enumerate(thin_posterior[parameter_names[0]]): + + drawn_params = {} + for par_name in parameter_names: + drawn_params[par_name] = thin_posterior[par_name][i].reshape(1,-1) + + I_interem = Andie_model(np.array(data.x_predict), model_params=drawn_params) + + # I_interem = AndieBackend._second_order_I_vs_T(np.array(data.x_predict), + # TN=thin_posterior_TN_dist[i].reshape(1,-1), + # M0=thin_posterior_M0_dist[i].reshape(1,-1), + # J=thin_posterior_J_dist[i].reshape(1,-1), + # BK=thin_posterior_BK_dist[i].reshape(1,-1)) + Thermal_CI_curves = np.concatenate((Thermal_CI_curves, I_interem), axis=1) + + # q = np.array([0.95, 0.05]) + # I_ci = np.quantile(Thermal_CI_curves, q, axis=1)#.reshape((-1,2)) + # Thermal_lower_curve = I_ci[1,:].reshape(-1,1) + # Thermal_upper_curve = I_ci[0,:].reshape(-1,1) + + Thermal_mean_of_posterior_curves = np.mean(Thermal_CI_curves, axis=1).reshape(-1,1) + Thermal_variance = np.var(Thermal_CI_curves, axis=1).reshape(-1,1) + # dictionary = {'Thermal Mean Posterior Parameter Curve': Thermal_Mean_Posterior_parameter_curve, + # 'Thermal Mean of Posterior Curves': Thermal_mean_of_posterior_curves, + # 'Thermal Variance': Thermal_variance, + # 'Thermal CI Curves': Thermal_CI_curves, + # 'Thermal lower Curve': Thermal_lower_curve, + # 'Thermal upper Curve': Thermal_upper_curve} + return Thermal_mean_of_posterior_curves.flatten(), Thermal_variance.flatten() + + @staticmethod + def sample(module, model, data): + print("sample", "-"*20) + off_flag = False + + # Thermal_next_sample = np.max(data.measuredTemperatureIdx) + 1 + # Thermal_next_sample = np.max(data.X_train[-1][1]) + 1 + Thermal_next_sample = int(data.extra_args['last_idx']) + 1 + + + above_TN = data.extra_args['above_TN'] + # above_TN = data.above_TN + T_start = data.bounds[0][0] + T_stop = data.bounds[0][1] + + T_step = data.extra_args['T_step'] + T2 = np.linspace(T_start, T_stop, int((T_stop-T_start)/T_step)+1).reshape(-1,1) + + TN_best = model.state.best()[0][3] + TN_std = model.dx[3] + TN_upper = TN_best + TN_std + + data.x_predict = T2 + Thermal_mean_of_posterior_curves, Thermal_variance = module.predict(model, data) + + # def find_next_temperature(T2, T_step, above_TN, Thermal_next_sample, Thermal_posterior_results, Thermal_predictions): + + if Thermal_next_sample >= T2.shape[0]: + print('You have reached your destination.') + off_flag = True + + + if (off_flag == False): + next_sample_acquired = False + while next_sample_acquired == False: + uncertainty = Thermal_variance[Thermal_next_sample] + Bravery_factor = 4.0 + if T2[Thermal_next_sample] > TN_upper: + #If the next temp is bigger than the upper confidence bound TN, take a big step + #Or if the isothermal inference did not find a peak, take a big step + above_TN += 1 + print('We are in the end game now' ) + step = np.round(5*above_TN**2.5) #Degrees to increase the temperature. + if Thermal_next_sample + int(step/T_step) < T2.shape[0]: + #If the index after the step would still be within the search space, choose that temp. + Thermal_next_sample = Thermal_next_sample + int(step/T_step) + + else: + print('You have reached your destination.') + off_flag = True + + next_sample_acquired = True + #Optional further condition to take more data points within confidence interval of the predicted TN + # elif T2[Thermal_next_sample] >= TN_lower: + # #If the next temp is bigger than the lower confidence bound of TN, take a small step + # x_step = 0.5 #Degrees to increase the temperature + # Thermal_next_sample = Thermal_next_sample + int(x_step/T_step) + # next_sample_acquired = True + # print('Theres nothing for it, you got to take the next step' ) + elif uncertainty < Bravery_factor*Thermal_mean_of_posterior_curves[Thermal_next_sample]/100: #With the instrument error scaling factor + Thermal_next_sample = Thermal_next_sample+1 + if Thermal_next_sample >= T2.shape[0]: + print('You have reached your destination.') + off_flag = True + next_sample_acquired = True + else: + print('Acquired Next Temp on Main path') + next_sample_acquired = True + + + + if off_flag: + Thermal_next_sample = T2.shape[0] - 1 + + print("*"*20, 'Next sample:', T2[Thermal_next_sample]) + + + ## these updates may not be necessary (they are on the server side) + ## the updates on the client side are from the outputs of this function + ## the client updates are handled in handle_next_points() in the clinet file + data.extra_args['above_TN'] = above_TN + data.extra_args['last_idx'] = Thermal_next_sample + + return [T2[Thermal_next_sample], above_TN, Thermal_next_sample] + + + diff --git a/src/dial_service/utilities/parametric_models.py b/src/dial_service/utilities/parametric_models.py new file mode 100644 index 0000000..74792f6 --- /dev/null +++ b/src/dial_service/utilities/parametric_models.py @@ -0,0 +1,108 @@ +import numpy as np +from scipy.optimize import root + + + +PARAMETRIC_MODELS = { + 'ANDIE_second_order_I_vs_T': 'ANDIE_second_order_I_vs_T' +} + + +def ANDIE_second_order_I_vs_T(inputs, model_params): + + ### Example usage: + ### Temperature = np.linspace(100,500,11) + ### ANDIE_second_order_I_vs_T(Temperature, {'TN':250, 'M0':75, 'J':18, "BK":0.6}) + + ### --- The following functions can be put into a model.py file --- ### + ### --- Brillouin_J(), Wiess(), Mag_solve(), and second_order_I_vs_T() --- ### + + # Define the Brillouin Function + def Brillouin_J(x, J): + # x is a matrix + # J is a row-vector + + f1 = (2*J+1)/(2*J) #factor 1 + f2 = 1/(2*J) #factor 2 + B_J = f1/(np.tanh(f1*x)) - f2/(np.tanh(f2*x)) + return B_J + + # Define the Wiess equation (re-arranged so that it always = 0) + def Wiess(t, m, J): + # t is a collum vector + # J is a row vector + # m is a matrix + + f1 = 3.0*J/(J +1.0) # factor 1 + f2 = 1.0/t #factor 2 + + m = np.ones_like(t)*m + x = f1*m*f2 + + + W = m - Brillouin_J(x, J) + return W.reshape(-1) + + #Define a Magnetization Solver for givin t and J: + def Mag_solve(t, J): + # t is a collum vector + # J is a row vector + + def fun(x): + x= x.reshape((t.shape[0], J.shape[1])) + return Wiess(t, x, J) + + X0 = np.ones((t.shape[0], J.shape[1]), dtype=np.float64) # matrix of every t,J combination + + solution = root(fun, X0, method ='hybr') + + m = np.array(solution.x).reshape((t.shape[0], J.shape[1])) + return m + + #Define the Intensity as function of Temperature for Second order phase transitions + def second_order_I_vs_T(Temperatures, TN=None, M0=None, J=None, BK=None): + # TN, M0, J, BK are row vectors + TN = np.array(TN).reshape(1,-1) + M0 = np.array(M0).reshape(1,-1) + J = np.array(J).reshape(1,-1) + BK = np.array(BK).reshape(1,-1) + # Temperatures is a collumn vector + Temperatures = Temperatures.reshape(-1,1) + + Inverse_TN = 1/TN + + t = Temperatures @ Inverse_TN #Create a matrix of every (temperature, TN) combination + + #Find where all the temperatures below the TN. + t_low_truth = t <1.0 + t_low_truth = t_low_truth.reshape(t.shape) + t_low = t*t_low_truth + + #Because "Mag_solve" is expensive, only pass the non-zero rows of the t_low matrix. + t_low_cut = t_low[t_low.sum(axis=1) !=0] + + #Calculate the reduced magnetization for below Tn\ + m_low = Mag_solve(t_low_cut, J) #returns a matrix where each collumn of m is associated with a collum of J + + #The mag is zero above the TN, construct that matrix + m_high = np.zeros((t.shape[0]-t_low_cut.shape[0], t.shape[1])) + #Combine the m matrixies for below and above the TN + m = np.concatenate((m_low,m_high), axis=0) + + #"t_low_cut" could have a raggid bottom edge, so only keep the m values below the TN + m = m*t_low_truth + + #Calcualte the Magnetization (not reduced) + M = m * M0 #row by row mulitplication of m by M0 + + ##Calculate the intensities + #Square of Magnetization + I = M**2 + + #Add the background + I = I + BK + + return I + + return second_order_I_vs_T(inputs, **model_params) + From 651bd996dc1880003676053d47086f194a0a3248 Mon Sep 17 00:00:00 2001 From: Paul Laiu Date: Mon, 18 May 2026 16:54:31 -0400 Subject: [PATCH 8/8] added andie-type approach for general parametric models, added client script andie_general_mock.py, the caveat is that the sampling method is still hardcoded since it depends on the specific problem --- scripts/andie_general_mock.py | 281 ++++++++++++++++++ src/dial_dataclass/dial_dataclass.py | 2 +- src/dial_service/backends/__init__.py | 1 + .../backends/andie_backend_general.py | 34 ++- 4 files changed, 313 insertions(+), 5 deletions(-) create mode 100644 scripts/andie_general_mock.py diff --git a/scripts/andie_general_mock.py b/scripts/andie_general_mock.py new file mode 100644 index 0000000..0d06411 --- /dev/null +++ b/scripts/andie_general_mock.py @@ -0,0 +1,281 @@ +import argparse +import json +import logging +import os +import sys +from pathlib import Path + +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np + +from intersect_sdk import ( + INTERSECT_JSON_VALUE, + HierarchyConfig, + IntersectClient, + IntersectClientCallback, + IntersectClientConfig, + IntersectDirectMessageParams, + default_intersect_lifecycle_loop, +) + +from dial_dataclass import ( + DialInputPredictions, + DialInputSingleOtherStrategy, + DialWorkflowCreationParamsClient, + DialWorkflowDatasetUpdate, +) + +mpl.use('agg') + +logging.basicConfig(level=logging.INFO) +logger = logging.getLogger(__name__) + +# Priors passed to andie_backend_general via extra_args. +# Structure: {param_name: {guess, std, limits}} — matches get_prior_distributions(). +_PRIORS = { + 'TN': {'guess': 250.0, 'std': 80.0, 'limits': [90.0, 500.0]}, + 'M0': {'guess': 18.0, 'std': 3.0, 'limits': [0.01, 35.0]}, + 'J': {'guess': 0.6, 'std': 0.5, 'limits': [0.01, 14.0]}, + 'BK': {'guess': 75.0, 'std': 50.0, 'limits': [0.0, 150.0]}, +} + +_MODEL_NAME = 'ANDIE_second_order_I_vs_T' + + +def peak_val_at_T(T: np.ndarray) -> np.ndarray: + result = 70.0 * 1.0/(1.0+np.exp((T-200)/7)) + 30.0 + logger.debug(result) + return result + + +class IntersectCallbackError(Exception): + def __init__(self, operation, payload): + message = f"Intersect callback error during operation '{operation}'. Payload: {payload}" + super().__init__(message) + +class IntersectCallbackEnd(Exception): + def __init__(self): + message = "Stopping Intersect Calls" + super().__init__(message) + +class ActiveLearningOrchestrator: + + def __init__(self, service_destination: str): + Temperature_Loops = 30 + + T_start = 100.0 # Kelvin + T_stop = 500.0 # Kelvin + T_step = 0.5 # Kelvin + + T_grid = np.linspace(T_start, T_stop, int((T_stop-T_start)/T_step)+1).reshape(-1, 1) + + self.T_step = T_step + self.above_TN = 0 + self.last_idx = 0 + + self.bounds = np.array([[T_start, T_stop]]) + self.num_dims = len(self.bounds) + + self.plot_results = True + + self.x_raw = np.array([[T_start]]) + self.x_test = np.array(T_grid) + self.y_raw = peak_val_at_T(self.x_raw) + + self.meshgrid_size = len(T_grid) + + self.dataset_x = self.x_raw.reshape(-1, 1).tolist() + self.dataset_y = self.y_raw.reshape(-1).tolist() + self.test_points = self.x_test.reshape(-1, 1).tolist() + + self.kernel = None + self.kernel_args = None + self.backend = 'andie_general' + self.backend_args = None + self.strategy = None + self.strategy_args = None + self.niter = 0 + self.max_iter = Temperature_Loops + self.x_next = None + + # Static extra_args stored in the workflow so every backend call has access + # to model_name and priors without re-sending them each request. + self._static_extra_args = { + 'model_name': _MODEL_NAME, + 'priors': _PRIORS, + } + + self.workflow_id = None + self.service_destination = service_destination + self.initialize_workflow_message = IntersectClientCallback(messages_to_send=[ + IntersectDirectMessageParams( + destination=self.service_destination, + operation='dial.initialize_workflow', + payload=DialWorkflowCreationParamsClient( + dataset_x=self.dataset_x, + dataset_y=self.dataset_y, + bounds=self.bounds, + kernel=self.kernel, + backend=self.backend, + preprocess_standardize=False, + y_is_good=True, + seed=20, + extra_args=self._static_extra_args, + )) + ]) + + def _dynamic_extra_args(self, include_T_step: bool = False) -> dict: + """Build the per-request extra_args carrying mutable state. + + These are merged server-side with the stored workflow extra_args (which + holds model_name and priors), so only the fields that change each + iteration need to be sent here. + """ + args = { + 'above_TN': self.above_TN, + 'last_idx': self.last_idx, + } + if include_T_step: + args['T_step'] = self.T_step + return args + + def callback(self, operation: str) -> IntersectClientCallback: + print("send", operation) + + next_payload = None + + if operation == 'dial.get_surrogate_values': + _points_to_predict = np.array(self.test_points).reshape(-1, self.num_dims) + next_payload = DialInputPredictions( + workflow_id=self.workflow_id, + points_to_predict=_points_to_predict, + extra_args=self._dynamic_extra_args(), + ) + + elif operation == 'dial.get_next_point': + next_payload = DialInputSingleOtherStrategy( + workflow_id=self.workflow_id, + strategy=self.strategy, + strategy_args=self.strategy_args, + bounds=self.bounds.tolist(), + extra_args=self._dynamic_extra_args(include_T_step=True), + ) + + elif operation == 'dial.update_workflow_with_data': + next_payload = DialWorkflowDatasetUpdate( + workflow_id=self.workflow_id, + next_x=self.dataset_x[-1], + next_y=self.dataset_y[-1], + ) + + else: + err_msg = f'Unknown operation received: {operation}' + raise Exception(err_msg) # noqa: TRY002 + + return IntersectClientCallback(messages_to_send=[ + IntersectDirectMessageParams( + destination=self.service_destination, + operation=operation, + payload=next_payload) + ]) + + def __call__(self, _source: str, operation: str, _has_error: bool, + payload: INTERSECT_JSON_VALUE) -> IntersectClientCallback: + print(operation) + if _has_error: + print('============ERROR==============', file=sys.stderr) + print(operation, file=sys.stderr) + print(payload, file=sys.stderr) + raise IntersectCallbackError(operation, payload) + + if operation == 'dial.initialize_workflow': + self.workflow_id = payload + print("\n", "--"*20, "\n") + self.niter += 1 + return self.callback('dial.get_next_point') + + elif operation == 'dial.get_next_point': + self.handle_next_points(payload) + return self.callback('dial.update_workflow_with_data') + + elif operation == 'dial.update_workflow_with_data': + self.niter += 1 + return self.callback('dial.get_next_point') + + else: + raise IntersectCallbackError(operation, payload) + + def handle_surrogate_values(self, payload): + self.variance_grid = np.array(payload[1]).reshape( + (self.meshgrid_size,) * self.num_dims) + self.mean_grid = np.array(payload[0]).reshape( + (self.meshgrid_size,) * self.num_dims) + + if self.niter > self.max_iter: + raise IntersectCallbackEnd() + + def handle_next_points(self, payload): + print(payload) + self.x_next = [payload[0]] + self.above_TN = payload[1] + self.last_idx = payload[2] + + print(f'Running simulation at ({self.x_next}): ', end='', flush=True) + y = peak_val_at_T(*self.x_next) + print(f'{y:.3f}') + print(f'Adding ({self.x_next}, {y}) to dataset') + + self.dataset_x.append(self.x_next) + self.dataset_y.append(y) + + optpos = np.argmax(self.dataset_y) + y_opt = self.dataset_y[optpos] + optimal_coords = self.dataset_x[optpos] + coord_str = ', '.join([f'{coord:.2f}' for coord in optimal_coords]) + print(f'Optimal simulated datapoint at ({coord_str}), y={y_opt:.3f}\n') + + if self.plot_results: + plt.figure() + x_plot = np.linspace(0, 600, 300) + plt.plot(x_plot, peak_val_at_T(x_plot)) + plt.plot(self.dataset_x, self.dataset_y, 'o') + plt.savefig('andie_general.png') + + if self.niter >= self.max_iter: + print("Maximum iteration reached") + raise IntersectCallbackEnd() + if self.last_idx == len(self.x_test) - 1: + print("Last grid point reached") + raise IntersectCallbackEnd() + + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description='Automated client for andie_general backend') + parser.add_argument( + '--config', + type=Path, + default=os.environ.get('DIAL_CONFIG_FILE', Path(__file__).parents[1] / 'remote-conf.json'), + ) + args = parser.parse_args() + + try: + with Path(args.config).open('rb') as f: + from_config_file = json.load(f) + except (json.decoder.JSONDecodeError, OSError) as e: + logger.critical('unable to load config file: %s', str(e)) + sys.exit(1) + + active_learning = ActiveLearningOrchestrator( + service_destination=HierarchyConfig( + **from_config_file['intersect-hierarchy'] + ).hierarchy_string('.') + ) + + config = IntersectClientConfig( + initial_message_event_config=active_learning.initialize_workflow_message, + **from_config_file['intersect'], + ) + + client = IntersectClient(config=config, user_callback=active_learning) + default_intersect_lifecycle_loop(client) diff --git a/src/dial_dataclass/dial_dataclass.py b/src/dial_dataclass/dial_dataclass.py index fc0fa56..f5b695f 100644 --- a/src/dial_dataclass/dial_dataclass.py +++ b/src/dial_dataclass/dial_dataclass.py @@ -6,7 +6,7 @@ PositiveIntType = Annotated[int, Field(ge=0)] -_POSSIBLE_BACKENDS = ('sklearn', 'gpax', 'andie', 'sable') +_POSSIBLE_BACKENDS = ('sklearn', 'gpax', 'andie', 'sable', 'andie_general') BackendType = Literal[_POSSIBLE_BACKENDS] diff --git a/src/dial_service/backends/__init__.py b/src/dial_service/backends/__init__.py index d78bb51..6bef782 100644 --- a/src/dial_service/backends/__init__.py +++ b/src/dial_service/backends/__init__.py @@ -10,6 +10,7 @@ 'sable': ('dial_service.backends.sable_backend', 'SABLEBackend'), 'sklearn': ('dial_service.backends.sklearn_backend', 'SklearnBackend'), 'andie': ('dial_service.backends.andie_backend', 'AndieBackend'), + 'andie_general': ('dial_service.backends.andie_backend_general', 'AndieBackend'), } _MODEL = TypeVar('_MODEL') diff --git a/src/dial_service/backends/andie_backend_general.py b/src/dial_service/backends/andie_backend_general.py index 86fed82..36609cc 100644 --- a/src/dial_service/backends/andie_backend_general.py +++ b/src/dial_service/backends/andie_backend_general.py @@ -1,5 +1,7 @@ """NOTE: This file should not be imported in application code except dynamically via the get_backend_module function in __init__.py .""" +import inspect + import numpy as np import scipy as sp from scipy.optimize import OptimizeResult @@ -97,19 +99,43 @@ def get_prior_distributions(data): @staticmethod def get_parametric_model(data): - model_name = data.model_name.lower() + model_name = data.extra_args['model_name'].lower() if model_name not in _MODELS_ANDIE: raise ValueError(f"Unknown startegy {model_name}") return _MODELS_ANDIE[model_name] + @staticmethod + def _make_bumps_callable(model_func, param_names): + """Wrap a dict-style parametric model into a bumps-compatible callable. + + bumps.Curve inspects the function signature to discover fit parameters. + Dict-style models have signature (inputs, model_params: dict), but bumps + requires (x, p1, p2, ...) with explicit named parameters. This builds a + wrapper whose __signature__ bumps can introspect, while delegating calls + to the dict-based interface. + """ + def wrapper(x, **kwargs): + return model_func(x, model_params=kwargs) + + params = [ + inspect.Parameter('x', inspect.Parameter.POSITIONAL_OR_KEYWORD), + *[ + inspect.Parameter(name, inspect.Parameter.POSITIONAL_OR_KEYWORD) + for name in param_names + ], + ] + wrapper.__signature__ = inspect.Signature(params) + return wrapper + @staticmethod def _Second_order_model(X_grid, Y_grid, model_to_be_optimized, guesses, priors): - - M = Curve(model_to_be_optimized, X_grid, Y_grid, + bumps_callable = AndieBackend._make_bumps_callable(model_to_be_optimized, list(guesses.keys())) + + M = Curve(bumps_callable, X_grid, Y_grid, dy = np.sqrt(Y_grid/100), - pars = guesses) + **guesses) for par_name, prior_dist in priors.items(): M.__getattribute__(par_name).dist = prior_dist