diff --git a/bg_net/archiver/skb_naming.py b/bg_net/archiver/skb_naming.py index 75225ad5e9a823177bd18f29ec63cacd27eb24ce..6f7ff08b3ec2b1345f36c23619ee5fcfc0098aae 100644 --- a/bg_net/archiver/skb_naming.py +++ b/bg_net/archiver/skb_naming.py @@ -703,6 +703,20 @@ SKB2_UNITS = { "CGLOPT:SLYT*:Y": "mm", } +def get_sensitivity_units(observable): + """ + Returns a dictionary with units for all sensitivies for the given observable hitrate. + """ + sens_units_dict = { + 'sens_inj_her': get_unit(observable), + 'sens_inj_ler': get_unit(observable), + 'sens_beam_gas_her': '{:s}/nPa/mA'.format(get_unit(observable)), + 'sens_touschek_her': '{:s}$\mu$m$^3$/mA'.format(get_unit(observable)), + 'sens_beam_gas_ler': '{:s}/nPa/mA'.format(get_unit(observable)), + 'sens_touschek_ler': '{:s}$\mu$m$^3$/mA'.format(get_unit(observable)), + 'sens_lumi': '{:s}cm$^2$s$10^{{-30}}$'.format(get_unit(observable)), + } + return sens_units_dict def get_unit(variable): """ diff --git a/bg_net/hitrate.py b/bg_net/hitrate.py index 9aa705bdd35381c9dabb89c70f83b306a8469790..8cac1b8c98cfd8ebbfd7201a9967cd5ce4b20920 100644 --- a/bg_net/hitrate.py +++ b/bg_net/hitrate.py @@ -38,7 +38,7 @@ from sklearn.preprocessing import RobustScaler, MinMaxScaler pd.set_option('mode.chained_assignment', 'raise') -# Component names for background sources for BGNET +# Column names for background components bg_component_names = [ 'pred_inj_her', 'pred_inj_ler', @@ -50,6 +50,17 @@ bg_component_names = [ 'pred_noise', ] +# Column names for background sensitivities +bg_sensitivity_names = [ + 'sens_inj_her', + 'sens_inj_ler', + 'sens_beam_gas_her', + 'sens_beam_gas_ler', + 'sens_touschek_her', + 'sens_touschek_ler', + 'sens_lumi', +] + # HER injection timing PVS her_inj_timing_pvs = [ "CGHINJ_BEAM_GATE_STATUS", @@ -160,7 +171,7 @@ def get_tensors_no_scaling( # rates, injection efficiency and steerings for beam # injection. - # Input tensor ZH for HER injection amplitude model + # Input tensor ZH for HER injection sensitivity model her_inj_pvs_expanded = bgnet.expand_lagged_pvs(her_inj_pvs) ZH = df[her_inj_pvs_expanded].copy() ZH = ZH.fillna(0) @@ -169,7 +180,7 @@ def get_tensors_no_scaling( ZHG = df[her_inj_timing_pvs].copy() ZHG = ZHG.fillna(0) - # Input tensor ZL for LER injection amplitude model + # Input tensor ZL for LER injection sensitivity model ler_inj_pvs_expanded = bgnet.expand_lagged_pvs(ler_inj_pvs) ZL = df[ler_inj_pvs_expanded].copy() ZL = ZL.fillna(0) @@ -535,8 +546,57 @@ def operateWithConstant(input_batch): tiled_constant = K.tile(tf_constant, (batch_size, 1)) return tiled_constant +def get_mlp_model(num_pvs, name, hidden_units, hidden_layers, activation='tanh'): + """ + Returns a MultiLayerPerceptron keras.models.Model. + The output activation function is always softmax. Only the hidden activation is + configurable. + + :num_pvs: int + number of input variable + + :name: str + name of model + + :hidden_units: int + number of hidden units + + :hidden_layers: int + number of hidden layers + + :activation: str + name of hidden activation function. + """ + inputs = Input(shape=num_pvs) + + output = Dense( + hidden_units, + use_bias=True, + activation=activation, + )(inputs) + + for _ in range(hidden_layers-1): + output = Dense( + hidden_units, + use_bias=True, + activation=activation, + )(output) + + output = Dense( + 1, + use_bias=True, + activation='softplus', + )(output) + + mlp_model = Model( + inputs=inputs, + outputs=output, + name=name + ) + + return mlp_model -def get_gated_injection_model(num_injection, name, hidden_units, hidden_layers): +def get_gated_injection_model(num_injection, name, hidden_units, hidden_layers, num_beam_gate_pvs): """ Returns gated keras.Model for predicting background from injected bunches @@ -544,40 +604,31 @@ def get_gated_injection_model(num_injection, name, hidden_units, hidden_layers): """ inj_input = Input(shape=num_injection) - inj_gate = Input(shape=9) + inj_gate = Input(shape=num_beam_gate_pvs) constant = Lambda(operateWithConstant)(inj_gate) attention = Dense( - 9, + num_beam_gate_pvs, use_bias=True, activation='softmax', name=name+'_gate_attention', )(constant) - inj_output = Dense( - hidden_units, - use_bias=True, - activation='tanh', - )(inj_input) - - for _ in range(hidden_layers): - inj_output = Dense( - hidden_units, - use_bias=True, - activation='tanh', - )(inj_output) + # Compute a predicted gate using attention + gate = Dot(axes=1, name=name+'_pred_gate')([inj_gate, attention]) - inj_output = Dense( - 1, - use_bias=True, - activation='softplus', - )(inj_output) + inj_sensitivity_model = get_mlp_model( + num_pvs=num_injection, + name=name+'_sensitivity', + hidden_units=hidden_units, + hidden_layers=hidden_layers, + ) - # Use attention to attend to the best gate status pv - gate = Dot(axes=1, name=name+'_pred_gate')([inj_gate, attention]) + # Compute sensitivity of injection background + inj_output = inj_sensitivity_model(inj_input) - # Modulate injection background with gate status + # Modulate injection background rate with gate status inj_output = Dot(axes=1)([inj_output, gate]) inj_model = Model( @@ -594,24 +645,15 @@ def get_storage_model(num_pvs, name, hidden_units, hidden_layers): sens_input = Input(shape=num_pvs) feature_input = Input(shape=1) - sens_output = Dense( - hidden_units, - use_bias=True, - activation='tanh', - )(sens_input) - - for _ in range(hidden_layers): - sens_output = Dense( - hidden_units, - use_bias=True, - activation='tanh', - )(sens_output) + sensitivity_model = get_mlp_model( + num_pvs=num_pvs, + name=name+'_sensitivity', + hidden_units=hidden_units, + hidden_layers=hidden_layers, + ) - sens_output = Dense( - 1, - use_bias=True, - activation='softplus', - )(sens_output) + # Model for sensitvity to storage background + sens_output = sensitivity_model(sens_input) storage_output = Dot(axes=1)([sens_output, feature_input]) @@ -681,7 +723,7 @@ def get_beam_lumi_model( sensitivity_pvs_ler_touschek, useDetNoise=False, sensitivity_net_hidden_units=8, - sensitivity_net_hidden_layers=1, + sensitivity_net_hidden_layers=2, const_her_beamgas=None, const_her_touschek=None, const_ler_beamgas=None, @@ -807,9 +849,9 @@ def get_model( obs, useDetNoise=False, sensitivity_net_hidden_units=8, - sensitivity_net_hidden_layers=1, + sensitivity_net_hidden_layers=2, injection_net_hidden_units=16, - injection_net_hidden_layers=1, + injection_net_hidden_layers=2, const_her_beamgas=None, const_her_touschek=None, const_ler_beamgas=None, @@ -855,12 +897,16 @@ def get_model( num_injection_her = len(bgnet.expand_lagged_pvs(injection_pvs_her)) num_injection_ler = len(bgnet.expand_lagged_pvs(injection_pvs_ler)) + num_beam_gate_pvs_her = len(her_inj_timing_pvs) + num_beam_gate_pvs_ler = len(ler_inj_timing_pvs) + # Create model for predicting her injection rate inj_her_model = get_gated_injection_model( num_injection_her, name="injection_her", hidden_units=injection_net_hidden_units, hidden_layers=injection_net_hidden_layers, + num_beam_gate_pvs=num_beam_gate_pvs_her, ) # Create model for predicting ler injection rate @@ -869,6 +915,7 @@ def get_model( name="injection_ler", hidden_units=injection_net_hidden_units, hidden_layers=injection_net_hidden_layers, + num_beam_gate_pvs=num_beam_gate_pvs_ler, ) # Create bgnet model with all parts @@ -879,8 +926,8 @@ def get_model( sens_ler_touschek_input = Input(shape=len(sensitivity_pvs_ler_touschek), name='input_sens_ler_touschek') inj_her_input = Input(shape=num_injection_her, name='input_injection_her') inj_ler_input = Input(shape=num_injection_ler, name='input_injection_ler') - inj_her_gate = Input(shape=9, name='input_injection_her_gate') - inj_ler_gate = Input(shape=9, name='input_injection_ler_gate') + inj_her_gate = Input(shape=num_beam_gate_pvs_her, name='input_injection_her_gate') + inj_ler_gate = Input(shape=num_beam_gate_pvs_ler, name='input_injection_ler_gate') y_bl = bl_model([bl_input, sens_her_beamgas_input, sens_her_touschek_input, sens_ler_beamgas_input, sens_ler_touschek_input]) y_ih = inj_her_model([inj_her_input, inj_her_gate]) @@ -1031,65 +1078,238 @@ def get_sample_weights(X, weights_single_beams=0.0, weights_both_beams=0.0): # CONVINIENCE FUNCTIONS -def get_predictions(model, x_scaled, scalers_x, scaler_y): +def get_predictions(model, inputs, scaler_y, index): + """ + Returns pandas.DataFrame with predictions for all background components using model. + + Parameters + ---------- + model : keras.models.Model + Trained BGNet model for hitrate of Belle II subdetector - scaler_X = scalers_x[0] - X, SHB, SHT, SLB, SLT, ZH, ZHG, ZL, ZLG = x_scaled + inputs : tuple(np.ndarray) + Tuple of preprocessed input tensors for BGNet model - # For masking unscaled inputs - mask = { - 'pred_beam_gas_her': np.array([1, 0, 0, 0, 0]), - 'pred_touschek_her': np.array([0, 1, 0, 0, 0]), - 'pred_beam_gas_ler': np.array([0, 0, 1, 0, 0]), - 'pred_touschek_ler': np.array([0, 0, 0, 1, 0]), - 'pred_lumi': np.array([0, 0, 0, 0, 1]), - } + scaler_y: sklearn.preprocessing.RobustScaler + Scaler for target detector hitrate + index : DateTimeIndex + Index for output DataFrame + + Returns + ------- + output: pandas.DataFrame + DataFrame with predicted components of background + """ # Dictionary with all output components output = {} - # Compute total predicted rate - yhat = model.predict([X, SHB, SHT, SLB, SLT, ZH, ZHG, ZL, ZLG]) + # Unpack scaled (preprocessed) inputs to model + X, SHB, SHT, SLB, SLT, ZH, ZHG, ZL, ZLG = inputs + # Compute total predicted background rate + yhat = model.predict([X, SHB, SHT, SLB, SLT, ZH, ZHG, ZL, ZLG]) output['pred'] = scaler_y.inverse_transform(yhat) + # Compute predicted background from HER injections yhat = model.get_layer('injection_her').predict([ZH, ZHG]) output['pred_inj_her'] = scaler_y.inverse_transform(yhat) + # Compute predicted background from LER injections yhat = model.get_layer('injection_ler').predict([ZL, ZLG]) output['pred_inj_ler'] = scaler_y.inverse_transform(yhat) - # Compute components of total rate as well - Xraw = scaler_X.inverse_transform(X) + # Compute predicted background from HER beam gas + yhat = model.get_layer('beam_lumi').get_layer('her_beam_gas').predict([SHB,X[:,0]]) + output['pred_beam_gas_her'] = scaler_y.inverse_transform(yhat) + + # Compute predicted background from HER Touschek + yhat = model.get_layer('beam_lumi').get_layer('her_touschek').predict([SHT,X[:,1]]) + output['pred_touschek_her'] = scaler_y.inverse_transform(yhat) + + # Compute predicted background from LER beam gas + yhat = model.get_layer('beam_lumi').get_layer('ler_beam_gas').predict([SLB,X[:,2]]) + output['pred_beam_gas_ler'] = scaler_y.inverse_transform(yhat) + + # Compute predicted background from LER Touschek + yhat = model.get_layer('beam_lumi').get_layer('ler_touschek').predict([SLT,X[:,3]]) + output['pred_touschek_ler'] = scaler_y.inverse_transform(yhat) - Xtemp = scaler_X.transform(np.zeros_like(Xraw)) - yhat = model.get_layer('beam_lumi').predict([Xtemp, SHB, SHT, SLB, SLT]) + # Compute predicted background from detector noise (implemented as bias of lumi model) + yhat = model.get_layer('beam_lumi').get_layer('lumi').predict(np.zeros_like(X[:,4])) output['pred_noise'] = scaler_y.inverse_transform(yhat) - for k in mask.keys(): - Xtemp = scaler_X.transform(mask[k]*Xraw) - yhat = model.get_layer('beam_lumi').predict([Xtemp, SHB, SHT, SLB, SLT]) - output[k] = scaler_y.inverse_transform(yhat) - output['pred_noise'] + # Compute predicted background from Lumi (need to subract detector noise) + yhat = model.get_layer('beam_lumi').get_layer('lumi').predict(X[:,4]) + output['pred_lumi'] = scaler_y.inverse_transform(yhat) - output['pred_noise'] + + for comp in output: + output[comp] = output[comp].squeeze(-1) + + output = pd.DataFrame(data=output, index=index) return output +def get_sensitivity_predictions(model, inputs, scalers_x, scaler_y, index): + """ + Returns pandas.DataFrame with sensitiviy predictions. Sensitivities are returned + in physical units. + + Parameters + ---------- + model : keras.models.Model + Trained BGNet model for hitrate of Belle II subdetector + + inputs : tuple(np.ndarray) + Tuple of preprocessed input tensors for BGNet model + + scalers_x: sklearn.preprocessing.RobustScaler + Scaler for inputs + + scaler_y: sklearn.preprocessing.RobustScaler + Scaler for target detector hitrate + + index : DateTimeIndex + Index for output DataFrame + + Returns + ------- + output: pandas.DataFrame + DataFrame with background sensitivities + """ + # Dictionary with all output components + output = {} + + # Unpack scaled (preprocessed) inputs to model + X, SHB, SHT, SLB, SLT, ZH, ZHG, ZL, ZLG = inputs + + sub_model = model.get_layer('injection_her').get_layer('injection_her_sensitivity') + yhat = sub_model.predict(ZH) + output['sens_inj_her'] = scaler_y.inverse_transform(yhat)*get_X_scale_factor(scalers_x,'injection_her') + + sub_model = model.get_layer('injection_ler').get_layer('injection_ler_sensitivity') + yhat = sub_model.predict(ZL) + output['sens_inj_ler'] = scaler_y.inverse_transform(yhat)*get_X_scale_factor(scalers_x,'injection_ler') + + sub_model = model.get_layer("beam_lumi").get_layer('her_beam_gas').get_layer('her_beam_gas_sensitivity') + yhat = sub_model.predict(SHB) + output['sens_beam_gas_her'] = scaler_y.inverse_transform(yhat)*get_X_scale_factor(scalers_x,'her_beam_gas') + + sub_model = model.get_layer("beam_lumi").get_layer('her_touschek').get_layer('her_touschek_sensitivity') + yhat = sub_model.predict(SHT) + output['sens_touschek_her'] = scaler_y.inverse_transform(yhat)*get_X_scale_factor(scalers_x,'her_touschek') + + sub_model = model.get_layer("beam_lumi").get_layer('ler_beam_gas').get_layer('ler_beam_gas_sensitivity') + yhat = sub_model.predict(SLB) + output['sens_beam_gas_ler'] = scaler_y.inverse_transform(yhat)*get_X_scale_factor(scalers_x, 'ler_beam_gas') + + sub_model = model.get_layer("beam_lumi").get_layer('ler_touschek').get_layer('ler_touschek_sensitivity') + yhat = sub_model.predict(SLT) + output['sens_touschek_ler'] = scaler_y.inverse_transform(yhat)*get_X_scale_factor(scalers_x,'her_touschek') + + sub_model = model.get_layer('beam_lumi').get_layer('lumi') + yhat = sub_model.predict(np.ones_like(X[:,4])) - sub_model.predict(np.zeros_like(X[:,4])) + output['sens_lumi'] = scaler_y.inverse_transform(yhat)*get_X_scale_factor(scalers_x,'lumi') + + for comp in output: + output[comp] = output[comp].squeeze(-1) + + output = pd.DataFrame(data=output, index=index) + + return output + +def get_predicted_gate_attention_weights(model): + """ + Returns tuple of np.ndarray with learned gate attention weights for HER and LER ring for BGNet model. + Note that attention weights are static. They only depend on the given Belle II subdetector whose hitrate + was used as training target. + + Parameters + ---------- + model : keras.models.Model + Trained BGNet model for hitrate of Belle II subdetector + + Returns + ------- + output: tuple(np.ndarray) + tuple of attention weights for HER and LER rings. + """ + + model_attention_her_gate = Model( + inputs=model.get_layer('injection_her').input, + outputs=model.get_layer('injection_her').get_layer('injection_her_gate_attention').output + ) + + num_her_inputs = model.get_layer('injection_her').input_shape[0][1] + num_her_gate_inputs = model.get_layer('injection_her').input_shape[1][1] + + ZH = np.zeros((1,num_her_inputs)) + ZHG = np.zeros((1,num_her_gate_inputs)) + w_her = model_attention_her_gate([ZH, ZHG]) + + model_attention_ler_gate = Model( + inputs=model.get_layer('injection_ler').input, + outputs=model.get_layer('injection_ler').get_layer('injection_ler_gate_attention').output + ) + + num_ler_inputs = model.get_layer('injection_ler').input_shape[0][1] + num_ler_gate_inputs = model.get_layer('injection_ler').input_shape[1][1] + + ZL = np.zeros((1,num_ler_inputs)) + ZLG = np.zeros((1,num_ler_gate_inputs)) + w_ler = model_attention_ler_gate([ZL, ZLG]) + + return w_her.numpy().squeeze(), w_ler.numpy().squeeze() + + +def get_X_scale_factor(scalers_x, name): + """ + Returns scale factor for features engineered heuristics in input tensor X. + Since each column in X uniquely belongs to a background component, we can + use the name of the bg components to query the scale factor. + + :scalers_x: sklearn.preprocessing.RobustScaler + Scaler of input tensor X + + :name: str + Name of bg component + """ + if name == 'her_beam_gas': + return scalers_x[0].transform([[1, 1, 1, 1, 1]])[0][0] + elif name == 'her_touschek': + return scalers_x[0].transform([[1, 1, 1, 1, 1]])[0][1] + elif name == 'ler_beam_gas': + return scalers_x[0].transform([[1, 1, 1, 1, 1]])[0][2] + elif name == 'ler_touschek': + return scalers_x[0].transform([[1, 1, 1, 1, 1]])[0][3] + elif name == 'lumi': + return scalers_x[0].transform([[1, 1, 1, 1, 1]])[0][4] + else: + return 1.0 + def get_predictions_from_data( - df, obs, model, config, repair=False, cleanup=False, + df, obs, model, config, repair=False, cleanup=False, extended=False, ): """ Return pandas.DataFrame with predictions for observable obs on input pandas.DataFrame data using trained model and scalers. - :param df: pandas.DataFrame with data from archive - :param obs: Name of observable which BGNet is predicting - :param model: The trained model, instance of tf.keras.Model - :param repair: Bool flag indicating if repair of df is needed - :param cleanup: Bool flag indicating if cleaning of df is needed + :param df: pandas.DataFrame + data from archive + :param obs: str + Name of observable which BGNet is predicting + :param model: tf.keras.models.Model + The trained model + :param repair: bool + Flag indicating if repair of df is needed + :param cleanup: bool + Flag indicating if cleaning of df is needed + :param extended: bool + Flag indicating if computing of sensitivities is needed """ - # ToDo data = df.copy() # Preprocess data @@ -1122,35 +1342,20 @@ def get_predictions_from_data( x_scaled = bgnet.scale_x(x, scalers_x) - decomp_bgnet = get_predictions( - model=model, x_scaled=x_scaled, scalers_x=scalers_x, scaler_y=scaler_y + bgnet_df = get_predictions( + model=model, inputs=x_scaled, scaler_y=scaler_y, index=data.index, ) - for comp in decomp_bgnet: - decomp_bgnet[comp] = decomp_bgnet[comp].squeeze(-1) - - decomp_bgnet = pd.DataFrame(data=decomp_bgnet, index=data.index) - - return decomp_bgnet - - -def get_sensitivity_predictions(pred_df, df): - """ - Returns DataFrame with sensitivities to storage and lumonosity backgrounds. + if extended: + sens_df = get_sensitivity_predictions( + model=model, inputs=x_scaled, scalers_x=scalers_x, scaler_y=scaler_y, index=data.index + ) + bgnet_df = pd.concat([bgnet_df, sens_df], axis=1) - Note: Storage sensitivities for HER (LER) will be NaN in case there is no beam in HER (LER) ring - Note: Luminosity sensitivty will be NaN in case there is no colliding beams - """ - s0 = pred_df['pred_beam_gas_her'].divide(prep_utils.get_storage_beam_feature_0(df)) - s1 = pred_df['pred_touschek_her'].divide(prep_utils.get_storage_beam_feature_1(df)) - s2 = pred_df['pred_beam_gas_ler'].divide(prep_utils.get_storage_beam_feature_2(df)) - s3 = pred_df['pred_touschek_ler'].divide(prep_utils.get_storage_beam_feature_3(df)) - s4 = pred_df['pred_lumi'].divide(prep_utils.get_storage_beam_feature_4(df)) - sens_df = pd.concat([s0, s1, s2, s3, s4], axis=1) - return sens_df + return bgnet_df -def get_scaled_sensitivty(sensitivity, model, scalers_x, scaler_y): +def get_scaled_sensitivity(sensitivity, model, scalers_x, scaler_y): """ Returns scaled sensitivity for model given bgnet scalers. @@ -1179,67 +1384,55 @@ map_sub_models_to_decomp = { 'ler_touschek': 'pred_touschek_ler', } +map_sub_models_to_sensitivity = { + 'injection_her': 'sens_inj_her', + 'injection_ler': 'sens_inj_ler', + 'her_beam_gas': 'sens_beam_gas_her', + 'her_touschek': 'sens_touschek_her', + 'ler_beam_gas': 'sens_beam_gas_ler', + 'ler_touschek': 'sens_touschek_ler', +} + index_map = { - 'injection_her': (5, 6), - 'injection_ler': (7, 8), - 'her_beam_gas': (1,), - 'her_touschek': (2,), - 'ler_beam_gas': (3,), - 'ler_touschek': (4,), + 'injection_her': 5, + 'injection_ler': 7, + 'her_beam_gas': 1, + 'her_touschek': 2, + 'ler_beam_gas': 3, + 'ler_touschek': 4, } feature_map = { - 'injection_her': ('her_inj_pvs',), - 'injection_ler': ('ler_inj_pvs',), - 'her_beam_gas': ('her_beamgas_sens_pvs', 'F_B_HER'), - 'her_touschek': ('her_touschek_sens_pvs', 'F_T_HER'), - 'ler_beam_gas': ('ler_beamgas_sens_pvs', 'F_B_LER'), - 'ler_touschek': ('ler_touschek_sens_pvs', 'F_T_LER'), + 'injection_her': 'her_inj_pvs', + 'injection_ler': 'ler_inj_pvs', + 'her_beam_gas': 'her_beamgas_sens_pvs', + 'her_touschek': 'her_touschek_sens_pvs', + 'ler_beam_gas': 'ler_beamgas_sens_pvs', + 'ler_touschek': 'ler_touschek_sens_pvs', } def get_interpretable_inputs(x, name): + i = index_map[name] + return x[i] - if name in ['injection_her', 'injection_ler']: - i1, i2 = index_map[name][0], index_map[name][1] - interpret_x = np.concatenate((x[i1], x[i2]), axis=1) - elif name in ["her_beam_gas", "ler_beam_gas", "her_touschek", "ler_touschek"]: - i = index_map[name][0] - feature = np.expand_dims(x[0][:, i], axis=1) - interpret_x = np.concatenate((x[i], feature), axis=1) - else: - raise ValueError("Unkown sub model name {}".format(name)) - - return interpret_x - -def get_intepretable_model(bgnet_model, num_inputs, name='injection_her'): - - # Clone bgnet injection model and set weights +def get_intepretable_model(model, num_inputs, name): + # Clone submodel and copy weights if name in ['injection_her', 'injection_ler']: - cloned_model = keras.models.clone_model(bgnet_model.get_layer(name)) - cloned_model.set_weights(bgnet_model.get_layer(name).get_weights()) + cloned_model = keras.models.clone_model(model.get_layer(name).get_layer(name+'_sensitivity')) + cloned_model.set_weights(model.get_layer(name).get_layer(name+'_sensitivity').get_weights()) elif name in ["her_beam_gas", "ler_beam_gas", "her_touschek", "ler_touschek"]: - cloned_model = keras.models.clone_model(bgnet_model.get_layer("beam_lumi").get_layer(name)) - cloned_model.set_weights(bgnet_model.get_layer("beam_lumi").get_layer(name).get_weights()) + cloned_model = keras.models.clone_model(model.get_layer("beam_lumi").get_layer(name).get_layer(name+'_sensitivity')) + cloned_model.set_weights(model.get_layer("beam_lumi").get_layer(name).get_layer(name+'_sensitivity').get_weights()) else: raise ValueError("Unkown sub model name {}".format(name)) # Input layer for interpretable model interpretable_inputs = keras.layers.Input(shape=num_inputs) - # If the cloned_model has more than 1 input feature, add adapter layers - if name in ['injection_her', 'injection_ler']: - timing_inputs = len(her_inj_timing_pvs) - x1 = keras.layers.Lambda(lambda x: x[:, 0:num_inputs-timing_inputs])(interpretable_inputs) - x2 = keras.layers.Lambda(lambda x: x[:, num_inputs-timing_inputs:num_inputs])(interpretable_inputs) - y = cloned_model([x1, x2]) - elif name in ["her_beam_gas", "ler_beam_gas", "her_touschek", "ler_touschek"]: - x1 = keras.layers.Lambda(lambda x: x[:, 0:num_inputs-1])(interpretable_inputs) - x2 = keras.layers.Lambda(lambda x: x[:, num_inputs-1:num_inputs])(interpretable_inputs) - y = cloned_model([x1, x2]) - else: - raise ValueError("Unkown sub model name {}".format(name)) + # Output layer for interpretable model + y = cloned_model(interpretable_inputs) # Create interpretable model interpret_model = keras.models.Model( @@ -1254,17 +1447,11 @@ def get_intepretable_model(bgnet_model, num_inputs, name='injection_her'): def get_interpretable_feature_names(config, name): if name in ['injection_her', 'injection_ler']: - feature_group = feature_map[name][0] + feature_group = feature_map[name] feature_names = bgnet.expand_lagged_pvs(config['preparation'][feature_group]) - if 'her' in name: - feature_names = feature_names + her_inj_timing_pvs - else: - feature_names = feature_names + ler_inj_timing_pvs elif name in ["her_beam_gas", "ler_beam_gas", "her_touschek", "ler_touschek"]: - feature_group = feature_map[name][0] - gate_feature = feature_map[name][1] + feature_group = feature_map[name] feature_names = config['preparation'][feature_group] - feature_names.append(gate_feature) else: raise ValueError("Unkown sub model name {}".format(name)) diff --git a/bg_net/make_explanations.py b/bg_net/make_explanations.py index 82253b1607d37485dc5d7ea843726cc61366340f..0e413ef472d4e247c342a1317a9fe0a38bea1396 100644 --- a/bg_net/make_explanations.py +++ b/bg_net/make_explanations.py @@ -33,51 +33,10 @@ import matplotlib.pyplot as plt import bg_net.hitrate as hitrate import bg_net.bgnet as bgnet +from bg_net.report_utils import gate_attention_plots pd.set_option('mode.chained_assignment', 'raise') - -def get_best_gate_status(attributions, feature_names,): - mean_abs_attr = np.mean(np.abs(attributions), axis=0) - max_order = np.argsort(mean_abs_attr) - index_order = max_order[::-1] - names_order = [feature_names[i] for i in index_order] - - for i, name in enumerate(names_order): - if 'BEAM_GATE_STATUS' in name: - if 'shift' in name: - delay = int(name.split('_')[-1]) - else: - delay = 0 - - return index_order[i], name, delay - - return None, None, None - - -def get_controlables(features): - """ - Returns a list of controlable features found among features. - """ - - controlables = [] - - for feature in features: - - if 'GATE' in feature: - continue - - if 'shift' in feature: - continue - - controlables.append(feature) - - if 'ECL_LUM_MON_ACC' not in controlables: - controlables.append('ECL_LUM_MON_ACC') - - return controlables - - def find_configs(local_dir): """ Returns a list of configs for validations runs of BGNet models in local_dir. @@ -109,11 +68,10 @@ def make_attribution_plots(sub_model_name, x_train, x_test, config, plot_path, p logging.info('Loading saved models from ' + config['training']['model_path']) # Load the model and metavariables model = tf.keras.models.load_model(config['training']['model_path'] + '/{}'.format(config['obs'])) - except Exception as e: logging.info('Error with {}. Reason: {}'.format(config['obs'], repr(e))) logging.info('Skipping explanation for {}'.format(config['obs'])) - return + return None, None # Create an interpretable model for sub_model interpret_model = hitrate.get_intepretable_model(model, num_inputs=len(feature_names), name=sub_model_name) @@ -149,6 +107,9 @@ def make_attribution_plots(sub_model_name, x_train, x_test, config, plot_path, p verbose=False, ) + # Scale attributions to physical units + attributions = scaler_y.inverse_transform(attributions)*hitrate.get_X_scale_factor(scalers_x,sub_model_name) + if len(feature_names) < plot_top_k: plot_top_k = len(feature_names) @@ -164,45 +125,13 @@ def make_attribution_plots(sub_model_name, x_train, x_test, config, plot_path, p # Close figure created by summary_plot plt.close('all') - # Making extra plots to understand delays of injections and mitigation of injections bg - if 'injection' in sub_model_name: - - # Extract month from data path for plot labels - month = '_'.join(config['preparation']['data_path'].split('_')[-2:]) - - Path(plot_path+'/delays_{}'.format(sub_model_name)).mkdir(parents=True, exist_ok=True) - - # Get name of best working GATE_STATUS variable. Also extracts delay value and index - best_gate_index, best_gate_name, best_delay = get_best_gate_status(attributions, feature_names,) - - # Predict injection background. Output is still scaled and has not yet physical units. - interpret_y = interpret_model.predict(interpret_x_test_scaled) + # Plot attention weights for Beam Gate Status PVs + w_her, w_ler = hitrate.get_predicted_gate_attention_weights(model) + gate_attention_plots(w_her, w_ler) + plt.savefig(os.path.join(plot_path, "attention_weights.pdf"), bbox_inches='tight') - # Scanning the delays, marking the one with highest attributions - # Here we use the predicted injection rate on y axis. GATE_STATUS off should give zeros when shift is correct. - - gate_status_pvs = [name for name in feature_names if 'BEAM_GATE_STATUS' in name] - for gate_status_pv in gate_status_pvs: - - gate_index = feature_names.index(gate_status_pv) - - plt.scatter( - x=np.arange(interpret_y.shape[0]), - y=interpret_y.squeeze(), - c=interpret_x_test[:, gate_index], - ) - - if gate_index == best_gate_index: - plt.title('Color coded by BEST {}'.format(feature_names[gate_index])) - else: - plt.title('Color coded by {} \n {} {}'.format(feature_names[gate_index], config['obs'], month)) - plt.xlabel('Data index') - plt.ylabel('Predicted injection rate / a.u.') - - plt.savefig(plot_path+'/delays_{}/{}.png'.format(sub_model_name, feature_names[gate_index])) - - # Close figure created by summary_plot - plt.close('all') + # Close figure created by gate_attention_plots + plt.close('all') return attributions, feature_names diff --git a/bg_net/make_prediction.py b/bg_net/make_prediction.py index 285f08632deea526179f41f221d40507410c68a9..7eaf4d4b20fec49fd366594a13059b54419db8ad 100644 --- a/bg_net/make_prediction.py +++ b/bg_net/make_prediction.py @@ -128,7 +128,7 @@ def predict_save(config, write_csv=False): # Process data bgnet_df = hitrate.get_predictions_from_data( - data, obs, model, model_config, repair=True, + data, obs, model, model_config, repair=True, extended=True, ) # Save hitrate prediction diff --git a/bg_net/make_validation.py b/bg_net/make_validation.py index 3e1b19c31a4d1bf0f31891e4a414537411a3545c..857efa8747507dffeb30f562291b8e56e0f27162 100644 --- a/bg_net/make_validation.py +++ b/bg_net/make_validation.py @@ -37,118 +37,63 @@ import bg_net.hitrate as hitrate pd.set_option('mode.chained_assignment', 'raise') +sens_labels = [ + 'B<sub>HER</sub>', + 'T<sub>HER</sub>', + 'B<sub>LER</sub>', + 'T<sub>LER</sub>', + 'C<sub>L</sub>', +] + +unit_obs = "unit(O)" + +sens_units = [ + '{:s}/nPa/mA'.format(unit_obs), + '{:s}μm<sup>3</sup>/mA'.format(unit_obs), + '{:s}/nPa/mA'.format(unit_obs), + '{:s}μm<sup>3</sup>/mA'.format(unit_obs), + '{:s}cm<sup>2</sup>s10<sup>-30</sup>'.format(unit_obs), +] + +sens_columns = [ + 'sens_beam_gas_her', + 'sens_touschek_her', + 'sens_beam_gas_ler', + 'sens_touschek_ler', + 'sens_lumi', +] + def get_decay_type(df): """ Returns integer decay type of dataframe df with raw PV data. + + Types ars: LER single beam decay (1), LER single beam decay (2), + Both beams decay (3) and any other decay (4). """ if len(df) < 40: - return 4 # Overflow option + return 4 delta_ler = df["SKB_LER_current"].iloc[20] - df["SKB_LER_current"].iloc[-10] delta_her = df["SKB_HER_current"].iloc[20] - df["SKB_HER_current"].iloc[-10] if ((df["SKB_HER_current"] < 2) & (df["SKB_LER_current"] > 50)).all(): - decay_type = 1 # LER single beam decay + decay_type = 1 elif ((df["SKB_LER_current"] < 2) & (df["SKB_HER_current"] > 50)).all(): - decay_type = 2 # HER single beam decay + decay_type = 2 elif delta_ler > 3 and delta_her > 3: - decay_type = 3 # Both beams decay + decay_type = 3 else: - decay_type = 4 # Overflow option + decay_type = 4 return decay_type -def print_master_formula(config, plot_path): - import matplotlib.pyplot as plt - - plt.figure(figsize=(10, 8), dpi=180) - - d = 0.065 - fontsize = 11 - plt.text( - 0.01, - 0.9, - r'$O = T_{H}(x_{1}) \cdot F_{T,H} + B_{H}(x_{2}) \cdot F_{B,H} + T_{L}(x_{3}) \cdot F_{T,L} + B_{L}(x_{4}) \cdot F_{B,L} + C_{L} \cdot L + D + INJ_{H}(x_{5}) + INJ_{L}(x_{6})$', - fontsize=fontsize) - plt.text(0.01, 0.9-1*d, r'$F_{{B,H}}=3 \cdot P_{1,H} \cdot I_{H}^{2} + P_{0,H} \cdot I_{H}$', fontsize=fontsize) - plt.text( - 0.01, - 0.9 - - 2 * - d, - r'$F_{T,H}=\frac{I_{H}^{2}}{\sigma_{X,H} \cdot \sigma_{Y,H} \cdot (\sigma_{Z,0,H} + \sigma_{Z,1,H} \cdot I_{H}/n_{b,H}) \cdot n_{b,H}}$', - fontsize=fontsize) - plt.text(0.01, 0.9-3*d, r'$F_{{B,L}}=3 \cdot P_{1,L} \cdot I_{L}^{2} + P_{0,L} \cdot I_{L}$', fontsize=fontsize) - plt.text( - 0.01, - 0.9 - - 4 * - d, - r'$F_{T,L}=\frac{I_{L}^{2}}{\sigma_{X,L} \cdot \sigma_{Y,L} \cdot (\sigma_{Z,0,L} + \sigma_{Z,1,L} \cdot I_{L}/n_{b,L}) \cdot n_{b,L}}$', - fontsize=fontsize) - - plt.text( - 0.01, - 0.9-5*d, - r'$\sigma_{{Z,0,H}} = {:.3g}$ $\mu$m'.format( - config['preparation']['sigmaZ_params']['HER_z0']), - fontsize=fontsize) - plt.text(0.01, - 0.9-6*d, - r'$\sigma_{{Z,1,H}} = {:.3g}$ $\mu$m/mA'.format(config['preparation']['sigmaZ_params']['HER_z1']), - fontsize=fontsize) - plt.text( - 0.01, - 0.9-7*d, - r'$\sigma_{{Z,0,L}} = {:.3g}$ $\mu$m'.format( - config['preparation']['sigmaZ_params']['LER_z0']), - fontsize=fontsize) - plt.text(0.01, - 0.9-8*d, - r'$\sigma_{{Z,1,L}} = {:.3g}$ $\mu$m/mA'.format(config['preparation']['sigmaZ_params']['LER_z1']), - fontsize=fontsize) - - fig = plt.gca() - fig.axes.get_xaxis().set_visible(False) - fig.axes.get_yaxis().set_visible(False) - plt.savefig(plot_path + '/hitrate_formula.png') - - return fig - - def make_sensitivity_overview_plots(obs, df, plot_path, show_figs=False, write_html=True, write_image=True): - sens_list = [ - 'B_HER', - 'T_HER', - 'B_LER', - 'T_LER', - 'C_L', - ] - - sens_title = [ - 'B<sub>HER</sub>', - 'T<sub>HER</sub>', - 'B<sub>LER</sub>', - 'T<sub>LER</sub>', - 'C<sub>L</sub>', - ] - unit_obs = "unit(O)" - units = [ - '{:s}/nPa/mA'.format(unit_obs), - '{:s}μm<sup>3</sup>/mA'.format(unit_obs), - '{:s}/nPa/mA'.format(unit_obs), - '{:s}μm<sup>3</sup>/mA'.format(unit_obs), - '{:s}cm<sup>2</sup>s10<sup>-30</sup>'.format(unit_obs), - ] - - for i, sens_label in enumerate(sens_list): + for i in range(len(sens_columns)): fig = go.Figure() - fig.add_trace(go.Scatter(x=df.index, y=df[i], - name=sens_title[i],), - ) + fig.add_trace(go.Scatter(x=df.index, y=df[sens_columns[i]], name=sens_labels[i],),) fig.update_layout( legend_title='', @@ -157,16 +102,16 @@ def make_sensitivity_overview_plots(obs, df, plot_path, show_figs=False, write_h title_font_size=22, xaxis_title='time', xaxis_title_font_size=18, - yaxis_title=sens_title[i] + ' ' + units[i], + yaxis_title=sens_labels[i] + ' ' + sens_units[i], yaxis_title_font_size=18, ) # Fine tune y range of figure, make sure not only rounding errors # are displayed - manual_min_range = (1-1.0e-04)*df[i].mean() - manual_max_range = (1+1.0e-04)*df[i].mean() + manual_min_range = (1-1.0e-04)*df[sens_columns[i]].mean() + manual_max_range = (1+1.0e-04)*df[sens_columns[i]].mean() - if df[i].min() > manual_min_range and df[i].max() < manual_max_range: + if df[sens_columns[i]].min() > manual_min_range and df[sens_columns[i]].max() < manual_max_range: fig.update_yaxes(range=(manual_min_range, manual_max_range)) if show_figs: @@ -174,56 +119,24 @@ def make_sensitivity_overview_plots(obs, df, plot_path, show_figs=False, write_h if write_html: fig.write_html( - plot_path + "/{}.html".format(sens_label) + plot_path + "/{}.html".format(sens_columns[i]) ) if write_image: fig.write_image( - plot_path + "/{}.png".format(sens_label) + plot_path + "/{}.png".format(sens_columns[i]) ) def make_sensitivity_summary_plots(obs, df, plot_path, tag='decay', show_figs=False, write_html=False, write_image=True): - sens_labels = [ - 'B_HER', - 'T_HER', - 'B_LER', - 'T_LER', - 'C_L', - ] - - sens_titles = [ - 'B<sub>HER</sub>', - 'T<sub>HER</sub>', - 'B<sub>LER</sub>', - 'T<sub>LER</sub>', - 'C<sub>L</sub>', - ] - unit_obs = "unit(O)" - sens_units = [ - '{:s}/nPa/mA'.format(unit_obs), - '{:s}μm<sup>3</sup>/mA'.format(unit_obs), - '{:s}/nPa/mA'.format(unit_obs), - '{:s}μm<sup>3</sup>/mA'.format(unit_obs), - '{:s}cm<sup>2</sup>s1e<sup>-30</sup>'.format(unit_obs), - ] - - sens_keys = [ - ('sens_her_beam_gas_mean',), - ('sens_her_touschek_mean',), - ('sens_ler_beam_gas_mean',), - ('sens_ler_touschek_mean',), - ('sens_lumi_mean',), - ] - - for i, sens_label in enumerate(sens_labels): + for i in range(len(sens_columns)): fig = go.Figure([ go.Scatter( - name=sens_titles[i], + name=sens_labels[i], x=df['start'], - y=df[sens_keys[i][0]], + y=df[sens_columns[i]], mode='lines+markers', line=dict(color='rgb(31, 119, 180)'), ), @@ -234,16 +147,16 @@ def make_sensitivity_summary_plots(obs, df, plot_path, tag='decay', show_figs=Fa title_font_size=22, xaxis_title='time', xaxis_title_font_size=18, - yaxis_title=sens_titles[i] + ' ' + sens_units[i], + yaxis_title=sens_labels[i] + ' ' + sens_units[i], yaxis_title_font_size=18, ) # Fine tune y range of figure, make sure not only rounding errors # are displayed - manual_min_range = (1-1.0e-04)*df[sens_keys[i][0]].mean() - manual_max_range = (1+1.0e-04)*df[sens_keys[i][0]].mean() + manual_min_range = (1-1.0e-04)*df[sens_columns[i]].mean() + manual_max_range = (1+1.0e-04)*df[sens_columns[i]].mean() - if df[sens_keys[i][0]].min() > manual_min_range and df[sens_keys[i][0]].max() < manual_max_range: + if df[sens_columns[i]].min() > manual_min_range and df[sens_columns[i]].max() < manual_max_range: fig.update_yaxes(range=(manual_min_range, manual_max_range)) if show_figs: @@ -251,11 +164,11 @@ def make_sensitivity_summary_plots(obs, df, plot_path, tag='decay', show_figs=Fa if write_html: fig.write_html( - plot_path + "/summary_{}_{}_{}.html".format(tag, obs, sens_keys[i][0]) + plot_path + "/summary_{}_{}_{}.html".format(tag, obs, sens_columns[i]) ) if write_image: fig.write_image( - plot_path + "/summary_{}_{}_{}.png".format(tag, obs, sens_keys[i][0]) + plot_path + "/summary_{}_{}_{}.png".format(tag, obs, sens_columns[i]) ) @@ -432,11 +345,7 @@ def make_pressure_plots(data, plot_path, show_figs=False, write_html=False, writ def make_window_plots( - obs, windows, data, pred, plot_path, tag='decay', show_figs=False, write_html=False, sens_df=None, - annotate=True, - ): - - summary = [] + obs, windows, data, pred, plot_path, tag='decay', show_figs=False, write_html=False, annotate=True): summary_labels = [ 'start', @@ -448,18 +357,16 @@ def make_window_plots( 'dI_HER', 'dI_LER', 'decay_type', - 'sens_her_beam_gas_mean', - 'sens_her_touschek_mean', - 'sens_ler_beam_gas_mean', - 'sens_ler_touschek_mean', - 'sens_lumi_mean', 'HER_p0', 'HER_p1', 'LER_p0', 'LER_p1', 'D', + *sens_columns ] + summary = [] + for start, end in windows: df = data[start:end].copy() @@ -468,10 +375,10 @@ def make_window_plots( if len(df) < 40: continue - if sens_df is not None: - sens_mean = list(sens_df[start:end].mean()) + if 'sens_lumi' in pred_df.columns: + sens_mean = list(pred_df[sens_columns][start:end].mean()) else: - sens_mean = [np.NaN, np.NaN, np.NaN, np.NaN, np.NaN] + sens_mean = [np.NaN for column in sens_columns] # Compute statistics for decay decay_type = get_decay_type(df) @@ -491,12 +398,12 @@ def make_window_plots( delta_her, delta_ler, decay_type, - *sens_mean, df["SKB_HER_P0_fit"].mean(), df["SKB_HER_P1_fit"].mean(), df["SKB_LER_P0_fit"].mean(), df["SKB_LER_P1_fit"].mean(), pred_df['pred_noise'].mean(), + *sens_mean, ]) fig = plot_utils.plot_pred_timeseries_only(obs, df, pred_df) @@ -508,15 +415,6 @@ def make_window_plots( ) if annotate: - unit_obs = "unit(O)" - sens_units = [ - '{:s}/nPa/mA'.format("unit(O)"), - '{:s}μm<sup>3</sup>/mA'.format(unit_obs), - '{:s}/nPa/mA'.format(unit_obs), - '{:s}μm<sup>3</sup>/mA'.format(unit_obs), - '10<sup>-30</sup>{:s}cm<sup>2</sup>/s'.format(unit_obs), - ] - fig.update_layout( annotations=[ go.layout.Annotation( @@ -607,13 +505,13 @@ def create_validation_config(config_pred, local_dir): return config -def add_to_resampled_data(resampled_data, data, bgnet_df, sens_df, mask, resample_period): +def add_to_resampled_data(resampled_data, data, bgnet_df, mask, resample_period): """ Returns an extended dataframe resampled_data after adding new data. """ # Concat data - concat_data = pd.concat([data, bgnet_df, sens_df], axis=1) + concat_data = pd.concat([data, bgnet_df], axis=1) # Optionally zero masked data if mask is not None: @@ -680,13 +578,10 @@ def validate_save(config, daily_plots=True, decay_plots=True, top_up_plots=True, logging.error("Cannot find predictions {}(pred) != {}(data)".format(bgnet_df.index[0], data.index[0]) ) continue - sens_df = hitrate.get_sensitivity_predictions(bgnet_df, data) - resampled_data = add_to_resampled_data( resampled_data, data, bgnet_df, - sens_df, obs_mask, resample_period=config['validation']['validation_resample_periode'], ) @@ -710,7 +605,6 @@ def validate_save(config, daily_plots=True, decay_plots=True, top_up_plots=True, tag='decay', plot_path=plot_path+'/decays', show_figs=False, - sens_df=sens_df, ) if os.path.isfile(plot_path+'/decays_summary.csv'): @@ -739,7 +633,6 @@ def validate_save(config, daily_plots=True, decay_plots=True, top_up_plots=True, tag='top_up', plot_path=plot_path+'/top_ups', show_figs=False, - sens_df=sens_df, ) if os.path.isfile(plot_path+'/top_ups_summary.csv'): @@ -763,7 +656,6 @@ def validate_save(config, daily_plots=True, decay_plots=True, top_up_plots=True, tag='study_times', plot_path=plot_path+'/study_times', show_figs=False, - sens_df=sens_df, write_html=True, annotate=False, ) @@ -786,8 +678,6 @@ def validate_save(config, daily_plots=True, decay_plots=True, top_up_plots=True, # Plot sensitivities for whole month make_sensitivity_overview_plots(obs, resampled_data, plot_path, show_figs=False, write_html=True, write_image=True) - print_master_formula(config, plot_path) - # Save config as well pickle.dump( config, diff --git a/bg_net/plot_utils.py b/bg_net/plot_utils.py index 5c50a7f479a472767688237556f352da6e9f00d0..50cb104e85aee19a8a6b2854ca8911ae18f38d16 100644 --- a/bg_net/plot_utils.py +++ b/bg_net/plot_utils.py @@ -19,41 +19,34 @@ pd.set_option('mode.chained_assignment', 'raise') def plot_pred_timeseries_only(obs, data, bgnet_df): """ - Returns plotly figure object with pred timeries + Returns plotly figure object with predicted timeries + + obs: str + Name of target hit rate + + data: pandas.DataFrame + DataFrama with target hit rate + + bgnet_df: pandas.DataFrame + DataFrame with bgnet predictions + + returns: plotly.graph_objects.Figure + Handle to time series plot """ - fig = go.Figure() + index = bgnet_df.index - fig.add_trace(go.Scatter(x=data.index, y=data[obs], - name='Observed',) - ) - fig.add_trace(go.Scatter(x=data.index, y=bgnet_df['pred'], - name='BGNet Total',) - ) - fig.add_trace(go.Scatter(x=data.index, y=bgnet_df['pred_inj_her'], - name='BGNet HER Inj.') - ) - fig.add_trace(go.Scatter(x=data.index, y=bgnet_df['pred_inj_ler'], - name='BGNet LER Inj.',) - ) - fig.add_trace(go.Scatter(x=data.index, y=bgnet_df['pred_touschek_ler'], - name='BGNet LER Touschek', ) - ) - fig.add_trace(go.Scatter(x=data.index, y=bgnet_df['pred_beam_gas_ler'], - name='BGNet LER Beam-gas',) - ) - fig.add_trace(go.Scatter(x=data.index, y=bgnet_df['pred_touschek_her'], - name='BGNet HER Touschek') - ) - fig.add_trace(go.Scatter(x=data.index, y=bgnet_df['pred_beam_gas_her'], - name='BGNet HER beam-gas') - ) - fig.add_trace(go.Scatter(x=data.index, y=bgnet_df['pred_lumi'], - name='BGNet Lumi', line_dash='dash',) - ) - fig.add_trace(go.Scatter(x=data.index, y=bgnet_df['pred_noise'], - name='BGNet Noise', line_dash='dash',) - ) + fig = go.Figure() + fig.add_trace(go.Scatter(x=index, y=data[obs], name='Observed')) + fig.add_trace(go.Scatter(x=index, y=bgnet_df['pred'],name='BGNet Total')) + fig.add_trace(go.Scatter(x=index, y=bgnet_df['pred_inj_her'],name='BGNet HER Inj.')) + fig.add_trace(go.Scatter(x=index, y=bgnet_df['pred_inj_ler'],name='BGNet LER Inj.',)) + fig.add_trace(go.Scatter(x=index, y=bgnet_df['pred_touschek_ler'],name='BGNet LER Touschek', )) + fig.add_trace(go.Scatter(x=index, y=bgnet_df['pred_beam_gas_ler'],name='BGNet LER Beam-gas',)) + fig.add_trace(go.Scatter(x=index, y=bgnet_df['pred_touschek_her'],name='BGNet HER Touschek')) + fig.add_trace(go.Scatter(x=index, y=bgnet_df['pred_beam_gas_her'],name='BGNet HER beam-gas')) + fig.add_trace(go.Scatter(x=index, y=bgnet_df['pred_lumi'],name='BGNet Lumi', line_dash='dash',)) + fig.add_trace(go.Scatter(x=index, y=bgnet_df['pred_noise'],name='BGNet Noise', line_dash='dash',)) fig.update_layout(legend=dict(itemsizing='constant')) diff --git a/bg_net/report_utils.py b/bg_net/report_utils.py index caf38ce620aa8b7181ad86ea9cc7237c9ecee2c3..8a2ebca665b6fd1e351f69b0e7068f0c68e89b5a 100644 --- a/bg_net/report_utils.py +++ b/bg_net/report_utils.py @@ -77,7 +77,7 @@ def hitrate_two_axis_plot( PV_pos, reference_data_times, test_data_times, - obs_unit="NA", + obs_unit, PV_unit="NA", look_back=7): """ @@ -92,14 +92,13 @@ def hitrate_two_axis_plot( par1 = host.twinx() host.ticklabel_format(useOffset=False) par1.ticklabel_format(useOffset=False) - host.set_ylabel("{} [{}]".format(obs, get_unit(obs))) + host.set_ylabel("{} [{}]".format(sub_model_pred_name, obs_unit)) host.set_xlabel("Date") par1.set_ylabel("{} [{}]".format(PV, get_unit(PV))) - p0, = host.plot(index, y, color='black', label=obs) - p1, = host.plot(index, sub_model_pred_y, color='blue', label=sub_model_pred_name) - p2, = par1.plot(index, x[:, PV_pos], color='red', label=PV) + p0, = host.plot(index, sub_model_pred_y, color='blue', label=sub_model_pred_name) + p1, = par1.plot(index, x[:, PV_pos], color='red', label=PV) for i in range(len(reference_data_times)): h1 = host.axvspan(pd.to_datetime(reference_data_times[i][0]), pd.to_datetime(reference_data_times[i][1]), @@ -107,11 +106,11 @@ def hitrate_two_axis_plot( for i in range(len(test_data_times)): h2 = host.axvspan(pd.to_datetime(test_data_times[i][0]), pd.to_datetime(test_data_times[i][1]), label="Test", color="crimson", alpha=0.3, ymin=0, ymax=100) - lns = [p0, p1, p2, h1, h2] + lns = [p0, p1, h1, h2] host.legend(handles=lns, loc='upper center', bbox_to_anchor=(0.5, 1.2), fancybox=True, shadow=True, ncol=3) host.yaxis.label.set_color(p0.get_color()) - par1.yaxis.label.set_color(p2.get_color()) + par1.yaxis.label.set_color(p1.get_color()) start_time = index[-1] - timedelta(days=look_back) end_time = index[-1] @@ -260,6 +259,24 @@ def make_default_plots( plot.close('all') +def gate_attention_plots(weights_her, weights_ler): + plt.rcParams.update({'font.size': 20}) + fig, host = plt.subplots(figsize=(16, 10)) # (width, height) in inches + host.ticklabel_format(useOffset=False) + plt.title("Beam Gate Status attention weights") + + host.set_xlabel("Time shift in seconds") + host.set_ylabel("Attention weights") + + shifts = np.arange(0, len(weights_her)) + host.plot(shifts, weights_her, color='blue', label="HER Beam Gate Status") + shifts = np.arange(0, len(weights_ler)) + host.plot(shifts, weights_ler, color='red', label="LER Beam Gate Status") + host.grid() + host.legend(loc='best') + return plt + + def hitrate_prediction_plot(obs, decomp_bgnet, reference_data_times, test_data_times, look_back=7): plt.rcParams['timezone'] = 'Japan' plt.rcParams.update({'font.size': 20}) @@ -271,7 +288,7 @@ def hitrate_prediction_plot(obs, decomp_bgnet, reference_data_times, test_data_t ax[0].set_xlabel("Date") ax[0].set_ylabel("{} [{}]".format(obs, get_unit(obs))) ax[0].xaxis.set_major_formatter(mdates.DateFormatter('%m-%d %H:%M')) - loc = plticker.MultipleLocator(base=0.25) + loc = plticker.MultipleLocator(base=(1/24*ceil(0.05*look_back*24))) ax[0].xaxis.set_major_locator(loc) ax[0].grid() ax[1].yaxis.set_label_position("right") @@ -335,6 +352,7 @@ def hitrate_prediction_plot(obs, decomp_bgnet, reference_data_times, test_data_t ) for stack, hatch in zip(stacks, hatches): stack.set_hatch(hatch) + stack.set_edgecolor([0.0, 0.0, 0.0, 1]) ax[0].plot(decomp_bgnet.index, decomp_bgnet[obs], diff --git a/bg_net/run_training.py b/bg_net/run_training.py index 03f1f32a28177665315f7e78962f96d562c02b04..e20d9e406099a1c2194c66d93cedf3148cfea464 100644 --- a/bg_net/run_training.py +++ b/bg_net/run_training.py @@ -135,7 +135,7 @@ def process_const_model_params(const_params, model, scalers): else: value = const_params[0] if const_params[1]: - value = hitrate.get_scaled_sensitivty(value, model=model, scalers_x=scalers_x, scaler_y=scaler_y) + value = hitrate.get_scaled_sensitivity(value, model=model, scalers_x=scalers_x, scaler_y=scaler_y) return value diff --git a/examples/BGNetTuning.ipynb b/examples/BGNetTuning.ipynb deleted file mode 100644 index 2c2c3c6f2ee1e16b8e0a3cad4098f0c9455bf982..0000000000000000000000000000000000000000 --- a/examples/BGNetTuning.ipynb +++ /dev/null @@ -1,412 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "id": "planned-project", - "metadata": {}, - "outputs": [], - "source": [ - "# BGNet (Background prediction with neural nets for Belle II at SuperKEKB)\n", - "# Author: The BGNet developers \n", - "# \n", - "# See git log for contributors and copyright holders. \n", - "# This file is licensed under MIT licence, see LICENSE.md. \n", - "\n", - "\n", - "# BGNetTuning: Jupyter notebook for Belle II background studies \n", - "#\n", - "# This notebook performs hyperparameter optimization for BGNet using ray tune. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "copyrighted-locator", - "metadata": {}, - "outputs": [], - "source": [ - "import pandas as pd\n", - "import numpy as np\n", - "import pickle\n", - "import pytz\n", - "import datetime\n", - "import shutil\n", - "import glob\n", - "import os\n", - "import tqdm\n", - "import yaml\n", - "import random\n", - "import csv\n", - "import sys\n", - "import matplotlib.pyplot as plt\n", - "import plotly.express as px\n", - "import plotly.graph_objects as go\n", - "from pathlib import Path\n", - "\n", - "from ray import tune\n", - "import ray\n", - "import tensorflow.keras as keras\n", - "from tensorflow.keras.callbacks import ModelCheckpoint\n", - "\n", - "import bg_net.archiver.skb_naming as skb_naming\n", - "import bg_net.archiver.indexer_utils as indexer_utils\n", - "import bg_net.plot_utils as plot_utils\n", - "import bg_net.run_training as run_training\n", - "import bg_net.bgnet as bgnet\n", - "\n", - "# This seeds the hyperparameter sampling.\n", - "np.random.seed(5) " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "classical-delaware", - "metadata": {}, - "outputs": [], - "source": [ - "def set_lag(items, new_lag): \n", - " return [ (pv, lag) if lag == 0 else (pv, new_lag) for pv, lag in items ]\n", - "def delete_rand_items(items,n):\n", - " to_delete = set(random.sample(range(len(items)),n))\n", - " return [x for i,x in enumerate(items) if not i in to_delete]\n", - "\n", - "\n", - "def add_rand_items(items, n, cand_items): \n", - " new_items = items + random.sample(cand_items,n)\n", - " return new_items\n", - "\n", - "def mutate_inj(inj_pvs, inj_cands, del_max=1, add_max=3, lag_max=7): \n", - " \"\"\"Mutate the specifications for preprocessing inputs for inj. nets\"\"\"\n", - " \n", - " # Remove features from default list\n", - " n_del = random.sample(range(del_max+1),k=1)[0]\n", - " n_add = random.sample(range(add_max+1),k=1)[0]\n", - " n_lag = 7 #random.sample(range(lag_max+1),k=1)[0]\n", - " \n", - " # Delete PVs from the current list\n", - " inj_pvs = delete_rand_items(inj_pvs, n=n_del)\n", - " # Add PVs from the candidate list\n", - " inj_pvs = add_rand_items(inj_pvs, n=n_add, cand_items=inj_cands)\n", - " # Adjust the lag\n", - " inj_pvs = set_lag(inj_pvs, n_lag)\n", - " \n", - " return inj_pvs\n", - "\n", - "class TuneReporterCallback(keras.callbacks.Callback):\n", - " \"\"\"Tune Callback for Keras.\n", - " \n", - " The callback is invoked every epoch.\n", - " \"\"\"\n", - "\n", - " def __init__(self, logs={}):\n", - " self.iteration = 0\n", - " super(TuneReporterCallback, self).__init__()\n", - "\n", - " def on_epoch_end(self, batch, logs={}):\n", - " self.iteration += 1\n", - " tune.report(keras_info=logs, mean_loss=logs.get(\"loss\"))\n", - "\n", - "\n", - "def tune_bgnet(config, checkpoint_dir=None): \n", - " # Run this to make pathes local to trial folder\n", - " config = run_training.create_train_config(config)\n", - " \n", - " # Run training with callbacks to tuning\n", - " run_training.train_save(config, callbacks=[TuneReporterCallback()]) \n", - "\n", - "def run_tuning(config, num_samples, experiment_name, local_dir):\n", - "\n", - " ray.shutdown() # Restart Ray defensively in case the ray connection is lost. \n", - " #ray.init(num_cpus=2)\n", - " ray.init()\n", - " \n", - " analysis = tune.run(\n", - " tune_bgnet, \n", - " verbose=1, \n", - " config=config,\n", - " num_samples=num_samples, \n", - " #resources_per_trial={\"cpu\": 2},\n", - " local_dir=local_dir,\n", - " name=experiment_name,\n", - " )\n", - " \n", - " # Will shutdown the ray cluster\n", - " ray.shutdown()\n", - " \n", - " return analysis " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "middle-thinking", - "metadata": {}, - "outputs": [], - "source": [ - "# Set a single observable for optimization\n", - "obs =\"CDC_CUR_AVERAGE\" \n", - "\n", - "# Set data path \n", - "data_path = os.path.expandvars('${B2_ARCHIVER}/skb_archiver_data_2020_06')\n", - "\n", - "if not os.path.isdir(data_path):\n", - " print('Cannot find data dir at ', data_path)\n", - " sys.exit(1)\n", - "\n", - "\n", - "# Put all logging of ray here \n", - "local_dir = \"./ray_tune_results\"\n", - "\n", - "# Experiment name\n", - "experiment_name = 'exp1'\n", - "\n", - "# Full path to logs \n", - "logdir = local_dir+'/'+experiment_name\n", - "\n", - "# We clean out the logs before running for a clean visualization later.\n", - "shutil.rmtree(logdir, ignore_errors=True)\n", - "\n", - "her_cands = [\n", - " (\"CGHOPT_IP_BETA_X\",0),\n", - " (\"CGHOPT_IP_BETA_Y\",0), \n", - " (\"CGLOPT_IP_BETA_X\",0),\n", - " (\"CGLOPT_IP_BETA_Y\",0), \n", - " (\"SKB_HER_P_avg\",0),\n", - " (\"SKB_LER_P_avg\",0),\n", - " ('HER_BUNCH_CURRENT',0),\n", - " (\"BMHXRM_BEAM_EMITTX\",4),\n", - " (\"BMHXRM_BEAM_EMITTY\",4),\n", - " (\"TM_EVR0_HER_INJ_EFF_BCM\",4),\n", - " (\"CGeBT_DPP0\",4),\n", - " (\"CGHINJ_SEPTUM_ANG_R\",4),\n", - " (\"CGHINJ_SEPTUM_POS_R\",4),\n", - " (\"CGHINJ_KICKER_HEIGHT_R\",4),\n", - " (\"CGHINJ_KICKER_JUMP_R\",4),\n", - " (\"CGHOPT_GATED_TUNE_X_MEAS\",4), \n", - " (\"CGHOPT_GATED_TUNE_Y_MEAS\",4), \n", - "]\n", - "\n", - "ler_cands = [\n", - " (\"CGLOPT_IP_BETA_X\",0),\n", - " (\"CGLOPT_IP_BETA_Y\",0), \n", - " (\"CGHOPT_IP_BETA_X\",0),\n", - " (\"CGHOPT_IP_BETA_Y\",0), \n", - " (\"SKB_LER_P_avg\",0),\n", - " (\"SKB_HER_P_avg\",0),\n", - " ('LER_BUNCH_CURRENT',0),\n", - " (\"BMLXRM_BEAM_EMITTX\",4),\n", - " (\"BMLXRM_BEAM_EMITTY\",4),\n", - " (\"TM_EVR0_LER_INJ_EFF_BCM\",4),\n", - " (\"CGpBT_DPP0\",4),\n", - " (\"CGLINJ_SEPTUM_ANG_R\",4),\n", - " (\"CGLINJ_SEPTUM_POS_R\",4),\n", - " (\"CGLINJ_KICKER_HEIGHT_R\",4),\n", - " (\"CGLINJ_KICKER_JUMP_R\",4),\n", - " (\"CGLOPT_GATED_TUNE_X_MEAS\",4),\n", - " (\"CGLOPT_GATED_TUNE_Y_MEAS\",4),\n", - "]\n", - "\n", - "# Get a default config file\n", - "with open('configs/training_default.yaml') as f:\n", - " config = yaml.load(f, Loader=yaml.FullLoader)\n", - " \n", - "# Make format conversions (str-> timestamps) \n", - "config = run_training.convert_config(config) \n", - "\n", - "# For training only storage nets\n", - "config['training']['two_stage_fit'] = False\n", - "\n", - "default_her_inj_pvs = config['preparation']['her_inj_pvs']\n", - "default_ler_inj_pvs = config['preparation']['ler_inj_pvs']\n", - "\n", - "# Tune feature selection\n", - "bgnet_features_space = {\n", - " #\"lag\" : 7, #tune.grid_search([3,4,5,6,7,8]),\n", - " #\"n_add\" : tune.choice([0,1,2]),\n", - " #\"n_del\" : tune.choice([0,1,2]),\n", - "}\n", - "\n", - "# Tune the preprocessing of inputs here\n", - "bgnet_preparation_space = {\n", - " #\"use_study_time_weighting\": tune.grid_search([True, False]), # defaults to True\n", - " #\"train_test_split\" : tune.choice([0.1, 0.3, 0.5]), #defaults to 0.5\n", - " #\"weights_scale\" : 1, #tune.choice([1, 2, 3]), #defaults to 3.0 but 1 seems better\n", - " #'her_inj_pvs' : tune.sample_from(\n", - " # lambda _: mutate_inj(default_her_inj_pvs, her_cands) \n", - " #),\n", - " #'ler_inj_pvs' : tune.sample_from(\n", - " # lambda _: mutate_inj(default_ler_inj_pvs, ler_cands) \n", - " #),\n", - "}\n", - "\n", - "# Tune the network arch or training loop here\n", - "bgnet_training_space = {\n", - " \"lr\": tune.loguniform(0.00001, 0.1), # 0.0004, \n", - "}\n", - "\n", - "num_samples = 4\n", - "\n", - "# Set config\n", - "config = run_training.create_train_config(config, obs=obs, data_path=data_path)\n", - "config['feature_selection'] = bgnet_features_space\n", - "config['training'].update(bgnet_training_space)\n", - "config['preparation'].update(bgnet_preparation_space)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "accompanied-karma", - "metadata": {}, - "outputs": [], - "source": [ - "# This is the meat \n", - "analysis = run_tuning(config, num_samples, experiment_name, local_dir)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "intimate-population", - "metadata": {}, - "outputs": [], - "source": [ - "ray.shutdown()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "opponent-group", - "metadata": {}, - "outputs": [], - "source": [ - "# Picke observable specific config with scalers \n", - "pickle.dump(analysis, open('./tune_analysis.pkl', 'wb'))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "universal-omega", - "metadata": {}, - "outputs": [], - "source": [ - "print(\"Best config: \", analysis.get_best_config(\n", - " metric=\"keras_info/val_loss\", mode=\"min\"))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "french-pierre", - "metadata": {}, - "outputs": [], - "source": [ - "dfs = analysis.trial_dataframes\n", - "\n", - "# Plot by epoch\n", - "ax = None # This plots everything on the same plot\n", - "for d in dfs.values():\n", - " ax = d['keras_info/val_loss'].plot(ax=ax, legend=False)\n", - "\n", - " \n", - "ax.set_xlabel(\"epoch\")\n", - "ax.set_ylabel(\"validation loss\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "mounted-neutral", - "metadata": {}, - "outputs": [], - "source": [ - "# Obtain the directory where the best model is saved.\n", - "logdir = analysis.get_best_logdir(\"keras_info/val_loss\", mode=\"min\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "fresh-victim", - "metadata": {}, - "outputs": [], - "source": [ - "print(logdir)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "danish-broadway", - "metadata": {}, - "outputs": [], - "source": [ - "all_configs = analysis.get_all_configs()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "eastern-subscription", - "metadata": {}, - "outputs": [], - "source": [ - "summary = []\n", - "\n", - "for k, d in dfs.items():\n", - " trial_config = all_configs[k]\n", - " val_loss = d['keras_info/val_loss'].min()\n", - " lr = trial_config['training']['lr']\n", - " her_inj_pvs= trial_config['preparation']['her_inj_pvs']\n", - " ler_inj_pvs = trial_config['preparation']['ler_inj_pvs']\n", - " \n", - " summary.append({'her_inj_pvs': str(her_inj_pvs),'ler_inj_pvs': str(ler_inj_pvs), 'val_loss': val_loss, 'lr':lr})\n", - " \n", - "summary = sorted(summary, key=lambda x: x.get('val_loss')) \n", - " \n", - "with open('dict.txt', 'w') as csv_file: \n", - " writer = csv.writer(csv_file)\n", - " #for key, value in mydict.items():\n", - " for trial in summary: \n", - " writer.writerow([trial['val_loss']])\n", - " writer.writerow([trial['lr']])\n", - " writer.writerow([trial['her_inj_pvs']])\n", - " writer.writerow([trial['ler_inj_pvs']])\n", - " writer.writerow([])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "protecting-implementation", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.8.8" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/examples/HitRate.ipynb b/examples/HitRate.ipynb index 74ff5ce8e71b9f4815e3822462093fe742f85083..20d7dd344b099bcccd4f4af6f836b689a01a8b64 100644 --- a/examples/HitRate.ipynb +++ b/examples/HitRate.ipynb @@ -10,19 +10,31 @@ "# Author: The BGNet developers \n", "# \n", "# See git log for contributors and copyright holders. \n", - "# This file is licensed under MIT licence, see LICENSE.md. \n", + "# This file is licensed under MIT licence, see LICENSE.md. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Neural network model for background prediction at Belle II\n", + "\n", + "This notebooks trains a BGNET model for predicting background hit rate and its decomposition into their\n", + "physical sources. One group of background sources are storage backgrounds for the electron (HER) and positron\n", + "(LER) rings. Storage backgrouds can be further differentiated by the loss mechanism of beam particles: Touschek \n", + "intra bunch scattering or Beam-gas scattering of a beam particles with residual gas in the ring. Another source\n", + "of backgrounds originate from top up injections into the rings needed to maintain beam currents over long \n", + "times of continous data taking. The final sources of backgrounds are luminosity backgrounds, i.e. QED \n", + "process at the interaction point of Belle II with high cross sections like radiative Bhabba scattering or two \n", + "photon processes. \n", "\n", + "For more information, see draft paper on overleaf https://www.overleaf.com/read/qmtdjjfvzjnz\n", "\n", - "# HitRate: Jupyter notebook for Belle II background studies \n", - "#\n", - "# This notebooks trains a BGNET model for predicting hit rates and computes feature attributions and interactions to\n", - "# give insigth into what the model has learned. \n", - "#\n", - "# This notebook is part of bg_net project: https://gitlab.desy.de/benjamin.schwenker/bg_net.git\n", - "# Access to background data can be provided by author on reasonable request. \n", - "#\n", - "# This notebook uses the path_explain project (https://github.com/suinleelab/path_explain) to compute feature \n", - "# attributions and plot them. " + "This notebook is part of bg_net project: https://gitlab.desy.de/benjamin.schwenker/bg_net.git\n", + "Access to background data can be provided by author on reasonable request. \n", + "\n", + "This notebook uses the path_explain project (https://github.com/suinleelab/path_explain) to compute feature \n", + "attributions and plot them. " ] }, { @@ -39,6 +51,7 @@ "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", + "import pandas as pd\n", "import seaborn as sns\n", "import yaml\n", "from path_explain import PathExplainerTF, summary_plot\n", @@ -49,6 +62,7 @@ "import bg_net.hitrate as hitrate\n", "import bg_net.run_training as run_training\n", "import bg_net.make_explanations as make_explanations\n", + "import bg_net.plot_utils\n", "from bg_net.archiver.indexer_utils import str_to_jst_datetime" ] }, @@ -74,18 +88,20 @@ "metadata": {}, "outputs": [], "source": [ - "# Set a single observable for optimization\n", + "# Select a hit rate of a Belle II subdetector as target varible for training\n", "obs = \"VXD_Rad_BP_BW_325_DoseRate\"\n", + "\n", + "# Set unit of target variable \n", "obs_unit = \"mRad/s\"\n", "\n", - "# Set data path\n", - "data_path = os.path.expandvars(\"${B2_ARCHIVER}\")\n", + "# Set path to data, adjust to your system\n", + "data_path = os.path.expandvars(\"/home/benjamin/b2_archiver_v6\")\n", "\n", - "# Time interval to look at\n", + "# Select time intervall used for trainig \n", "start = str_to_jst_datetime(\"{:04d}-{:02d}-{:02d} 00:00:00\".format(2021, 6, 1))\n", "end = str_to_jst_datetime(\"{:04d}-{:02d}-{:02d} 23:59:59\".format(2021, 6, 30))\n", "\n", - "# Get a default config file\n", + "# Get a default training config file\n", "with open('configs/training_default.yaml') as f:\n", " config = yaml.load(f, Loader=yaml.FullLoader)\n", " \n", @@ -260,16 +276,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Interpreting the Model(s)" + "# Interpreting the Model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Path explain API requires that the model to be interpreted has a single input tensor. The plotting needs a list of all input feature names.\n", - "\n", - "We focus on explaining the injection and storage BG sub models of our freshly trained hitrate model. Since these models take a pair of input tensors, we must concatenate those tensors and feed the result into an interpretable model that internally splits the input tensor intp a pair of tensors for the original hitrate sub model. This procedure does not add any new trainable parameters and does not change any prediction result. " + "Path explain API requires that the model to be interpreted has a single input tensor and a scalar output. The plotting needs a list of all input feature names. We explain the injection and storage BG sub models of our freshly trained hitrate model. In order to get the sub models and their inputs, we use some helper functions." ] }, { @@ -278,8 +292,15 @@ "metadata": {}, "outputs": [], "source": [ - "sub_model_names = ['injection_ler', 'injection_her', ]\n", - "sub_model_name = sub_model_names[1]" + "sub_model_names = [\n", + " 'injection_her',\n", + " 'injection_ler', \n", + " 'her_beam_gas', \n", + " 'her_touschek', \n", + " 'ler_beam_gas', \n", + " 'ler_touschek',\n", + "]\n", + "sub_model_name = sub_model_names[0]" ] }, { @@ -290,12 +311,16 @@ "source": [ "npoints = 500 \n", "\n", - "interpret_x_train = hitrate.get_interpretable_inputs(x_train, sub_model_name, size=npoints) \n", - "interpret_x_train_scaled = hitrate.get_interpretable_inputs(x_train_scaled, sub_model_name, size=npoints)\n", - "interpret_x_test = hitrate.get_interpretable_inputs(x_test, sub_model_name, size=npoints)\n", - "interpret_x_test_scaled = hitrate.get_interpretable_inputs(x_test_scaled, sub_model_name, size=npoints)\n", + "interpret_x_train_scaled = hitrate.get_interpretable_inputs(x_train_scaled, sub_model_name)\n", + "interpret_x_test = hitrate.get_interpretable_inputs(x_test, sub_model_name)\n", + "interpret_x_test_scaled = hitrate.get_interpretable_inputs(x_test_scaled, sub_model_name)\n", "\n", - "interpret_model = hitrate.get_intepretable_model(model, num_inputs=interpret_x_train.shape[1], name=sub_model_name)\n", + "interpret_x_train_scaled = bgnet.downsample_data_to_npoints(interpret_x_train_scaled, npoints)\n", + "interpret_x_test = bgnet.downsample_data_to_npoints(interpret_x_test, npoints)\n", + "interpret_x_test_scaled = bgnet.downsample_data_to_npoints(interpret_x_test_scaled, npoints)\n", + "\n", + "\n", + "interpret_model = hitrate.get_intepretable_model(model, num_inputs=interpret_x_train_scaled.shape[1], name=sub_model_name)\n", "\n", "feature_names = hitrate.get_interpretable_feature_names(config, sub_model_name) " ] @@ -361,23 +386,18 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Understanding delays for non synchronous data " + "# Understanding the timing of injections" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "A binary input feature called `BEAM_GATE_STATUS` typically figures among the most important inputs. This variable tells us when the injection gate is open/close and injections into the ring could happen. A closed gate means we can expect no top up injections into the riung and no injection background rate at the Belle II detector. \n", - "\n", - "It is a good sign that the `BEAM_GATE_STATUS` is ranked highly. \n", + "An important input variable is the `BEAM_GATE_STATUS` for the HER (LER) ring. This binary variable reports when the injection gate is open (1) and injections into the ring are performed. A closed (0) gate reports that no injection are performed and no additional hitrate from injections is expected. \n", "\n", - "The injection networks not only get the `BEAM_GATE_STATUS` at time t as input, but also values at previous time t-1s, t-2s, ..., t-Ns. The feature of the `BEAM_GATE_STATUS` with lag i is called `BEAM_GATE_STATUS_shift_i`. If all variables are perfectly syncronized, there is no need for lagged inputs and the network would have all reason to ignore these inputs. However, it is observed that the measured rate of a Belle II detector follows the gate status with a delay. Learning this delas is critical for predicting injection backgrounds. \n", + "The computation of online hit rates (used as training targets) gives rise to a delay or time shift of the hit rate relative to the `BEAM_GATE_STATUS`. Even when the gate is opened and top up injection starts at time t, the injected bunches give rise to an increase in hit rate only a few seconds later. Empirically, the delay depends on the Belle II subdetector but is very stable otherwise. \n", "\n", - "For diamonds we see that `BEAM_GATE_STATUS_shift_1` is the most relavant beam gate status. For other Belle II detector this will be different. The most likely reason for delays are processing delays caused by computing one second aggregated hit rates by the subdetector readout. \n", - "\n", - "\n", - "We can easily demonstrate this concept by plotting the model predictions and data with different delays. If our understanding is correct, we expect zero predicted injection background rate when the correctly delayed gate status says closed (0). If we look at a wrong delay, there can be non zero injection backgrounds. " + "The injection network, a submodel of BGNet, learns these time delays. This is implemented by offering time shifted versions of the `BEAM_GATE_STATUS` as input. The feature of the `BEAM_GATE_STATUS` with lag i is called `BEAM_GATE_STATUS_shift_i`. BGNet uses an attention mechanism to compute weights `w_i` for the different time shifts. We display the learned weights below. " ] }, { @@ -386,7 +406,10 @@ "metadata": {}, "outputs": [], "source": [ - "best_gate_index, best_gate_name, best_delay = make_explanations.get_best_gate_status(attributions, feature_names,)" + "weights_her, weights_ler = hitrate.get_predicted_gate_attention_weights(model)\n", + "\n", + "print('Attention weights for HER: ', weights_her)\n", + "print('Attention weights for LER: ', weights_ler)" ] }, { @@ -395,7 +418,18 @@ "metadata": {}, "outputs": [], "source": [ - "interpret_y = interpret_model.predict(interpret_x_test_scaled)" + "from bg_net.report_utils import gate_attention_plots\n", + "\n", + "gate_attention_plots(weights_her, weights_ler)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Testing the decomposition of backgrounds \n", + "\n", + "Now we can use the BGNet model to select a historic time window from the archive, predict the decomposition of background rates and plot the results" ] }, { @@ -404,27 +438,58 @@ "metadata": {}, "outputs": [], "source": [ - "# Scanning the delays, marking the one with highest attributions\n", - "# Here we use the predicted injection rate on y axis. GATE_STATUS off should give zeros when shift is correct. \n", - "import numpy as np\n", + "test_start = str_to_jst_datetime(\"{:04d}-{:02d}-{:02d} 20:00:00\".format(2021, 6, 16))\n", + "test_end = str_to_jst_datetime(\"{:04d}-{:02d}-{:02d} 20:5:00\".format(2021, 6, 16))\n", "\n", - "for gate_index in range(best_gate_index-best_delay, best_gate_index-best_delay+8):\n", + "x_test, y_test, index_test = hitrate.load_data_from_path(\n", + " config,\n", + " data_path=data_path,\n", + " start_time=test_start,\n", + " end_time=test_end,\n", + ")\n", "\n", - " plt.scatter(\n", - " x=np.arange(interpret_y.shape[0]),\n", - " y=interpret_y.squeeze(), \n", - " c=interpret_x_test[:,gate_index ]\n", - " )\n", + "x_scaled_test = bgnet.scale_x(x_test, scalers_x)\n", "\n", - " if gate_index == best_gate_index:\n", - " print('BEST!!')\n", - " \n", - " plt.title('Color coded by {}'.format( feature_names[gate_index] ) )\n", - " plt.xlabel('Data index' )\n", - " plt.ylabel('Predicted injection rate / a.u.')\n", - " \n", - " plt.show()\n", - " " + "bgnet_df = hitrate.get_predictions(\n", + " model=model, inputs=x_scaled_test, scaler_y=scaler_y, index=index_test,\n", + ")\n", + "\n", + "data = pd.DataFrame({obs: y_test[:, 0]}, index=index_test)\n", + "\n", + "fig = bg_net.plot_utils.plot_pred_timeseries_only(obs, data, bgnet_df)\n", + "\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Sensitivity to backgrounds\n", + "\n", + "The intermediat layers of the network can also predict the sensitivity of the Belle II subdetector to different backhground types. \n", + "\n", + "For the injection backgrounds, the sensitivity is the height of the injection spike seen in 1Hz timeseries.\n", + "For storage backgrounds, the sensitivity is equals the size of the components divided by a heuristic feature ($F_B$, $F_T$ or luminosity $\\mathcal{L}$):\n", + "\n", + "\\begin{equation}\n", + " \\mathcal{O}_{\\rm beam-gas} = B \\times IP_{\\rm eff} = B \\times F_{B},\n", + " \\label{eq1}\n", + "\\end{equation}\n", + "\n", + "\\begin{equation}\n", + " \\mathcal{O}_{\\rm Touschek} = T \\times \\frac{I^{2}}{n_{\\rm b}\\sigma_{\\rm x}\\sigma_{\\rm y}\\sigma_{\\rm z}} = T \\times F_{T},\n", + " \\label{eq2}\n", + "\\end{equation}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sens_df = hitrate.get_sensitivity_predictions(model, x_scaled_test, scalers_x, scaler_y, index_test)" ] }, { @@ -432,7 +497,9 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "sens_df.head()" + ] } ], "metadata": { @@ -452,6 +519,11 @@ "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.12" + }, + "vscode": { + "interpreter": { + "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1" + } } }, "nbformat": 4, diff --git a/examples/configs/templates/daily_report_hitrate_template.tex b/examples/configs/templates/daily_report_hitrate_template.tex index fe6f42313cf5595322b833d85035d25eca32877c..629bc49346618d892822078e6aa2a1c44403db15 100644 --- a/examples/configs/templates/daily_report_hitrate_template.tex +++ b/examples/configs/templates/daily_report_hitrate_template.tex @@ -52,14 +52,14 @@ \begin{frame}{Neural network model for hitrate} \small - $O = T_{H}(x_{1}) \cdot F_{T,H} + B_{H}(x_{2}) \cdot F_{B,H} + T_{L}(x_{3}) \cdot F_{T,L} + B_{L}(x_{4}) \cdot F_{B,L} + C_{L} \cdot L + D + INJ_{H}(x_{5}) + INJ_{L}(x_{6})$ \\ + $O = T_{H}(x_{1}) \cdot F_{T,H} + B_{H}(x_{2}) \cdot F_{B,H} + T_{L}(x_{3}) \cdot F_{T,L} + B_{L}(x_{4}) \cdot F_{B,L} + INJ_{H}(x_{5}) + INJ_{L}(x_{6}) + C_{L} \cdot L + D$\\ \vspace{5pt} \normalsize - $F_{{B,H}}=3 \cdot P_{1,H} \cdot I_{H}^{2} + P_{0,H} \cdot I_{H}$\\ + $F_{{B,H}}= I_{H} \cdot P_{{H,eff}}$\\ \vspace{5pt} $F_{T,H}=\frac{I_{H}^{2}}{\sigma_{X,H} \cdot \sigma_{Y,H} \cdot (\sigma_{Z,0,H} + \sigma_{Z,1,H} \cdot I_{H}/n_{b,H}) \cdot n_{b,H}}$\\ \vspace{5pt} - $F_{{B,L}}=3 \cdot P_{1,L} \cdot I_{L}^{2} + P_{0,L} \cdot I_{L}$\\ + $F_{{B,L}}= I_{L} \cdot P_{{L,eff}}$\\ \vspace{5pt} $F_{T,L}=\frac{I_{L}^{2}}{\sigma_{X,L} \cdot \sigma_{Y,L} \cdot (\sigma_{Z,0,L} + \sigma_{Z,1,L} \cdot I_{L}/n_{b,L}) \cdot n_{b,L}}$\\ \vspace{5pt} @@ -73,9 +73,10 @@ \vspace{5pt} $\sigma_{{Z,1,L}} = 1.7642$ mm/mA\\ \vspace{10pt} - $O$ denotes the hitrate as 1Hz time series. \\ - Note that $T_{H,L}$ are $B_{H,L}$ are neural networks modeling the sensitivity to storage bg.\\ - Note that $INJ_{H,L}$ are neural networks modeling the injection bg.\\ + $O$ is the predicted hitrate for OBS_VAR.\\ + $T_{H,L}$ are $B_{H,L}$ are neural networks modeling the sensitivity to storage bg.\\ + $INJ_{H,L}$ are neural networks modeling the injection bg in HER or LER.\\ + See backup slides for listing of input variables for tensors $x_{i}$.\\ \end{frame} \begin{frame}{OBS_VAR} @@ -86,7 +87,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: HER Inj Summary} +\begin{frame}{OBS_VAR: HER Injection Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/feature_attribution_CHANNEL_VAR_injection_her.pdf} @@ -94,7 +95,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: HER Inj Summary} +\begin{frame}{OBS_VAR: HER Injection Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank0_CHANNEL_VAR_injection_her.pdf} @@ -102,7 +103,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: HER Inj Summary} +\begin{frame}{OBS_VAR: HER Injection Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank1_CHANNEL_VAR_injection_her.pdf} @@ -110,7 +111,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: HER Inj Summary} +\begin{frame}{OBS_VAR: HER Injection Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank2_CHANNEL_VAR_injection_her.pdf} @@ -118,8 +119,7 @@ \end{figure} \end{frame} - -\begin{frame}{OBS_VAR: LER Inj Summary} +\begin{frame}{OBS_VAR: LER Injection Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/feature_attribution_CHANNEL_VAR_injection_ler.pdf} @@ -127,7 +127,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: LER Inj Summary} +\begin{frame}{OBS_VAR: LER Injection Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank0_CHANNEL_VAR_injection_ler.pdf} @@ -135,7 +135,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: LER Inj Summary} +\begin{frame}{OBS_VAR: LER Injection Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank1_CHANNEL_VAR_injection_ler.pdf} @@ -143,7 +143,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: LER Inj Summary} +\begin{frame}{OBS_VAR: LER Injection Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank2_CHANNEL_VAR_injection_ler.pdf} @@ -151,7 +151,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: HER Beam Gas Summary} +\begin{frame}{OBS_VAR: HER Beam Gas Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/feature_attribution_CHANNEL_VAR_her_beam_gas.pdf} @@ -159,7 +159,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: HER Beam Gas Summary} +\begin{frame}{OBS_VAR: HER Beam Gas Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank0_CHANNEL_VAR_her_beam_gas.pdf} @@ -167,7 +167,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: HER Beam Gas Summary} +\begin{frame}{OBS_VAR: HER Beam Gas Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank1_CHANNEL_VAR_her_beam_gas.pdf} @@ -175,7 +175,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: HER Beam Gas Summary} +\begin{frame}{OBS_VAR: HER Beam Gas Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank2_CHANNEL_VAR_her_beam_gas.pdf} @@ -183,7 +183,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: LER Beam Gas Summary} +\begin{frame}{OBS_VAR: LER Beam Gas Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/feature_attribution_CHANNEL_VAR_ler_beam_gas.pdf} @@ -191,7 +191,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: LER Beam Gas Summary} +\begin{frame}{OBS_VAR: LER Beam Gas Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank0_CHANNEL_VAR_ler_beam_gas.pdf} @@ -199,7 +199,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: LER Beam Gas Summary} +\begin{frame}{OBS_VAR: LER Beam Gas Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank1_CHANNEL_VAR_ler_beam_gas.pdf} @@ -207,7 +207,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: LER Beam Gas Summary} +\begin{frame}{OBS_VAR: LER Beam Gas Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank2_CHANNEL_VAR_ler_beam_gas.pdf} @@ -215,7 +215,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: HER Touschek Summary} +\begin{frame}{OBS_VAR: HER Touschek Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/feature_attribution_CHANNEL_VAR_her_touschek.pdf} @@ -223,7 +223,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: HER Touschek Summary} +\begin{frame}{OBS_VAR: HER Touschek Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank0_CHANNEL_VAR_her_touschek.pdf} @@ -231,7 +231,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: HER Touschek Summary} +\begin{frame}{OBS_VAR: HER Touschek Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank1_CHANNEL_VAR_her_touschek.pdf} @@ -239,7 +239,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: HER Touschek Summary} +\begin{frame}{OBS_VAR: HER Touschek Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank2_CHANNEL_VAR_her_touschek.pdf} @@ -247,7 +247,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: LER Touschek Summary} +\begin{frame}{OBS_VAR: LER Touschek Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/feature_attribution_CHANNEL_VAR_ler_touschek.pdf} @@ -255,7 +255,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: LER Touschek Summary} +\begin{frame}{OBS_VAR: LER Touschek Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank0_CHANNEL_VAR_ler_touschek.pdf} @@ -263,7 +263,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: LER Touschek Summary} +\begin{frame}{OBS_VAR: LER Touschek Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank1_CHANNEL_VAR_ler_touschek.pdf} @@ -271,7 +271,7 @@ \end{figure} \end{frame} -\begin{frame}{OBS_VAR: LER Touschek Summary} +\begin{frame}{OBS_VAR: LER Touschek Sensitivity} \begin{figure} \centering \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/rank2_CHANNEL_VAR_ler_touschek.pdf} @@ -279,6 +279,14 @@ \end{figure} \end{frame} +\begin{frame}{OBS_VAR: Beam Gate Status attention weights} + \begin{figure} + \centering + \includegraphics[width=\linewidth]{CHANNEL_VAR_plots/attention_weights.pdf} + \caption{Attention weights for different time shifts of Beam Gate Status} + \end{figure} +\end{frame} + INPUTS_VAR \end{document} diff --git a/examples/configs/training_beam_life.yaml b/examples/configs/training_beam_life.yaml index d49fcc2fb93b7ff08edfc4a5c3525044e8b2b3cc..9432a77c4219de971996e8bfdc4feeede910d2a1 100644 --- a/examples/configs/training_beam_life.yaml +++ b/examples/configs/training_beam_life.yaml @@ -81,10 +81,10 @@ preparation: "BEAST_OPT_HER_hor_betak", "DQM_Beam_IP_Y_Mean", - "IGPF_HERDH_SHIFTGAIN", - "IGPF_HERDV_SHIFTGAIN", - "IGPF_HERUH_SHIFTGAIN", - "IGPF_HERUV_SHIFTGAIN", + #"IGPF_HERDH_SHIFTGAIN", # not available June 2021 + #"IGPF_HERDV_SHIFTGAIN", + #"IGPF_HERUH_SHIFTGAIN", + #"IGPF_HERUV_SHIFTGAIN", "RFLNC_D05_DPS_ROOMPHASE_ZERO", "RFLNC_D05_ROOMPHASE", @@ -170,10 +170,10 @@ preparation: "BEAST_OPT_LER_hor_betak", "DQM_Beam_IP_Y_Mean", - "IGPF_LERDH_SHIFTGAIN", - "IGPF_LERDV_SHIFTGAIN", - "IGPF_LERUH_SHIFTGAIN", - "IGPF_LERUV_SHIFTGAIN", + #"IGPF_LERDH_SHIFTGAIN", # not available June 2021 + #"IGPF_LERDV_SHIFTGAIN", + #"IGPF_LERUH_SHIFTGAIN", + #"IGPF_LERUV_SHIFTGAIN", "RFLNC_D05_DPS_ROOMPHASE_ZERO", "RFLNC_D05_ROOMPHASE", diff --git a/examples/configs/training_default.yaml b/examples/configs/training_default.yaml index b55bfa4607e3f9148c0e4fd81a8c0e25515cb98a..27be554203e1c9b6780aa3eee69837db0f7472d9 100644 --- a/examples/configs/training_default.yaml +++ b/examples/configs/training_default.yaml @@ -318,11 +318,11 @@ training: # Number of hidden units for sensitivity net sensitivity_net_hidden_units : 8 # Number of hidden layers for sensitivity net - sensitivity_net_hidden_layers: 1 + sensitivity_net_hidden_layers: 2 # Number of hidden units for injection net injection_net_hidden_units : 32 # Number of hidden layers for injection net - injection_net_hidden_layers: 2 + injection_net_hidden_layers: 3 # Set training validation split validation_split : 0.2 # Verbosity of model.fit() diff --git a/examples/configs/training_inj_duration.yaml b/examples/configs/training_inj_duration.yaml index ed72929050eee1bab2b54314317e88484f546b70..3ad2f5ff88a9f234b4baaef34dc48598fe7143f5 100644 --- a/examples/configs/training_inj_duration.yaml +++ b/examples/configs/training_inj_duration.yaml @@ -83,10 +83,10 @@ preparation: "BEAST_OPT_HER_hor_betak", "DQM_Beam_IP_Y_Mean", - "IGPF_HERDH_SHIFTGAIN", - "IGPF_HERDV_SHIFTGAIN", - "IGPF_HERUH_SHIFTGAIN", - "IGPF_HERUV_SHIFTGAIN", + #"IGPF_HERDH_SHIFTGAIN", # not available June 2021 + #"IGPF_HERDV_SHIFTGAIN", + #"IGPF_HERUH_SHIFTGAIN", + #"IGPF_HERUV_SHIFTGAIN", "RFLNC_D05_DPS_ROOMPHASE_ZERO", "RFLNC_D05_ROOMPHASE", @@ -173,10 +173,10 @@ preparation: "BEAST_OPT_LER_hor_betak", "DQM_Beam_IP_Y_Mean", - "IGPF_LERDH_SHIFTGAIN", - "IGPF_LERDV_SHIFTGAIN", - "IGPF_LERUH_SHIFTGAIN", - "IGPF_LERUV_SHIFTGAIN", + #"IGPF_LERDH_SHIFTGAIN", # not available June 2021 + #"IGPF_LERDV_SHIFTGAIN", + #"IGPF_LERUH_SHIFTGAIN", + #"IGPF_LERUV_SHIFTGAIN", "RFLNC_D05_DPS_ROOMPHASE_ZERO", "RFLNC_D05_ROOMPHASE", diff --git a/examples/configs/training_lumi.yaml b/examples/configs/training_lumi.yaml index 41ef06db728b5c503b771eefac4955160ad12ef9..dae18cc182b6126ffe989b3d05234a140871dfdf 100644 --- a/examples/configs/training_lumi.yaml +++ b/examples/configs/training_lumi.yaml @@ -67,14 +67,14 @@ preparation: "RFLNC_D08_ROOMPHASE_SET", # Feedback system: - "IGPF_HERDH_SHIFTGAIN", - "IGPF_HERDV_SHIFTGAIN", - "IGPF_HERUH_SHIFTGAIN", - "IGPF_HERUV_SHIFTGAIN", - "IGPF_LERDH_SHIFTGAIN", - "IGPF_LERDV_SHIFTGAIN", - "IGPF_LERUH_SHIFTGAIN", - "IGPF_LERUV_SHIFTGAIN", + #"IGPF_HERDH_SHIFTGAIN", # not available June 2021 + #"IGPF_HERDV_SHIFTGAIN", + #"IGPF_HERUH_SHIFTGAIN", + #"IGPF_HERUV_SHIFTGAIN", + #"IGPF_LERDH_SHIFTGAIN", + #"IGPF_LERDV_SHIFTGAIN", + #"IGPF_LERUH_SHIFTGAIN", + #"IGPF_LERUV_SHIFTGAIN", # IP knob scan PVs "CGLOPT_IPr1", diff --git a/examples/goBGNET.py b/examples/goBGNET.py index d12899d895ed52714a406ce61e1e14246bb9fde1..ca71d60d18ce3c2594cec5b8dd5c34f8e06cb899 100644 --- a/examples/goBGNET.py +++ b/examples/goBGNET.py @@ -42,11 +42,7 @@ from bg_net.archiver.indexer_utils import str_to_jst_datetime def get_archives(): """List of tuples (starttime, endtime) for some calendar months to be processed.""" - archives = [ - get_month(year=2020, month=5), - get_month(year=2020, month=6), - get_month(year=2020, month=12), get_month(year=2021, month=6), ] diff --git a/examples/goReporter.py b/examples/goReporter.py index fa4269241d9454c81bb8afc9d6d72080494bf70c..08100dde17f841b5c5d65bec59c3c4c5fe1602a2 100644 --- a/examples/goReporter.py +++ b/examples/goReporter.py @@ -44,7 +44,7 @@ for feature attribution for special reports. """ -__version__ = "0.4.1" +__version__ = "0.5.0" from bg_net.archiver.indexer_utils import str_to_jst_datetime import bg_net.bgnet as bgnet @@ -57,7 +57,7 @@ import bg_net.report_utils import bg_net.report_to_elog from bg_net.report_to_webserver import ReportsToWebServer import bg_net.archiver.indexer_utils as indexer_utils -from bg_net.archiver.skb_naming import get_unit +from bg_net.archiver.skb_naming import get_unit, get_sensitivity_units import bg_net.archiver.connect_vpn as connect_vpn import math @@ -77,7 +77,7 @@ from argparse import ArgumentParser from datetime import datetime, timedelta from calendar import monthrange -from bg_net.report_utils import hitrate_prediction_plot, hitrate_two_axis_plot, complete_index +from bg_net.report_utils import hitrate_prediction_plot, hitrate_two_axis_plot, complete_index, gate_attention_plots import matplotlib.pyplot as plt from path_explain import summary_plot import numpy as np @@ -625,20 +625,24 @@ def report_on_hitrate( x_scaled = bgnet.scale_x(x, scaler_x) - # Create plots - shutil.rmtree(save_path, ignore_errors=True) - Path(save_path).mkdir(parents=True, exist_ok=True) - decomp_bgnet = hitrate.get_predictions( - model=model, x_scaled=x_scaled, scalers_x=scaler_x, scaler_y=scaler_y + model=model, inputs=x_scaled, scaler_y=scaler_y, index=index, ) - for comp in decomp_bgnet: - decomp_bgnet[comp] = decomp_bgnet[comp].squeeze(-1) - decomp_bgnet[obs] = y[:, 0] - decomp_bgnet = pd.DataFrame(data=decomp_bgnet, index=index) + + sens_df = hitrate.get_sensitivity_predictions( + model=model, inputs=x_scaled, scalers_x=scaler_x, scaler_y=scaler_y, index=index, + ) + + # Add column for measured rates + ydf = pd.DataFrame({obs: y[:, 0]}, index=index) + decomp_bgnet = pd.concat([decomp_bgnet, ydf], axis=1) + decomp_bgnet = decomp_bgnet[~decomp_bgnet.index.duplicated(keep='first')] decomp_bgnet = decomp_bgnet.reindex(pd.date_range(start=decomp_bgnet.index[0], end=decomp_bgnet.index[-1], freq="1s")) + sens_df = sens_df[~sens_df.index.duplicated(keep='first')] + sens_df = sens_df.reindex(pd.date_range(start=sens_df.index[0], end=sens_df.index[-1], freq="1s")) + if dump_data: # Dump 1Hz timeseries here report_day = "{:04d}-{:02d}-{:02d}".format(year, month, day) @@ -648,6 +652,14 @@ def report_on_hitrate( # Resample time series decomp_bgnet = decomp_bgnet.resample(resample_frequency).mean() + sens_df = sens_df.resample(resample_frequency).mean() + + ######################### + # Create plots for report + ######################### + + shutil.rmtree(save_path, ignore_errors=True) + Path(save_path).mkdir(parents=True, exist_ok=True) hitrate_prediction_plot( obs, @@ -657,6 +669,7 @@ def report_on_hitrate( look_back=look_back) plt.savefig(os.path.join(save_path, "prediction_{}.pdf".format(channel)), bbox_inches='tight') plt.close('all') + for sub_model_name in hitrate.sub_models: # Create list of names for all input features @@ -698,7 +711,7 @@ def report_on_hitrate( if attributions.shape[0] > 0: # To have unit of attributions same as unit of y - attributions = scaler_y.inverse_transform(attributions) + attributions = scaler_y.inverse_transform(attributions)*hitrate.get_X_scale_factor(scaler_x,sub_model_name) if len(feature_names) < top_k: plot_top_k = len(feature_names) @@ -721,8 +734,7 @@ def report_on_hitrate( bbox_inches='tight') plt.close('all') - sub_model_pred_name = hitrate.map_sub_models_to_decomp[sub_model_name] - sub_model_pred_y = decomp_bgnet[sub_model_pred_name] + sub_model_sens_name = hitrate.map_sub_models_to_sensitivity[sub_model_name] mean_abs = np.mean(np.abs(attributions), axis=0) ranking = [x for (y, x) in sorted(zip(mean_abs.tolist(), feature_names), key=lambda pair: pair[0], reverse=True)] plotting_y, _ = complete_index(y, index, freq="1s", resample_frequency=resample_frequency) @@ -732,26 +744,27 @@ def report_on_hitrate( for i in range(plot_top_k): PV_name = ranking[i] PV_pos = feature_names.index(PV_name) + sens_unit = get_sensitivity_units(obs)[sub_model_sens_name] - if "GATE" not in PV_name and count_two_axis_plots <= 3: + if count_two_axis_plots <= 3: plot = hitrate_two_axis_plot( obs, - sub_model_pred_name=sub_model_pred_name, + sub_model_pred_name=sub_model_sens_name, index=plotting_index, x=interpret_x, y=plotting_y, - sub_model_pred_y=sub_model_pred_y, + sub_model_pred_y=sens_df[sub_model_sens_name], PV=PV_name, PV_pos=PV_pos, reference_data_times=ref_times, test_data_times=test_times, - obs_unit=obs_unit, + obs_unit=sens_unit, PV_unit="NA", look_back=look_back,) plot.title( - "Attribution rank: {} Mean attribution value: {:.3f} {}".format( - i, mean_abs[PV_pos], obs_unit), y=1.2) + "Attribution rank: {} Mean attribution value: {:.4g} {}".format( + i, mean_abs[PV_pos], sens_unit), y=1.2) plot.tight_layout() plot.savefig( os.path.join( @@ -768,8 +781,14 @@ def report_on_hitrate( logging.info('Empty attributions ', attributions.shape) logging.info('Skip attribution plots!') - report_day = get_report_day(year, month, day) + # Plot attention weights for Beam Gate Status PVs + w_her, w_ler = hitrate.get_predicted_gate_attention_weights(model) + gate_attention_plots(w_her, w_ler) + plt.savefig(os.path.join(save_path, "attention_weights.pdf"), bbox_inches='tight') + plt.close('all') + # Compile the report + report_day = get_report_day(year, month, day) compile_latex_report( channel=channel, obs=obs, @@ -786,7 +805,7 @@ def report_on_hitrate( config, name=sub_model) for sub_model in hitrate.sub_models}, report_path=report_path, - ) + ) return