from urllib.request import install_opener
import numpy as np
import matplotlib.pyplot as plt
import copy

import tensorflow as tf
from tensorflow.keras.models import Model


# ------------------------------------------------------------------------------------------------------
# The class for the specifications of the model:
class SpecsClass:
  """
      n_ins: Number of inputs
      n_lays: Number of layers
      n_outs: Number of outputs
  """
  def __init__(self, num_ins, num_lays, num_outs):
    self.n_ins = num_ins
    self.n_lays = num_lays
    self.n_outs = num_outs

# ------------------------------------------------------------------------------------------------------
# Function to retrieve the information of the model:
def study_model(main_model):
  """
      Loop through the layers of the base model, sequentially create new models,
      get all the info of the activation funcs, get the arrays of weights and biases:
  """

  # The main model:
  main_model_layers = main_model.layers

  # The way to obtain the intermediate outputs of the network is via having sub-modules for each layer:
  arr_submodels = [ Model(layer.input , layer.output) for i, layer in enumerate(main_model_layers) ]

  # Number of layers, number of inputs, number of outputs:
  num_lays = len(main_model_layers)
  num_ins = main_model_layers[0].input_shape[1]
  num_outs = main_model_layers[-1].output_shape[1]

  # Initialize the specs object:
  specs = SpecsClass(num_ins, num_lays, num_outs)

  # Array with the string names of the activation functions:
  arr_types_activfuncs = [ tf.keras.activations.serialize(main_model_layers[i].activation) for i in range(specs.n_lays) ]

  # Get the arrays of weights and biases:
  arr_ws = np.array([ layer.get_weights()[0].transpose() for layer in main_model_layers ], dtype=object)
  arr_bs = np.array([ layer.get_weights()[1] for layer in main_model_layers ], dtype=object)

  return specs, arr_types_activfuncs, arr_ws, arr_bs


# ------------------------------------------------------------------------------------------------------
# Get the rowvec_m_x_w:
# rowvec_m_x_w is a vector of slopes (m) times (product symbol x) the weights (w):
def apply_m_x_w(mat_slopes, colvec_weights):
  return np.sum(mat_slopes * colvec_weights, 0)

# Get the val_n_x_w:
def apply_n_x_w(vec_indep_terms, colvec_weights):
  return np.sum(vec_indep_terms * colvec_weights)

# Function to apply the product:
def apply_uv(rowvec_m_x_w, indep_term, u, v):
    return u * rowvec_m_x_w, u * indep_term + v

# ------------------------------------------------------------------------------------------------------
# Derivatives of Activation Functions, returns u (dep. term) and v (indep. term) of the derivative (except softmax):


# Elu:
def af_elu(af_in, af_out, alpha=1):
  if af_in<0:
    u = af_out + alpha
    v = af_out - u*af_in
  else:
    u, v = 1, 0
  return u, v


# Exponential:
def af_exp(af_in, af_out):
  u = af_out
  v = af_out - u*af_in
  return u, v


# Gelu:
def af_gelu(af_in, af_out):
  u = af_out/af_in + (af_in/(np.sqrt(2*np.pi))) * np.exp( - (af_in**2)/2 )
  v = af_out - u*af_in
  return u, v


# Hard Sigmoid:
def af_hardsigmoid(af_in, af_out):
  if af_in<-2.5:
    u, v = 0, 0
  elif 2.5<af_in:
    u, v = 0, 1
  else:
    u, v = 0.2, 0.5
  return u, v


# Linear:
def af_linear(af_in, af_out):
  u, v = 1, 0
  return u, v


# Relu:
def af_relu(af_in, af_out, alpha=0, max_value=np.inf, threshold=0):
  if af_out >= max_value:
    return 0, max_value

  if threshold < af_in:
    u, v = 1, 0
  else:
    u, v= alpha, -alpha*threshold

  return u, v


