The k-index measures the condition of the magnetosphere. It is usually averaged over three hour, so each day has 8 measurements.

  • The planetary k-index (kp-index) is an average of the measures taken in 15 locations in the world,
  • while the estimated kp-index is an average taken over only 8 locations. The data is available here:


If you take a look at their data for the last 30 days:


You can see that they give the k-index measured in two different locations (Fredericksburg and College), as well as the estimated kp-index. (You can ignore the columns marked "A".)

A single training point is the 8-value vector of the k-index in Frederickburg on a certain day, and its label is the 8-value vector of the estimated kp-index on that same day. Having trained the algorithm on a large number of sample points, a new query is the k-index in Frederickburg on a different day (with no label), and the algorithm must produce an 8-value vector which is its guess for the estimated kp-index on that day. For example, if the query vector was

2 2 1 2 1 2 1 2

(the k-index vector for Fredericksburg yesterday), then the answer

2.00 0.67 0.67 0.67 0.67 1.33 1.00 1.33

would have zero error, as this is exactly the estimated kp-index for that day. Of course, we can't really expect the algorithm to produce a solution with zero error.

My Intention:

An interesting question is can the k-index in either of the two individual locations be used to predict the estimated kp-index?

My Thoughts:

My thought is that this can be done by training a learning model (for example SVR) on some of the data (where a point is an 8-vector from either Fredericksburg or College representing some day, and its label is an 8-vector that is the estimated kp-index for the say day), and then testing its error on the rest of the data. I am programming in python, I might want to use the sklearn.multioutput.MultiOutputRegressor function.

My answer:

Given dataset.txt

#Importing necessary libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.svm import SVR
from sklearn.multioutput import MultiOutputRegressor
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
from sklearn.preprocessing import StandardScaler
import warnings

#Suppressing warnings for a cleaner 

# Load the dataset
dataset_path = "dataset.txt"
dataset = pd.read_csv(dataset_path, sep="\t", header=None)

# Split the space-separated data into columns
data = dataset[0].str.split(expand=True)

# Exclude unwanted columns: 4,5,6,7 and the last 5 columns
data = data.drop(columns=[3, 4, 5, 6] + list(range(data.shape[1]-5, data.shape[1])))

# Extracting date, kp-index, and Fredericksburg columns
date_columns = data.iloc[:, :3]
kp_index_columns = data.iloc[:, 3:11]
federickburg_columns = data.iloc[:, 11:19]

# Convert individual date columns to a single date column
date_columns['date'] = date_columns[0] + '-' + date_columns[1] + '-' + date_columns[2]

# Transforming Fredericksburg columns to float type
federickburg_transformed = federickburg_columns.astype(float)

# Defining the features (X) and target (y)
X = federickburg_transformed
y = kp_index_columns.astype(float)

# Splitting the data into training and testing sets (80% training, 20% testing)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Initialize the SVR model with a linear kernel
svr = SVR(kernel='linear', C=1)

# Using MultiOutputRegressor to handle multiple target columns
model = MultiOutputRegressor(svr).fit(X_train, y_train)

# Predict on the test data
y_pred = model.predict(X_test)

# Calculate the Mean Squared Error (MSE) and R^2 score for evaluation
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)

# Displaying the MSE and R^2 scores 
mse, r2 #(1.6405386761259917, 0.04202137156755004)

# Selecting a random sample from our features
random_index = np.random.randint(0, len(X))
sample_fredericksburg_random = X.iloc[random_index].values.reshape(1, -1)

# Predicting the kp-index for the randomly chosen Federicksburg k-index vector
predicted_kp_index_random = model.predict(sample_fredericksburg_random)

# Retrieving the actual kp-index for the chosen sample
actual_kp_index_random = y.iloc[random_index].values.reshape(1, -1)

# Scaling the features and target data for better performance during grid search
scaler_X = StandardScaler().fit(X_train)
scaler_y = StandardScaler().fit(y_train)
X_train_scaled = scaler_X.transform(X_train)
y_train_scaled = scaler_y.transform(y_train)
X_test_scaled = scaler_X.transform(X_test)
y_test_scaled = scaler_y.transform(y_test)

