Feature Request: Automatic K optimization for Fourier Terms

If I wish to fit a regression with Fourier terms then to find the optimal K I need to do something like this:

library(fable)
library(dplyr)
library(tidyr)

mbl = tsibbledata::ansett %>%
  tsibble::fill_gaps() %>%
  model(arima1 = ARIMA(Passengers ~ fourier(K = 1) + PDQ(0,0,0)),
        arima2 = ARIMA(Passengers ~ fourier(K = 2) + PDQ(0,0,0)),
        arima3 = ARIMA(Passengers ~ fourier(K = 3) + PDQ(0,0,0)))

metrics = mbl %>%
  glance()

mbl_best = metrics %>%
  select(Airports, Class, .model, AICc) %>%
  group_by(Airports, Class) %>%
  slice(which.min(AICc)) %>%
  left_join(mbl %>%
              gather('.model', 'model', -Airports, -Class),
            by = c('.model', 'Airports', 'Class')) %>%
  as_mable(key = c('Airports', 'Class'), models = 'model')

It would be more convenient for K to be automatically determined through something like this:

model(arima = ARIMA(Passengers ~ Fourier(K = 1:3) + PDQ(0,0,0)

On that note, when I look at the source code for ARIMA it appears that when fitting a regression + ARIMA the number of differences is determined after the regression. Because of this, it seems entirely possible that the arima1, arima2 and arima3 models I fit could potentially have a different number of differencing. If this is indeed the case perhaps determining K through cross validation is better?

Thanks!

The payload names no files or tests; begin at the ARIMA entry point and existing fourier(K = ...) handling. Resolve whether automatic K selection should use information criteria or cross-validation and how differencing interacts with it, then verify the agreed behavior using the ansett example.