# Selu:
def af_selu(af_in, af_out):

  # alpha and scale are pre-defined constants:
  alpha=1.67326324
  scale=1.05070098

  if af_in<0:
    u = scale*alpha*np.exp(af_in)
    v = af_out - u*af_in
  else:
    u, v = scale, 0
  return u, v


# Sigmoid:
def af_sigmoid(af_in, af_out):
  u = af_out*(1-af_out)
  v = af_out - u*af_in
  return u, v


# Softmax (not normal af, uses info of all the activation functions in the layer):
def afs_softmax(mat_m_x_w, rowvec_afs_outs, rowvec_x):

  # Transform the outs of the activation functions to a column vector:
  colvec_afs_outs = rowvec_afs_outs.reshape((-1, 1))

  # Get the matrix of slopes for the outputs:
  sum_of_prod_1 = np.sum(np.multiply(mat_m_x_w, colvec_afs_outs), axis =0)
  diff_1 = np.subtract(mat_m_x_w, sum_of_prod_1)
  new_mat_slopes = np.multiply(diff_1, colvec_afs_outs)

  # Get the vector of independent terms for the outputs:
  sum_of_prod_2 = np.sum(np.multiply(new_mat_slopes, rowvec_x), axis=1).reshape((-1, 1))
  new_colvec_indeps = np.subtract(colvec_afs_outs, sum_of_prod_2)

  return new_mat_slopes, new_colvec_indeps


# Softplus:
def af_softplus(af_in, af_out):
  u = 1/(1+np.exp(-af_in))
  v = af_out - u*af_in
  return u, v


# Softsign:
def af_softsign(af_in, af_out):
  u = (1+np.abs(af_in))**(-2)
  v = af_out - u*af_in
  return u, v


# Swish:
def af_swish(af_in, af_out):
  u = (af_out/af_in) * (1+af_in-af_out)
  v = af_out - u*af_in
  return u, v


# Tanh:
def af_tanh(af_in, af_out):
  u = (np.cosh(af_in))**(-2)
  v = af_out - u*af_in
  return u, v


# ------------------------------------------------------------------------------------------------------
# Function to decide which activation function to use:
def get_uv(af_type, af_in, af_out):

  if af_type == 'elu':
    return af_elu(af_in, af_out) # , alpha)

  elif af_type=='exponential':
    return af_exp(af_in, af_out)

  elif af_type=='gelu':
    return af_gelu(af_in, af_out)

  elif af_type=='hard_sigmoid':
    return af_hardsigmoid(af_in, af_out)

  elif af_type=='linear':
    return af_linear(af_in, af_out)

  elif af_type=='relu':
    return af_relu(af_in, af_out) # , alpha, max_value, threshold)

  elif af_type=='selu':
    return af_selu(af_in, af_out)

  elif af_type=='sigmoid':
    return af_sigmoid(af_in, af_out)

  elif af_type=='softplus':
    return af_softplus(af_in, af_out)

  elif af_type=='softsign':
    return af_softsign(af_in, af_out)

  elif af_type=='swish':
    return af_swish(af_in, af_out)

  elif af_type=='tanh':
    return af_tanh(af_in, af_out)

  # Allow users to have their own activation functions if they provide the code and the differentiation?
  # For now its okay.
  else:
    print(f"The activation function {af_type} was not found")