# Defining the hyperparameter grid for SVR
param_grid = {
 'estimator__C': [0.1, 1, 10, 100],
 'estimator__epsilon': [0.01, 0.1, 1],
 'estimator__kernel': ['linear', 'rbf']

# Performing grid search to find the best hyperparameters
grid_search = GridSearchCV(MultiOutputRegressor(SVR()), param_grid, cv=5, scoring='neg_mean_squared_error')
grid_search.fit(X_train_scaled, y_train_scaled)

# Using the best estimator found during grid search
best_svr = grid_search.best_estimator_

# Training the model with the best hyperparameters
best_svr.fit(X_train_scaled, y_train_scaled)

# Making predictions on the scaled test data
y_pred_scaled = best_svr.predict(X_test_scaled)

# Transforming the predicted data back to the original scale
y_pred_best = scaler_y.inverse_transform(y_pred_scaled)

# Calculating evaluation metrics for the optimized model
mse_best = mean_squared_error(y_test, y_pred_best)
r2_best = r2_score(y_test, y_pred_best)

# Displaying the improvements in MSE and R^2
mse_improvement = mse - mse_best
r2_improvement = r2_best - r2

# Displaying the results
mse_best, r2_best, mse_improvement, r2_improvement #(0.2357857371499288, 0.8438206360172018, 1.404752938976063, 0.8017992644496518)

# Predicting the kp-index for the previously chosen Federicksburg k-index vector using the optimized model
sample_fredericksburg_scaled = scaler_X.transform(sample_fredericksburg_random)
predicted_kp_index_scaled_optimized = best_svr.predict(sample_fredericksburg_scaled)

# Transforming the predicted data back to the original scale
predicted_kp_index_optimized = scaler_y.inverse_transform(predicted_kp_index_scaled_optimized)

# Displaying the optimized prediction and actual kp-index
predicted_kp_index_optimized, actual_kp_index_random 
#(array([[3.31721741, 2.61251221,1.63346665, 2.87117847,1.03796742,
#0.98798582, 1.65333422, 0.67988594]]),
#array([[3.333, 2.667, 1.667, 2.667, 1. , 1. , 1.667, 0.667]]))

# Extracting data for the last 30 days
date_last_30 = date_columns['date'].tail(30)
X_last_30 = X.tail(30)
y_last_30 = y.tail(30)

# Predicting KP-index for the last 30 days
y_pred_last_30_scaled = best_svr.predict(scaler_X.transform(X_last_30))
y_pred_last_30 = scaler_y.inverse_transform(y_pred_last_30_scaled)

# Visualizing the actual vs. predicted KP-index for the last 30 days
plt.figure(figsize=(20, 15))

 # Looping through each KP-index column
for i in range(8):
    plt.subplot(4, 2, i+1)
    #Plotting the actual values
    plt.plot(date_last_30, y_last_30.iloc[:, i].values, 'r', label='Actual')
    #Plotting the predicted values
    plt.plot(date_last_30, y_pred_last_30[:, i], 'b--', label='Predicted')
    plt.title(f'KP-Index Column {i+1}')
    if i == 0:

plt.suptitle('Actual vs Predicted KP-Index for Last 30 Days with Optimized Model', fontsize=16, y=1.05)

# Initializing the dataframe with the 'date' column
df_compact = pd.DataFrame()
df_compact['date'] = date_last_30.values

# Storing Fredericksburg values, predicted KP-index values, and actual KP-index values as lists
df_compact['federickburg'] = X_last_30.values.tolist()
df_compact['predicted_kp_index'] = y_pred_last_30.tolist()
df_compact['actual_kp_index'] = y_last_30.values.tolist()

# Rounding the values in 'predicted_kp_index' and 'actual_kp_index' columns to 2 decimal places
df_compact['predicted_kp_index'] = df_compact['predicted_kp_index'].apply(lambda x: [round(i, 2) for i in x])
df_compact['actual_kp_index'] = df_compact['actual_kp_index'].apply(lambda x: [round(i, 2) for i in x])

# Exporting the df_compact dataframe to a CSV file
csv_path = "pred_vs_act_last_30_days.csv"
df_compact.to_csv(csv_path, index=False)
df_compact.head(10) # Displaying the first few rows of the compact dataframe

# Calculating the evaluation metrics for the last 30 days' predictions
mae = mean_absolute_error(y_last_30, y_pred_last_30)
mse = mean_squared_error(y_last_30, y_pred_last_30)
r2 = r2_score(y_last_30, y_pred_last_30)
# Displaying the metrics
mae, mse, r2 #(0.06567479284845214, 0.010493205685711132, 0.9934506422188698)

My output images are:

enter image description here

enter image description here My question:

Is it possible to improve conciseness and complexity of my code? Or is there any alternative approach for this problem? Although I get the accuracy is almost 100% using the above algorithm.

  • \$\begingroup\$ You've made 10 revisions in a short amount of time. Do you consider the question finished now or are more incoming? \$\endgroup\$
    – Mast
    Sep 20 at 20:06
  • \$\begingroup\$ @Mast I don't change one word till now except spelling mistaken words. \$\endgroup\$
    – S. M.
    Sep 20 at 20:12


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.