text: import os
  import math 
  from mp.paths import JIP_dir
  import numpy as np
  import matplotlib.pyplot as plt
  from sklearn.preprocessing import StandardScaler
  from sklearn.svm import SVR
  from sklearn.linear_model import Ridge
  from sklearn.neural_network import MLPRegressor
  from sklearn.svm import SVC 
  from sklearn.linear_model import LogisticRegression
  from sklearn.metrics import confusion_matrix, plot_confusion_matrix
  
  #set environmental variables
  #for data_dirs folder, nothing changed compared to Simons version 
  os.environ["WORKFLOW_DIR"] = os.path.join(JIP_dir, 'data_dirs')
  os.environ["OPERATOR_IN_DIR"] = "input"
  os.environ["OPERATOR_OUT_DIR"] = "output"
  os.environ["OPERATOR_TEMP_DIR"] = "temp"
  os.environ["OPERATOR_PERSISTENT_DIR"] = os.path.join(JIP_dir, 'data_dirs', 'persistent')
  
  # preprocessing dir and subfolders 
  os.environ["PREPROCESSED_WORKFLOW_DIR"] = os.path.join(JIP_dir, 'preprocessed_dirs')
  os.environ["PREPROCESSED_OPERATOR_OUT_SCALED_DIR"] = "output_scaled"
  os.environ["PREPROCESSED_OPERATOR_OUT_SCALED_DIR_TRAIN"] = "output_scaled_train"
  
  #dir where train data for intensites is stored (this only needs to be trains_dirs, but since i have more 
  # datasets, another subfolder is here)
  os.environ["TRAIN_WORKFLOW_DIR"] = os.path.join(JIP_dir, 'train_dirs')
  
  #ignore
  ##below is for christian only, used for older data structures where models are trained on
  os.environ["TRAIN_WORKFLOW_DIR_GT"] = os.path.join('Covid-RACOON','All images and labels')
  os.environ["TRAIN_WORKFLOW_DIR_PRED"] = os.path.join('Covid-RACOON','All predictions')
  #os.environ["TRAIN_WORKFLOW_DIR_GT"] = os.path.join('gt_small')
  #os.environ["TRAIN_WORKFLOW_DIR_PRED"] = os.path.join('pred_small')
  
  #which mode is active either 'train' or 'inference' 
  os.environ["INFERENCE_OR_TRAIN"] = 'train'
  
  #ignore
  # the ending of the image files in train_dir is only for older datasets
  os.environ["INPUT_FILE_ENDING"] = 'nii.gz'
  
  from mp.utils.preprocess_utility_functions import basic_preprocessing
  from mp.utils.preprocess_utility_functions import extract_features_all_data,compute_all_prediction_dice_scores
  from mp.utils.feature_extractor import Feature_extractor
  
  # all steps in 1. are not used anymore
  ## 2. function to split data into splits depending on task
  # 2.1 (SET splits for training and eval) set the filter function, such that the splits are as wanted 
  def filter_feature_extr(id,model):
      '''function to decide which split a prediction should be put in.
      Is infered from the id(giving the dataset of the image) and model
      (the dataset the predicting model was trained on)
      Args:
          id(str):giving the dataset of the image
          model(str):the dataset the predicting model was trained on
  
      Returns(str):the split for the pair of id and model, can be other, which can be 
          category of unused data'''
      # # train data 
      # if model[:7] in ['Task740'] and id[:7] in ['Task740','Task541']:
      #     return 'train'
      # # id data
      # if model[:7] in ['Task740'] and id[:7] in ['Task741']:
      #     return 'id_challenge'
      # if model[:7] in ['Task740'] and id[:7] in ['Task542']:
      #     return 'id_inhouse'
      # # od data 
      # if model[:7] in ['Task740'] and id[:7] in ['Task200','Task201']:
      #     return 'ood_mosmed'
  
      # # train data 
      # if model[:7] in ['Task740'] and id[:7] in ['Task740','Task541']:
      #     return 'id_train'
      # # id data
      # if model[:7] in ['Task740'] and id[:7] in ['Task741','Task542']:
      #     return 'id_test'
      # # od data 
      # if model[:7] in ['Task740'] and id[:7] in ['Task200','Task201','Task100','Task101']:
      #     return 'ood_test'
  
      # train data 
      if model[:7] in ['Task740'] and id[:7] in ['Task740']:
          return 'train'
      if model[:7] in ['Task740'] and id[:7] in ['Task541']:
          return 'train'
      # id data
      if model[:7] in ['Task740'] and id[:7] in ['Task741']:
          return 'id_challenge'
      if model[:7] in ['Task740'] and id[:7] in ['Task542']:
          return 'id_inhouse'
      # od data 
      if model[:7] in ['Task740'] and id[:7] in ['Task200','Task201']:
          return 'ood_mosmed'
      if model[:7] in ['Task740'] and id[:7] in ['Task100','Task101']:
          return 'other'
  
      return 'other'
  
  # 2.2 extract the feature
  CLASS = False
  def extract_features_train_id_od(filter,splits,used_feat=[0,1,2,3,4,5]):
      '''extracts the features from the preprocess dictionary, and returns a  list,
      in which each entry corresponds to the features/labels of one split
      Args:
          filter(function: (str,str)->str): should be modified filter_feature_extr
          splits(list(str)):the name of the splits
          used_feat(list(int)):the used features, see in main() for more depth
  
      Returns(list): list, in which each entry corresponds to the features/labels of one split'''
      X = []
      y = []
      paths_to_pred = []
      for i in range(len(splits)):
          X.append([])
          y.append([])
          paths_to_pred.append([])
  
      feat_extr = Feature_extractor()
      work_path = os.path.join(os.environ["PREPROCESSED_WORKFLOW_DIR"],os.environ["PREPROCESSED_OPERATOR_OUT_SCALED_DIR_TRAIN"])
      for id in os.listdir(work_path):
          all_pred_path = os.path.join(work_path,id,'pred')
          if os.path.exists(all_pred_path):
              for model in os.listdir(all_pred_path):
                  split = filter(id,model)
                  for i,s in enumerate(splits):
                      if split == s:
                          feature_path = os.path.join(all_pred_path,model,'features.json')
                          label_path = os.path.join(all_pred_path,model,'dice_score.json')
                          feat_vec = feat_extr.read_feature_vector(feature_path)
                          feat_vec = [feat_vec[index] for index in used_feat]
                          label = feat_extr.read_prediction_label(label_path)
                          if CLASS:
                              label = get_class(label)
                          if np.isnan(np.sum(np.array(feat_vec))) or feat_vec[0]>100000:
                              pass 
                          else:
                              X[i].append(feat_vec)
                              y[i].append(label)
                              paths_to_pred[i].append(os.path.join(all_pred_path,model))
                              
      return X,y,paths_to_pred
  
  def get_class(label):
      if label < 0.2:
          return 1
      elif label < 0.4:
          return 2 
      elif label < 0.6:
          return 3
      elif label < 0.8:
          return 4
      else:
          return 5
  
  def get_float(classes):
      return np.array([(clas-1)/4 for clas in classes])
      
  # 3. loss functions 
  def l2_loss(pred,truth):
      n = len(pred)
      return (1/n)*(np.sum((pred - truth)**2))
  
  def l1_loss(pred,truth, std = False):
      losses = np.absolute(pred-truth)
      mean = np.mean(losses)
      std = np.std(losses)
      if std:
          return mean,std
      return mean
  
  def l1_loss_overestimation(pred,truth,std=False):
      '''computes the error, but underpredictions 
      are counted as correct predictions '''
      losses = pred-truth
      losses_pos = np.array([ max(loss,0) for loss in losses])
      mean = np.mean(losses_pos)
      std = np.std(losses_pos)
      if std:
          return mean,std
      return mean
  
  def get_l1_losses(truth,pred):    
      return np.absolute(truth-pred)    
         
  def confidence_bootstrap(X_train,y_train,X,y,svm,lr,scaler,splits):
      '''This is based on 
      https://machinelearningmastery.com/confidence-intervals-for-machine-learning/
      I simply resample from the training data and then evaluate the performance of the model
      on each boostrap'''
      from numpy.random import randint
      from numpy import percentile
      scores = [[] for _ in range(len(splits))] #scores for smv per run
      scores_lr = [[] for _ in range(len(splits))]# scores for lr per run
      for _ in range(10000): # i conudct 10000 runs 
          indices = randint(0,len(X_train),len(X_train)) #choose the indices in data to use 
          
          X_train_b = X_train[indices]
          y_train_b = y_train[indices]
          svm = SVC(class_weight='balanced')
          lr = LogisticRegression(class_weight='balanced')
          svm.fit(X_train_b,y_train_b) 
          lr.fit(X_train_b,y_train_b) 
  
          for i,_ in enumerate(splits):
          
              X_eval = scaler.transform(X[i])
              y_eval = y[i]
  
              y_svm = svm.predict(X_eval)
              y_lr = lr.predict(X_eval)
  
              conf_ma = conf_matrix(y_eval,y_svm)
              conf_ma_lr = conf_matrix(y_eval,y_lr)
  
              scores[i].append(specificity(conf_ma))
              scores_lr[i].append(specificity(conf_ma_lr))
      lower_p = 2.5
      upper_p = 97.5
      # Now searach for the 95 percent intervalls done like in the link in
      # the description 
      for i in range(len(splits)):
          print(splits[i])
          lower = max(0.0, percentile(scores[i], lower_p))
          print('%.1fth percentile = %.3f' % (lower_p, lower))
          upper = min(1.0, percentile(scores[i], upper_p))
          print('%.1fth percentile = %.3f' % (upper_p, upper))
      print()
      print('LR')
      for i in range(len(splits)):
          print(splits[i])
          lower = max(0.0, percentile(scores_lr[i], lower_p))
          print('%.1fth percentile = %.3f' % (lower_p, lower))
          upper = min(1.0, percentile(scores_lr[i], upper_p))
          print('%.1fth percentile = %.3f' % (upper_p, upper))
  
  def conf_matrix(truth,pred): 
      
      conf_ma =  confusion_matrix(np.array(truth),np.array(pred),labels=[1,2,3,4,5])
      #print(conf_ma)
      ma = np.zeros((2,2))
      ma [0,0] = np.sum(conf_ma[:3,:3])
      ma [0,1] = np.sum(conf_ma[:3,3:])
      ma [1,0] = np.sum(conf_ma[3:,:3])
      ma [1,1] = np.sum(conf_ma[3:,3:])
      return ma
  
  def sensitivity(ma):
      found = int(ma[0,0])
      not_f = int(ma[0,1])
      if found == 0:
          return 0
      else:
          return found/(found+not_f)
  
  def p_value_sensitivity(ma):
      found = int(ma[0,0])
      not_f = int(ma[0,1])
      return find_p_p_test(found,not_f)
  
  def specificity(ma):
      tn = int(ma[1,1])
      fp = int(ma[1,0])
      return tn/(tn+fp)
  
  def p_value_specificity(ma):
      tn = int(ma[1,1])
      fp = int(ma[1,0])
      return find_p_p_test(tn,fp)
  
  def p_test(p,y,n):
      '''takes a p value for the null hypothesis, the number of 
      successfull results and the number of unseccessfull results'''
      prob = 0
      combined = y + n 
      for times in np.arange(y,combined+1,1):
          term = p**times * (1-p)**(combined-times)
          bin = math.comb(combined,times)
          prob += term * bin
      return prob
  
  def find_p_p_test(y,n):
      for p in np.arange(1,0,-0.01):
          # print(p,p_test(p,y,n))
          if p_test(p,y,n) <= 0.05:
              return p, p_test(p,y,n)
  
  def micro_sens(conf_ma):
      conf_ma = np.array(conf_ma)
      tp = 0
      p = np.sum(conf_ma)
      for i in range(0,5):
          tp += conf_ma[i,i]
      return tp/p
  
  def micro_spec(conf_ma):
      conf_ma = np.array(conf_ma)
      tn_cum = 0
      fp_cum = 0
      all = np.sum(conf_ma)
      for i in range(5):
          fp = np.sum(conf_ma[:,i])-conf_ma[i,i]
          tp = conf_ma[i,i]
          fn = np.sum(conf_ma[i,:])-conf_ma[i,i]
          tn = all - (fp+tp+fn)
          tn_cum += tn
          fp_cum += fp
      return tn_cum/(tn_cum+fp_cum)
       
  
  
  # 4. plots for further analysis.
  # (OPTIONAL SET the saving path for the bar plots)
  def l1_loss_bins(pred,truth,split,by_sign=True):
      '''computes and plots a bar plot for the errors of a model, sorted by bins depending on the truth, 
      e.g. the error is taken for all predictions, where the ground trouth is in [0.5,0.6]. The errors 
      are divided into over and under estimation. Additionally the number of data points and the std 
      per bin is displayed
      Args:
          pred: the vector of predictions 
          truth: the vector of ground truths
          split(str): the name of the split
          by_sign(bool): whether the errors should be divided into over and underestimation, setting to 
              false is missing some normally displayed features'''
      save_path = os.path.join('storage','Results','Setup 5','all','Bins_{}'.format(split))
      if not by_sign:
          loss_by_bin = [[] for _ in range(10)]
          bins = np.array([])
          bins_std = np.array([])
          for i in range(10):
              for pre,tru in zip(pred,truth):
                  if i*0.1 <= tru and (i+1)*0.1 > tru:
                      loss_by_bin[i].append(abs(pre-tru))
              if loss_by_bin[i]:
                  bins = np.append(bins,np.mean(loss_by_bin[i]))
                  bins_std = np.append(bins_std,np.std(loss_by_bin[i]))
              else:
                  bins = np.append(bins,0)
                  bins_std = np.append(bins_std,0)
          plt.bar(0.05+np.arange(10)*0.1,height = bins, width = 0.05,yerr = bins_std)
          plt.plot(0.05+np.arange(10)*0.1,[0.1 for i in range(10)],'-r')
          plt.title(split)
          plt.xlabel('bin (e.g. 0.1 to 0.2)')
          plt.ylabel('l1_error')
      else : 
          widths = [[] for _ in range(10)]
          loss_by_bin = [[[],[]] for _ in range(10)]
          bins = [[] for _ in range(10)]
          bins_std = [[] for _ in range(10)]
          for i in range(10):
              for pre,tru in zip(pred,truth):
                  if i*0.1 <= tru and (i+1)*0.1 > tru:
                      err = pre-tru
                      if err >= 0:
                          loss_by_bin[i][1].append(err)
                      else:
                          loss_by_bin[i][0].append(err)
              for j in range(2):
                  if loss_by_bin[i][j]:
                      bins[i].append(np.mean(loss_by_bin[i][j]))
                      bins_std[i].append(np.std(loss_by_bin[i][j]))
                  else:
                      bins[i].append(0)
                      bins_std[i].append(0)
              unders = len(loss_by_bin[i][0])
              overs = len(loss_by_bin[i][1])
              total = unders+overs
              if total > 0:
                  if unders/total > 0.66:
                      widths[i]=[0.05,0.03]
                  elif overs/total >0.66:
                      widths[i]=[0.03,0.05]
                  else: 
                      widths[i]=[0.05,0.05]
              else:
                  widths[i]=[0.05,0.05]
          bins = np.array(bins)
          bins_std = np.array(bins_std)
          widths = np.array(widths)
          plt.bar(0.05+np.arange(10)*0.1,height = bins[:,0], width = widths[:,0],yerr = bins_std[:,0])
          plt.bar(0.05+np.arange(10)*0.1,height = bins[:,1], width = widths[:,1],yerr = bins_std[:,1])
          for i,x in enumerate(0.05+np.arange(10)*0.1):
              len_bin = len(loss_by_bin[i][0])
              if len_bin > 0:
                  plt.text(x,bins[i,0],str(len_bin))
          for i,x in enumerate(0.05+np.arange(10)*0.1):
              len_bin = len(loss_by_bin[i][1])
              if len_bin > 0:
                  plt.text(x,bins[i,1],str(len(loss_by_bin[i][1])))
          plt.plot(0.05+np.arange(10)*0.1,[0.1 for i in range(10)],'-r')
          plt.plot(0.05+np.arange(10)*0.1,[-0.1 for i in range(10)],'-r')
          plt.xlabel('y_true bin ')
          plt.ylabel('approximation error')
      plt.title(split)
      plt.savefig(save_path)
      plt.show()
  
  # (OPTIONAL SET the saving path for the plots)
  def plot_variable_influence(X_train,X,y,splits):
      '''for every feature in the data gives a scatter plot of all predicted segmentations,
      where the y axis is the dice score of the prediction and the x axis is the feature. The 
      different splits are also visualised. 
      Additionally trains a ridge regression, support vector regression and an MLP predictor 
      on the data, who are also displayed
      Args:
          X_train: the data to train the models on 
          X,y: The output of extract_features_train_id_od (so all features/labels sorted into
              groups by splits)
          splits(list): the name of the splits in a list
      '''
      #per variable 
      for i in range(len(X[0][0])):
          
          save_path = os.path.join('storage','Results','Setup 5','Variable {}'.format(i))
          plt.title('Feature {} vs dice'.format(i))
  
          # Plot the models predictions based on single features 
          X_train = [[X[0][k][i]] for k in range(len(X_train))]
          y_train = y[0]
  
          ridge = Ridge(normalize=False)
          svr = SVR()
          mlp = MLPRegressor((50,100,100,50))
  
          scaler = StandardScaler()
          X_train_s = scaler.fit_transform(X_train)
  
          ridge.fit(X_train_s,y_train)
          svr.fit(X_train_s,y_train)
          mlp.fit(X_train_s,y_train)
  
          x_min,x_max = np.min(X_train), np.max(X_train)
          x_span = x_max - x_min
          # take input range to plot model and scale it 
          points = np.reshape(np.arange(x_min,x_max,step=x_span/100),(-1,1))
          points_s = scaler.transform(points)
  
          y_ridge = ridge.predict(points_s)
          y_svr = svr.predict(points_s)
          y_mlp = mlp.predict(points_s)
                 
          plt.plot(points,y_ridge, label='ridge')
          plt.plot(points,y_svr, label='svr')
          plt.plot(points,y_mlp, label='mlp')
  
          #plot the data points
          for j,split in enumerate(splits):
              filt_X = [[X[j][k][i]] for k in range(len(X[j]))]
              plt.scatter(filt_X,y[j],label=split)
  
          plt.legend(loc='upper left')
          plt.xlabel('feature')
          plt.ylabel('dice score of prediction')
          #plt.savefig(save_path)
          plt.show()
  
  def paper_figures_by_split(X_train,X,y,splits):
      #per variable 
      #plt.rcParams.update({'font.size': 20})
      
      for i in range(len(X[0][0])):
          
          title,leg_loc,xlabel,legend_prop = figure_namesand_settings(i)
          save_path = os.path.join('storage','Results','for_paper','Variable {} _lesslabels'.format(i))
          # plt.title(title)
  
          for j,split in enumerate(splits):
              filt_X = [[X[j][k][i]] for k in range(len(X[j]))]
              plt.scatter(filt_X,y[j],label=split)
          if i == 0:
              pass #plt.legend(loc=leg_loc,prop = legend_prop)
          #plt.xlabel(xlabel)
          #plt.ylabel('Dice Score of segmentation')
          plt.xticks(fontsize=14)
          plt.yticks(fontsize=14)
          #plt.savefig(save_path)
          plt.show()
  
  def figure_namesand_settings(feat_nr):
      if feat_nr == 0:
          title = 'Segmentation smoothness vs Dice Score prediction'
          leg_loc = 'upper left'
          xlabel = 'Segmentation smoothness'
          legend_prop = {'size':15}
      if feat_nr == 1:
          title = 'Number Connected Components vs Dice Score prediction'
          leg_loc = 'upper right'
          xlabel = 'Number Connected Components'
          legend_prop = {'size':8}
      if feat_nr == 2:
          title = 'Intensity Mode vs Dice Score prediction'
          leg_loc = 'upper right'
          xlabel = 'Intensity Mode'
          legend_prop = {'size':8}
      return title,leg_loc,xlabel,legend_prop
          
  def find_ex_pred(X,y,paths):
      paths_ex = []
      seen_paths = ['storage\\JIP\\preprocessed_dirs\\output_scaled_train\\Task100_RadiopediaTrain_0007\\pred\\Task740_ChallengeTrainF4',
                      'storage\\JIP\\preprocessed_dirs\\output_scaled_train\\Task200_MosmedTrain_0027\\pred\\Task740_ChallengeTrainF4',
                      'storage\\JIP\\preprocessed_dirs\\output_scaled_train\\Task200_MosmedTrain_0013\\pred\\Task740_ChallengeTrainF4',
                      'storage\\JIP\\preprocessed_dirs\\output_scaled_train\\Task100_RadiopediaTrain_0008\\pred\\Task541_FrankfurtTrainF4']
      for i in range(len(X[0][0])):
          thresh = 0.4
          if i == 0:
              smallest_val = 1
              smallest_ind = [0,0]
              for j in range(len(X)):
                  filt_X = np.array([[X[j][k][i]] for k in range(len(X[j]))])
                  for k in range(len(filt_X)):
                      if filt_X[k] < smallest_val and y[j][k] <= thresh and paths[j][k] not in seen_paths: 
                          smallest_val = filt_X[k]
                          smallest_ind = [j,k]
              paths_ex.append(paths[smallest_ind[0]][smallest_ind[1]])
              print(smallest_val)
              print(y[smallest_ind[0]][smallest_ind[1]])
          if i == 1:
              smallest_val = 1
              smallest_ind = [0,0]
              for j in range(len(X)):
                  filt_X = np.array([[X[j][k][i]] for k in range(len(X[j]))])
                  for k in range(len(filt_X)):
                      if filt_X[k] > smallest_val and y[j][k] <= thresh and paths[j][k] not in paths_ex and paths[j][k] not in seen_paths:
                          smallest_val = filt_X[k]
                          smallest_ind = [j,k]
              paths_ex.append(paths[smallest_ind[0]][smallest_ind[1]])
              print(smallest_val)
              print(y[smallest_ind[0]][smallest_ind[1]])
          if i == 2:
              smallest_val = 0
              smallest_ind = [0,0]
              for j in range(len(X)):
                  filt_X = np.array([[X[j][k][i]] for k in range(len(X[j]))])
                  for k in range(len(filt_X)):
                      if filt_X[k] > smallest_val and y[j][k] <= thresh and paths[j][k] not in paths_ex and paths[j][k] not in seen_paths:
                          smallest_val = filt_X[k]
                          smallest_ind = [j,k]
              paths_ex.append(paths[smallest_ind[0]][smallest_ind[1]])
              print(smallest_val)
              print(y[smallest_ind[0]][smallest_ind[1]])
          if i == 3: 
              smallest_val = 1
              smallest_ind = [0,0]
              for j in range(len(X)):
                  filt_X = np.array([[X[j][k][i]] for k in range(len(X[j]))])
                  for k in range(len(filt_X)):
                      if filt_X[k] < smallest_val and y[j][k] <= thresh and paths[j][k] not in seen_paths: 
                          smallest_val = filt_X[k]
                          smallest_ind = [j,k]
              paths_ex.append(paths[smallest_ind[0]][smallest_ind[1]])
              print(smallest_val)
              print(y[smallest_ind[0]][smallest_ind[1]])
      print(paths_ex)
      return paths_ex
  
  def main(used_feat=[0,1,2,3,4,5],preprocessing=True,feature_extraction=True,
              extract_dice_scores=True,model_train=True,label=1):
      '''
      Args:
          used_feat(list(ints)): a list of ints to tell which features to use for model training:
                                  1: density_distance
                                  2: slice dice (smoothnes of segmentation)
                                  3: slice dice diff (basically the first derivation of 2 (not really used)
                                  4: connected components
                                  5: mean of a gaussian fitted to intensity values
                                  6: variance of a gaussain fitted to intenstiy values
          preprocessing(bool): whether to copy data from a directory into the working directoy
              includes things such as resizing and scaling of images
          train_density(bool): whether to retrain the density for intensity values (time consuming)
          feature_extraction(bool): weather to extract features from the data, not needed, when features.json 
              are already in the working directory(in e.g. subfolder seg for an id)
          extract_dice_scores(bool): whether to extract the dice scores of the predictions, not needed
              when there are already dice_score.json in the working directory(in e.g. subfolder seg for an id)
          model_train(bool): whether to train and evaluate the models '''
      if preprocessing:
          basic_preprocessing(label)
      if feature_extraction:
          extract_features_all_data()
      if extract_dice_scores:
          compute_all_prediction_dice_scores()
      if model_train:
  
          scaler = StandardScaler()
          splits = ['train','id_challenge','id_inhouse','ood_mosmed'] #['train_gc','train_fra','gc_gc','gc_frank','gc_mosmed','gc_radio']['id_train','id_test','ood_test']['train','gc_gc','gc_frank','gc_mosmed','gc_radio']['id_train','id_test','ood_test']['other'] 
          X,y,paths_pred = extract_features_train_id_od(filter_feature_extr,splits,used_feat)
  
          X_train = scaler.fit_transform(X[0])
          y_train = np.array(y[0])
          # #START of experiments with classification task
  
  # #################conf_ma class#################################################################################################################################        
  #         svm = SVC(class_weight='balanced')
  #         lr = LogisticRegression(class_weight='balanced')
  #         models = [svm,lr]
  #         model_names = ['SVM','Logistic_Regression']
  
  #         for model, model_name in zip(models,model_names):
  #             model.fit(X_train,y_train)
  #             for i,split in enumerate(splits):
  
  #                 X_eval = scaler.transform(X[i])
  #                 y_eval = y[i]
  
  #                 plot_confusion_matrix(model,X_eval,y_eval,labels=[1,2,3,4,5])
  #                 plt.title('Model {} on split {}'.format(model_name,split))
  #                 plt.savefig(os.path.join('storage','Results','confusion_matrixes',model_name+' '+split))
  #                 plt.show()
  #                 # y_pred = model.predict(X_eval)
  
  #                 # conf_matrix = confusion_matrix(y_eval,y_pred,labels=[1,2,3,4,5])
  #                 # if i != 0:
  #                 #     print(model_name,split)
  #                 #     print(conf_matrix)
  # #################conf_ma class#################################################################################################################################  
  
  # ##################3 feat models lr################################################################################################################################
  #         lr = LogisticRegression(class_weight='balanced')
  #         lr.fit(X_train,y_train)
  #         for i,split in enumerate(splits):
  #             X_eval = scaler.transform(X[i])
  #             y_eval = y[i]
  #             y_pred = lr.predict(X_eval)
  #             ma = conf_matrix(y_eval,y_pred)
  #             sens = sensitivity(ma)
  #             spec = specificity(ma)
  #             print(split)
  #             print(sens, spec)
  # ##################3 feat models lr################################################################################################################################
          # # confidence_bootstrap(X_train,y_train,X,y,svm,lr,scaler,splits)
          # model = svm
          # model.fit(X_train,[get_class(tru) for tru in y_train])
          # print(model)
          # for i,split in enumerate(splits):
          
          #     X_eval = scaler.transform(X[i])
          #     y_eval = y[i]
          #     y_gt = y_eval
          #     y_eval = np.array([get_class(tru) for tru in y_eval])
  
          #     y_pred =model.predict(X_eval)
          #     y_pred = get_float(y_pred)
              
          #     print(split)
  
          #     print('err :{}'.format(l1_loss(y_gt,y_pred,std=False)))
          #     print('overest errr :{}'.format(l1_loss_overestimation(y_gt,y_pred)))
          # model = lr
          # model.fit(X_train,y_train)
          # print(model)
          # for i,split in enumerate(splits):
          
          #     X_eval = scaler.transform(X[i])
          #     y_eval = y[i]
  
          #     y_pred =model.predict(X_eval)
          #     print(split)
          #     conf_ma = confusion_matrix(y_eval,y_pred,labels=[1,2,3,4,5])
          #     # print('sensitivity :{}'.format(micro_sens(conf_ma)))
          #     print('specificity :{}'.format(micro_spec(conf_ma)))
          # print('oberved, p-values')
          # print('LR')
          # lr.fit(X_train,y_train)
          # for i,split in enumerate(splits):
          
          #     X_eval = scaler.transform(X[i])
          #     y_eval = y[i]
  
          #     y_lr =lr.predict(X_eval)
          #     print(split)
          #     ma = conf_matrix(y_eval,y_lr)
          #     # print(ma)
          #     print('sensitivity :{},  {}'.format(sensitivity(ma), p_value_sensitivity(ma)))
          #     print('specificity :{},  {}'.format(specificity(ma), p_value_specificity(ma)))
          # print()
          # print('SVM')
          # svm.fit(X_train,y_train)
          # for i,split in enumerate(splits):
          
          #     X_eval = scaler.transform(X[i])
          #     y_eval = y[i]
  
          #     y_svm =svm.predict(X_eval)
          #     print(split)
          #     ma = conf_matrix(y_eval,y_svm)
          #     #print(ma)
          #     print('sensitivity :{},  {}'.format(sensitivity(ma), p_value_sensitivity(ma)))
          #     print('specificity :{},  {}'.format(specificity(ma), p_value_specificity(ma)))
          ##  END of classification tests 
  
          # #START of part for regression task
          # plot_variable_influence(X_train,X,y,splits)
          # paper_figures_by_split(X_train,X,y,splits)
          # find_ex_pred(X,y,paths_pred)
  
  ############conf_ma regression####################################################################################################       
          from sklearn.metrics._plot.confusion_matrix import ConfusionMatrixDisplay
  
          ridge = Ridge(normalize=False)
          svr = SVR()
          mlp = MLPRegressor((50,100,100,50))
          models = [ridge,svr,mlp]
          model_names = ['RidgeRegression','SVR','MLP']
  
          for model, model_name in zip(models,model_names):
              model.fit(X_train,y_train)
              for i,split in enumerate(splits):
                 
                  X_eval = scaler.transform(X[i])
                  y_eval = y[i]
  
                  y_pred_cont = model.predict(X_eval)
  
                  #In order to get a confusion matrix, we target and prediction values need to be binned:
                  y_eval = [get_class(truth) for truth in y_eval]
                  y_pred = [get_class(pred) for pred in y_pred_cont]
                  conf_matrix = confusion_matrix(y_eval,y_pred,labels=[1,2,3,4,5])
                  print(conf_matrix)
                  disp = ConfusionMatrixDisplay(conf_matrix,display_labels=[1,2,3,4,5])
                  disp.plot()
                  plt.title('Model {} on split {}'.format(model_name,split))
                  plt.savefig(os.path.join('storage','Results','confusion_matrixes',model_name+' '+split))
                  plt.show()
                  # print(model_name,split)
                  # print(conf_matrix)
  # ###################################################################################################################################
          
  
          # ridge.fit(X_train,y_train)
          # svr.fit(X_train,y_train)
          # mlp.fit(X_train,y_train)
          # for i,split in enumerate(splits):
          #     print()
          #     print(split)
          #     X_eval = scaler.transform(X[i])
          #     y_eval = y[i]
          #     print(split,'{:.3f} {:.3f} '.format(np.mean(y_eval),np.std(y_eval)))
  
              # y_ridge = ridge.predict(X_eval)
              # y_svr = svr.predict(X_eval)
              # y_mlp = mlp.predict(X_eval)
              
              # # # sens_r = sensitivity(conf_matrix([get_class(truth) for truth in y_eval],[get_class(pred) for pred in y_ridge]))
              # # # sens_s = sensitivity(conf_matrix([get_class(truth) for truth in y_eval],[get_class(pred) for pred in y_svr]))
              # # # sens_m = sensitivity(conf_matrix([get_class(truth) for truth in y_eval],[get_class(pred) for pred in y_mlp]))
  
              # print('minimal pred values {},{},{}'.format(np.min(y_ridge),np.min(y_svr),np.min(y_mlp)))
              # # ma_r = confusion_matrix([get_class(truth) for truth in y_eval],[get_class(pred) for pred in y_ridge],labels=[1,2,3,4,5])
              # # ma_s = confusion_matrix([get_class(truth) for truth in y_eval],[get_class(pred) for pred in y_svr],labels=[1,2,3,4,5])
              # # ma_m = confusion_matrix([get_class(truth) for truth in y_eval],[get_class(pred) for pred in y_mlp],labels=[1,2,3,4,5])
  
              # # # for matrix, model in zip([sens_r,sens_s,sens_m],['ridge','svr','mlp']):
              # # #     print(model)
              # # #     print('sensitivity :{}'.format(matrix))
                  
  
  
          #     ridge_err,ridge_std = l1_loss(y_ridge,y_eval, std = True)
          #     svr_err, svr_std  = l1_loss(y_svr,y_eval,std = True)
          #     mlp_err, mlp_std = l1_loss(y_mlp,y_eval,std = True)
  
          #     ridge_err_over,ridge_std_over = l1_loss_overestimation(y_ridge,y_eval, std = True)
          #     svr_err_over, svr_std_over  = l1_loss_overestimation(y_svr,y_eval,std = True)
          #     mlp_err_over, mlp_std_over = l1_loss_overestimation(y_mlp,y_eval,std = True)
  
          #     # l1_loss_bins(y_svr,y_eval,split)
  
          #     for i,std in enumerate([ridge_std,svr_std,mlp_std]):
          #         stds_of_splits[i].append(std)
          #     for i,std in enumerate([ridge_std_over,svr_std_over,mlp_std_over]):
          #         stds_of_splits_over[i].append(std)
          #     for i,errors in enumerate([get_l1_losses(y_ridge,y_eval),get_l1_losses(y_svr,y_eval),get_l1_losses(y_mlp,y_eval)]):
          #         for err in errors:  
          #             all_errors[i].append(err)
          #     for i in range(3):
          #         all_errors_over[i] = [max(err,0) for err in all_errors[i]]
  
          #     #a vector that predicts the mean value of y_train, is a baseline
          #     # y_mean = np.mean(y_train)*np.ones(np.shape(y_eval))
          #     print('{}    :        ridge               svr             mlp  '.format(split))
          #     print(u"error          {:.3f} \u00B1 {:.3f},   {:.3f} \u00B1 {:.3f},  {:.3f} \u00B1 {:.3f}".format(ridge_err,ridge_std,svr_err,svr_std,mlp_err,mlp_std))
          #     print(u"error overest. {:.3f} \u00B1 {:.3f},   {:.3f} \u00B1 {:.3f},  {:.3f} \u00B1 {:.3f}".format(ridge_err_over,ridge_std_over,svr_err_over,svr_std_over,mlp_err_over,mlp_std_over))
  
          #     # print('Using mean of train values, has error of {} and std of {}'.format(l1_loss(y_mean,y_train,std=True)[0],l1_loss(y_mean,y_train,std=True)[1]))
          #     print()
  
          # print('TOTAL:                       ridge     svr      mlp  ')
          # print('mean of split stds          {:.3f},   {:.3f},  {:.3f}'.format(np.mean(stds_of_splits[0]),np.mean(stds_of_splits[1]),np.mean(stds_of_splits[2])))
          # print('std of all errors           {:.3f},   {:.3f},  {:.3f}'.format(np.std(all_errors[0]),np.std(all_errors[1]),np.std(all_errors[2])))
          # print('mean of split stds overest. {:.3f},   {:.3f},  {:.3f}'.format(np.mean(stds_of_splits_over[0]),np.mean(stds_of_splits_over[1]),np.mean(stds_of_splits_over[2])))
          # print('std of all errors overest.  {:.3f},   {:.3f},  {:.3f}'.format(np.std(all_errors_over[0]),np.std(all_errors_over[1]),np.std(all_errors_over[2])))
          # # END regression task
  
  if __name__ == "__main__":
      #(SET all params, depending on desired train procedurey)
      main([0,1,2,3],preprocessing=False,feature_extraction=False, extract_dice_scores=False,model_train=True)
  
  
  ## Examples for images to show 
  # for dices : copies: img:  Task100_RadiopediaTrain_0009 pred:  Task740_ChallengeTrainF4'
  # for int mode : not lung tissue: img : Task200_MosmedTrain_0033 pred : Task740_ChallengeTrainF4'
  # for conn comp : too muh and many: img : Task100_RadiopediaTrain_0008 pred : Task541_FrankfurtTrainF4'
  #                                or img : Task100_RadiopediaTrain_0007 pred : Task541_FrankfurtTrainF4'
  # seg in lung 