# ------------------------------------------------------------------------------------------------------
# This is the fundamental function:
def core_algorithm(specs,  arr_afs, arr_afs_ins, arr_afs_outs, arr_ws, arr_bs, rowvec_x):
  """
      Bias:  what the nn model has in each neuron, trained variables.
      Independent term: what is carried from before, part of the explanation.

      specs: Contains the number of input dimensions, layers and output dimensions.

      arr_afs: Types of activation functions.
      arr_afs_ins: Same as the arr_entries_to_afs, inputs to the activation functions, after applying weights.
      arr_afs_outs: The outputs of the activation functions, i.e., outputs of the layers.
      arr_ws: Same as arr_weights
      arr_bs: Same as arr_biases

      rowvec_x: the input.
  """

  ranger_layers = range(specs.n_lays)

  # USE DICTIONARIES FOR arr_slopes and arr_indeps TO MAKE THINGS FASTER!!
  # DO THE SAME IN ALL THE FILE!!!

  # Initialize:
  arr_slopes = [ [] for _ in range(specs.n_lays+1)]
  arr_indeps = [ [] for _ in range(specs.n_lays+1)]

  # Initialize the first entry:
  arr_slopes[0] = np.array(np.matrix(np.eye(specs.n_ins)))
  arr_indeps[0] = np.array([0 for _ in range(specs.n_ins) ]).reshape(specs.n_ins,1)

  # Go through the layers:
  for layer, rowvec_afs_ins, rowvec_afs_outs, mat_weights, rowvec_biases in zip(ranger_layers, arr_afs_ins, arr_afs_outs, arr_ws, arr_bs):

    mat_slopes = arr_slopes[layer]
    vec_indep_terms = arr_indeps[layer]

    # If the activation function is not softmax:
    if arr_afs[layer] != 'softmax':

      # Loop through the neurons:
      for af_in, af_out, rowvec_weights, bias in zip(rowvec_afs_ins, rowvec_afs_outs, mat_weights, rowvec_biases):

        # Convert it to a column vector:
        colvec_weights = np.array(rowvec_weights).reshape((len(rowvec_weights),1))

        # Weight the slopes:
        rowvec_m_x_w = apply_m_x_w(mat_slopes, colvec_weights)

        # Weight the independent term:
        val_n_x_w = apply_n_x_w(vec_indep_terms, colvec_weights)
        # Add the bias:
        indep_and_bias = bias + val_n_x_w

        # Differentiate the activation function and get u,v:
        u, v = get_uv(arr_afs[layer], af_in, af_out)

        # Apply the u and v values from the differentiation of the activation function:
        new_rowvec_slopes, new_indep_term = apply_uv(rowvec_m_x_w, indep_and_bias, u, v)

        # Save:
        arr_slopes[layer+1].append( new_rowvec_slopes )
        arr_indeps[layer+1].append( new_indep_term )

      # Convert to np arrays:
      arr_slopes[layer+1] = np.array(np.matrix(arr_slopes[layer+1]))
      arr_indeps[layer+1] = np.array(arr_indeps[layer+1]).reshape((len(arr_indeps[layer+1]),1))

    # If the activation function is softmax:
    else:

      # Modify the mat_slopes.
      # The rows of mat_slopes represent the linear funcs of each neuron of the previous layer.
      # The rows of new_mat_slopes represent the linear funcs of each neuron of the current layer.

      # Reshape to make it a 3D mat:
      curr_shp = mat_weights.shape
      rshp_tp_mat_weights = np.reshape(mat_weights, (curr_shp[0], curr_shp[1], 1))
      # Multiply the rshp_tp_mat_weights times the mat_slopes:
      prod_m_x_w = np.multiply(mat_slopes, rshp_tp_mat_weights)
      # Get the sum of the slices:
      mat_m_x_w = np.sum(prod_m_x_w, axis=1)

      # Get the output slopes matrix and the output indep terms vec:
      new_mat_slopes, new_colvec_indeps = afs_softmax(mat_m_x_w, rowvec_afs_outs, rowvec_x)

      # Update arr_slopes[layer+1] and arr_indeps[layer+1]:
      arr_slopes[layer+1] = new_mat_slopes
      arr_indeps[layer+1] = new_colvec_indeps


  # No need to return the entire array of intermediate slopes and indep terms.
  # Only the array with last layer's neurons (there might be more than one output):
  return arr_slopes[-1], arr_indeps[-1]


# ------------------------------------------------------------------------------------------------------
# Function to get the contributions of the inputs:
def get_contributions(input_i, mat_LF_dep_terms_i, specs):

  # Create an array with each individual absolute value of slope*input:
  arr_abs_slopes_x_input = np.absolute(np.multiply(input_i, mat_LF_dep_terms_i))

  # Get the sum:
  aux_sum_of_arr = np.sum(arr_abs_slopes_x_input, axis=1).reshape(specs.n_outs,1)

  # Finally the contributions for all the outputs:
  expl_contribs = arr_abs_slopes_x_input/aux_sum_of_arr

  return expl_contribs


