Data adapter error for keras functional api in R - multi input time series data

[toastershock]

I want to run a keras model in R with two time series inputs and one output. But when I am trying to fit the model, I get an error which is related to the multi input and using the function timeseries_dataset_from_array function.

Error in py_call_impl(callable, call_args$unnamed, call_args$named) : 
  ValueError: Failed to find data adapter that can handle input: (<class 'list'> containing values of types {"<class 'tensorflow.python.data.ops.batch_op._BatchDataset'>"}), <class 'NoneType'>
Run `reticulate::py_last_error()` for details.
> reticulate::py_last_error()

── Python Exception Message ────────────────────────────────────────────────────────────────────
Traceback (most recent call last):
  File "C:\Users\tburing\DOCUME~1\VIRTUA~1\R-TENS~1\lib\site-packages\keras\src\utils\traceback_utils.py", line 70, in error_handler
    raise e.with_traceback(filtered_tb) from None
  File "C:\Users\tburing\DOCUME~1\VIRTUA~1\R-TENS~1\lib\site-packages\keras\src\engine\data_adapter.py", line 1105, in select_data_adapter
    raise ValueError(
ValueError: Failed to find data adapter that can handle input: (<class 'list'> containing values of types {"<class 'tensorflow.python.data.ops.batch_op._BatchDataset'>"}), <class 'NoneType'>

── R Traceback ─────────────────────────────────────────────────────────────────────────────────
    ▆
 1. ├─model %>% ...
 2. ├─generics::fit(...)
 3. └─keras:::fit.keras.engine.training.Model(...)
 4.   ├─base::do.call(object$fit, args)
 5.   └─reticulate (local) `<python.builtin.method>`(...)
 6.     └─reticulate:::py_call_impl(callable, call_args$unnamed, call_args$named)

Here is the reproducable example:

library(keras)

data <- list(
  array(rnorm(1152*52*53*2), dim = c(1152,52,53,2)),
  array(rnorm(1152*22*26), dim = c(1152,22,26))
)
# Create a random sample
indices<-1:n
n <- nrow(data[[1]]) # assuming all your data lists have the same number of rows

train_indices <- 1:((n/6)*4)
val_indices <- ((n/6)*4+1):((n/6)*5)
test_indices <- ((n/6)*5+1):n

# Split the data
initial_phase_train <- data[[1]][train_indices, , , ]
initial_phase_val <- data[[1]][val_indices, , ,]
initial_phase_test <- data[[1]][test_indices, , ,]

previous_3weeks_train <- data[[1]][train_indices, , ,]
previous_3weeks_val <- data[[1]][val_indices, , ,]
previous_3weeks_test <- data[[1]][test_indices, , ,]

target_data_train <- data[[2]][train_indices, , ]
target_data_val <- data[[2]][val_indices, , ]
target_data_test <- data[[2]][test_indices, , ]


initial_phase_train   <-  timeseries_dataset_from_array(initial_phase_train, 
                                                        target_data_train, 
                                                        sampling_rate = 1,
                                                        sequence_length = 25,
                                                        shuffle = FALSE,
                                                        batch_size = 32)

initial_phase_val     <-  timeseries_dataset_from_array(initial_phase_val, 
                                                        target_data_val, 
                                                        sampling_rate = 1,
                                                        sequence_length = 25,
                                                        shuffle = FALSE,
                                                        batch_size = 32)

previous_3weeks_train <- timeseries_dataset_from_array(previous_3weeks_train, 
                                                       target_data_train, 
                                                       sampling_rate = 1,
                                                       sequence_length = 3,
                                                       shuffle = FALSE,
                                                       batch_size = 32)

previous_3weeks_val   <- timeseries_dataset_from_array(previous_3weeks_val, 
                                                       target_data_val, 
                                                       sampling_rate = 1,
                                                       sequence_length = 3,
                                                       shuffle = FALSE,
                                                       batch_size = 32)


# Define the model with two inputs
input_initial <- layer_input(shape = c(25, 52, 53, 2), name =  "Initial Phase")
input_previous <- layer_input(shape = c(3, 52, 53, 2), name =  "2nd Phase (prev. 3 weeks)")

# LSTM layers for each input
inputs <- layer_concatenate(list(input_initial, input_previous), axis = -4)

outputs <- inputs %>%
  layer_conv_lstm_2d(
    filters = 32, 
    kernel_size = c(3, 3), 
    padding = "same",
    return_sequences = TRUE
  ) %>%
  layer_dropout(0.5) %>%
  layer_conv_lstm_2d(
    filters = 32, 
    kernel_size = c(3, 3), 
    padding = "same",
    return_sequences = FALSE
  ) %>%
  layer_dropout(0.5) %>%
  layer_max_pooling_2d(pool_size = c(2, 2)) %>%
  layer_flatten() %>%
  layer_dense(units = 22 * 26 * 1) %>%
  layer_reshape(target_shape = c(1, 22, 26, 1)) %>%
  layer_conv_2d(
    filters = 1, 
    kernel_size = c(3, 3), 
    padding = "same",
    activation = 'linear'
  )

# Create the model
model <- keras_model(inputs = c(input_initial, input_previous), outputs = outputs)

# Compile the model
model %>% compile(
  optimizer = "adam",# optimizer_adam(learning_rate = 1e-4),
  loss = 'mean_absolute_error',
  metrics = "mean_absolute_error"
)

# Fit the model
history <- model %>% fit(
  x = list(initial_phase_train, previous_3weeks_train),
  validation_data = list(initial_phase_val, previous_3weeks_val),
  epochs = 4
)

Any help is very much appreciated!
Many thanks!