# ------------------------------------------------------------------------------------------------------
# Function to get the output predicted with the LF:
def get_LF_output(input_i, mat_LF_dep_terms_i, mat_LF_indep_term_i):

    # Get the output with all its dimensions. Reshape colvec mat_LF_indep_term_i so its a row vector.
    return np.sum(input_i*mat_LF_dep_terms_i, axis=1) + mat_LF_indep_term_i.flatten()


# ------------------------------------------------------------------------------------------------------
# Function to retrieve the adjective within the ranges record using a value:
def get_adj_from_val(val, ranges):

  # A A A A this shouldn't be done, store it as a class so you have it already:
  locs = [ entry[0] for entry in ranges]
  names =[ entry[1] for entry in ranges]


  # Get where it fits the val:
  for ii in range(len(locs)-1):
      if locs[ii]<=val and val<locs[ii+1]:
          return names[ii]

  # If we reach this point is the last:
  return names[-1]

# ------------------------------------------------------------------------------------------------------
# Function to get the final text explanation:
def get_text_expl(inst, arr_expl_slopes, arr_expl_contribs):

  ranger = range(len(inst))

  # Our own baselines, rethink how they are saved A A A A A:
  norm_var_baseline = [[0, 'Very Low'], [0.2, 'Low'], [0.4, 'Average'], [0.6, 'High'], [0.8, 'Very High']]
  scld_var_baseline = [[-np.Inf, 'Very Low'], [-3, 'Low'], [-1, 'Average'], [1, 'High'], [3, 'Very High']]

  slopes_baseline = [[-np.Inf, 'Very Low'], [-3, 'Low'], [-1, 'Average'], [1, 'High'], [3, 'Very High']]
  contribs_baseline = [[0, 'Very Low'], [0.2, 'Low'], [0.4, 'Average'], [0.6, 'High'], [0.8, 'Very High']]

  # We construct for each input, using the information from the user (here we use the baselines):
  ranges_inputs = [ norm_var_baseline for _ in ranger ]
  ranges_slopes = [ slopes_baseline for _ in ranger ]
  ranges_contribs = [ contribs_baseline for _ in ranger ]
  # Also the names of the variables:
  names_inputs = [ ' Variable ' + str(i+1) +' ' for i in ranger ]

  # The start of the 3 sentences:
  expl_txt_slopes = " For this instance, the output variable depends on: "
  expl_txt_contribs = " For this instance, the contributions of the variables are: "
  expl_txt_input = " Because the input has values of: "

  # Initialize the text_expl variable, each entry refers to one output dimension:
  # Inside each of these output dimensions, we have the 3 explanations.
  arr_text_expl = []

  for expl_slopes, expl_contribs in zip(arr_expl_slopes, arr_expl_contribs):

    # Go through the slopes:
    for val, ranges, name in zip(expl_slopes, ranges_slopes, names_inputs):
      txt_to_add = get_adj_from_val(val, ranges)
      expl_txt_slopes = expl_txt_slopes + name + txt_to_add + ','

    # Go through the contributions:
    for val, ranges, name in zip(inst, ranges_contribs, names_inputs):
      txt_to_add = get_adj_from_val(val, ranges)
      expl_txt_contribs = expl_txt_contribs + name + txt_to_add + ','

    # Go through the input:
    for val, ranges, name in zip(expl_contribs, ranges_inputs, names_inputs):
      txt_to_add = get_adj_from_val(val, ranges)
      expl_txt_input = expl_txt_input + name + txt_to_add + ','

    # Replace with a dot the last entry:
    expl_txt_input = expl_txt_input[:-1]+'.'
    expl_txt_contribs = expl_txt_contribs[:-1]+'.'
    expl_txt_input = expl_txt_input[:-1]+'.'

    # Append this entry:
    arr_text_expl.append([expl_txt_slopes, expl_txt_contribs, expl_txt_input])

  return arr_text_expl


# ------------------------------------------------------------------------------------------------------
# The class for the output explanations of the request:
class ExplClass:
  """
      The class for the output explanations of the request
  """
  def __init__(self):
    # They will all be dictionaries:
    self.input_values = {}
    self.output_pred_NN = {}
    self.output_pred_LF = {}
    self.LF_dep_terms = {}
    self.LF_indep_term = {}
    self.input_contribs = {}
    self.text_explanation = {}
    self.reference_id = {}

  # When adding a new instance processed:
  def add_info(self, index_i, input_i, output_pred_NN_i, output_pred_LF_i,
               mat_LF_dep_terms_i, mat_LF_indep_term_i, mat_input_contribs_i, arr_text_explanation_i, reference_id):
    # Assing:
    # In the assignation we correct the shape of arr_LF_indep_term_i, mat_LF_indep_term_i, mat_input_contribs_i
    # to respect the criterion: rows instances, cols dimensions (both the input and outputs are row vecs)
    self.input_values[index_i] = input_i
    self.output_pred_NN[index_i] = output_pred_NN_i.tolist()
    self.output_pred_LF[index_i] = output_pred_LF_i.tolist()
    self.LF_dep_terms[index_i] = mat_LF_dep_terms_i.tolist()
    self.LF_indep_term[index_i] = mat_LF_indep_term_i.tolist()
    self.input_contribs[index_i] = mat_input_contribs_i.tolist()
    self.text_explanation[index_i] = arr_text_explanation_i
    self.reference_id[index_i] = reference_id


# ------------------------------------------------------------------------------------------------------
# Get the explanations function:
def run_explanations(X, model_l, specs, arr_afs, arr_ws, arr_bs):

  # Initialize the structure:
  generated_explanations = ExplClass()

  # The computationally most efficient approach:
  #    -  A single model that outputs all intermediate outputs at once.
  #    -  Calling this model only once for all the instances of the request.

  # Get the reference for the outputs of all the layers:
  ref_outs_all_layers = [layer.output for layer in model_l.layers]
  # Construct the new model that returns the outputs of all the layers:
  model_of_outs_all_layers  = Model(inputs=model_l.input, outputs=ref_outs_all_layers)
  # Feed the model with all the instances of the request:
  # Is vital to process all instances together. Much less computational time than calling the model iteratively:
  outs_all_layers = model_of_outs_all_layers.predict(X)


  # Loop through the instances of the submission:
  for i, input_i in enumerate(X):

    # We need the interm. outs of every instance separately. We use 4 arrays with different formats.
    # These 3 arrays have repeated info in the middle (which represents the interm. outs),
    # the left (lft) & right (rgt) sides tell us if we consider the input (lft) or the model output (rgt).

    # outs_mid_n_rgt: Mid and right, an array that has this instance's outputs of all the layers.
    # outs_mid_n_lft: Mid and left, an array that has this instance's outputs of interm layers and the input.
    # outs_mid: Mid, an array that has only this instance's outputs of interm layers.

    outs_mid_n_rgt = [outs_all_layers[lay][i] for lay in range(specs.n_lays)]
    outs_mid = outs_mid_n_rgt[:-1]
    outs_mid_n_lft = [input_i] + outs_mid


    # Array of entries for the activation functions:
    # ***** TODO: If the layer is not dense is different:
    arr_afs_ins_i = [ np.sum(arr_ws[j] * outs_mid_n_lft[j], 1) + arr_bs[j] for j in range(specs.n_lays)]

    # Note about:  mat_LF_dep_terms_i,  mat_LF_indep_term_i,  mat_input_contribs_i
    #   Their shape is different from the criterion we are following, for faster numpy ops and visibility.
    #   mat_LF_dep_terms_i is a 2D mat with rows that identify output dimensions, and cols the terms.
    #   The criterion is that the input is a row vector, and the output is also a row vector.
    #   However, when we save in the ExplClass, we will convert to lists, so we will respect criterion.

    # Run the frontpropagation algorithm, mat_LF_dep_terms_i & mat_LF_indep_term_i are both arrays
    # of arrays with the terms for all the outputs (one array of terms per output).
    mat_LF_dep_terms_i, mat_LF_indep_term_i = core_algorithm(specs, arr_afs, arr_afs_ins_i, outs_mid_n_rgt,
                                                             arr_ws, arr_bs, input_i)

    # Get the mat of arrays with the input contributions for each output:
    mat_input_contribs_i = get_contributions(input_i, mat_LF_dep_terms_i, specs)

    # Get the vector of outputs predicted with the LFs for each output dimension:
    output_pred_LF_i = get_LF_output(input_i, mat_LF_dep_terms_i, mat_LF_indep_term_i)

    # Get the text explanations:
    arr_text_explanation_i = get_text_expl(input_i, mat_LF_dep_terms_i, mat_input_contribs_i)

    # Save the information of the instance:
    generated_explanations.add_info(i, input_i, outs_mid_n_rgt[-1], output_pred_LF_i,
                                    mat_LF_dep_terms_i, mat_LF_indep_term_i, mat_input_contribs_i,
                                    arr_text_explanation_i, i)

  return generated_explanations



# ------------------------------------------------------------------------------
# Visualization of the linear function prediction vs NN prediction:
# Low level function to compare the predictions of a linear function prediction vs NN prediction:
def plot_LF_vs_NN_predictions(inst, arr_LF_deps, arr_LF_indep, arr_out_pred_LF, arr_out_pred_NN,
                              main_model, specs, M, proxim_thres, flag_normalization, which_out_dims):
  """
      inst: The instance of the request that we want to evaluate.

      arr_ : Each of the array entries represent a dimension.

      arr_LF_deps, arr_LF_indep: Dep and indep terms of the Linear Function approximation.

      arr_out_pred_LF, arr_out_pred_NN: Predicted output of the instance with the LF and the NN respectively.

      M: Number of points desired to plot in the vicinity of the instance.

      proxim_thres: Allowed maximum difference in every dimension to generate the sample.

      flag_normalization: Boolean to know if the input vector is normalized. If so, we make sure
                         that the random points generated are within the limits of 0 and 1.

      which_out_dims: To know which output dimesnions we want to plot. Could be "all" or
                      an array with the indexes of the out dims desired to plot.
                      For example [0,3,5]. Indexes start at 0.

  """

  # Convert to numpy and reshape for manipulation:
  inst = np.array(inst)
  arr_LF_deps = np.array(arr_LF_deps).reshape((specs.n_outs, 1, specs.n_ins))
  arr_LF_indep = np.array(arr_LF_indep).reshape((1, specs.n_outs))
  arr_out_pred_LF = np.array(arr_out_pred_LF)
  arr_out_pred_NN = np.array(arr_out_pred_NN)

  # Convert the variable which_out_dims to an array that serves as a flag:
  # If the user wants all:
  if which_out_dims =="all":
    arr_flags_plot_out_dims = [ True for a in range(specs.n_outs) ]
  # which_out_dims is an array with the dimensions:
  else:
    # Initialize:
    arr_flags_plot_out_dims = [ False for a in range(specs.n_outs) ]
    # Loop:
    for ind in which_out_dims:
      # Chack that it belongs to the possible outputs:
      if 0<=ind<specs.n_outs:
        # Assing:
        arr_flags_plot_out_dims[ind] = True


  # Generation of the random points in the neighborhood of the instance:

  # Initialize the random number generator:
  rng = np.random.default_rng()

  # Get the random normal sample of M new points. Rows, points, columns dimensions:
  numbers = rng.normal(size=(M, specs.n_ins))

  # Normalize each column/dimension:
  numbers = numbers/np.max(np.abs(numbers), axis=0)

  # Get the differences for the dimension:
  diffs = numbers * proxim_thres

  # Check if we are using normalized inputs.
  # If that is the case, we need to make sure we are not exceeding the limits:
  if flag_normalization:

    # Add instance values to max and min differences of each dimension (top and bottom values):
    tops = np.max(diffs, axis=0) + inst
    bots = np.min(diffs, axis=0) + inst

    # Get the two filters of being off-limits:
    filter_top = tops>1.0
    filter_bot = bots<0.0

    # See where they exceed the normalized values in either tops or bots:
    where_offlimits = set( np.where(filter_top)[0].tolist() + np.where(filter_bot)[0].tolist() )

    # Go through those columns where we are off limits:
    for col in where_offlimits:

      # Now get the positivies and the negatives:
      filter_positives = diffs[:,col]>0
      filter_negatives = ~ filter_positives
      positives = diffs[:,col][filter_positives]
      negatives = diffs[:,col][filter_negatives]

      # Check if you exceed the top:
      if filter_top[col]:
        # Denorm and renorm the positives so that the max is 1-instance[col]
        positives = (positives/np.max(positives))*(1-inst[col])

      if filter_bot[col]:
        # Denorm and renorm the negatives so that the min is -instance[col]
        negatives = (negatives/abs(np.min(negatives))) * inst[col]

      # Now join them, shuffle:
      new_diffs_of_this_col = np.concatenate((positives, negatives))
      np.random.shuffle( new_diffs_of_this_col )

      # Assign it to the column in diffs:
      diffs[:,col] = new_diffs_of_this_col


  # Get the euclidean distance directly from the differences:
  distances = np.mean( np.square(diffs), axis=1).reshape((M, 1))

  # Apply these differences to generate the random points in the neighborhood:
  random_points = (diffs + inst).reshape((1, M, specs.n_ins))

  # Get the predictions with the NN for the random points. [0] needed for shape adjustment:
  Yrandom_pred_NN = main_model.predict(random_points)[0]

  # Get the predictions with the LF approximation of the instance for the random points:
  aux_1 = random_points*arr_LF_deps
  aux_2 = np.sum(aux_1, axis=2).transpose()
  Yrandom_pred_LF = aux_2 + arr_LF_indep

  # Loop to create the figures:
  for j in range(specs.n_outs):

    # Check if you want to plot this one:
    if arr_flags_plot_out_dims[j]:

      # Figure:
      plt.figure()
      # Perform the scatterplot, should be a curve tangent to the 1-1 line:
      # Color indicates the distance from the instance:
      scat_im = plt.scatter(Yrandom_pred_NN[:,j], Yrandom_pred_LF[:,j], c = distances, cmap="jet")
      clb = plt.colorbar(scat_im)
      clb.ax.set_title('Eucl. Dist.')
      # Get current limits of plot, and the ref lims:
      x_lims = plt.gca().get_xlim()
      y_lims = plt.gca().get_ylim()
      ref_lims = [np.min((x_lims[0], y_lims[0])), np.max((x_lims[1], y_lims[1])) ]
      # Show the 1-1 line:
      plt.plot( ref_lims, ref_lims, color='grey')
      # Plot two dashed lines to mark the location of the instance:
      plt.plot([ref_lims[0], arr_out_pred_NN[j]], [arr_out_pred_LF[j], arr_out_pred_LF[j]], linestyle='dashed', color='grey')
      plt.plot([arr_out_pred_NN[j], arr_out_pred_NN[j]], [ref_lims[0], arr_out_pred_LF[j]], linestyle='dashed', color='grey')
      # Show the instance:
      plt.scatter(arr_out_pred_NN[j], arr_out_pred_LF[j], s=25, c="black", marker="*")
      # Set the limits to the ref lims:
      plt.xlim(ref_lims)
      plt.ylim(ref_lims)
      # Details:
      plt.xlabel(f"Output {j+1} Predicted with NN")
      plt.ylabel(f"Output {j+1} Predicted with LF")
      plt.title(f"Rand. Points Near Instance - Out Dim. {j+1}")

  return











