deepSSF Training

Author

Affiliation

Scott Forrest

Queensland University of Technology, CSIRO

Abstract

In this script, we will train a deepSSF model on the training data. In this case the training data was generated using the deepSSF_data_prep_id.qmd script, which crops out local images for each step of the observed telemetry data.

Model overview

Refer to the paper and the Model Overview section for a conceptual overview of the model.

Detect computing environment

If using Google Colab, mount the drive and set the base directory to the working folder. If using local, set the base directory.

import os       # Operating system utilities
import sys

# Detect environment
def is_colab():
    """Returns True if running in Google Colab, False otherwise."""
    try:
        import google.colab
        return True
    except ImportError:
        return False

# Set up environment-specific configurations
if is_colab():

    # Colab-specific setup
    !pip install rasterio
    from google.colab import drive
    drive.mount('/content/drive')
    sys.path.append('/content/drive/MyDrive/GitHub/deepSSF/Python')

    # for saving plots etc
    base_path = '/content/drive/MyDrive/GitHub/deepSSF'
    print("Running in Google Colab environment")

else:

    # Local environment setup
    base_path = '..'
    print("Running in local environment")

# Now you can use base_path regardless of environment
print(f"Using base path: {base_path}")

Running in local environment
Using base path: ..

Import packages

print(sys.version)  # Print Python version in use

import numpy as np                                      # Array operations
import matplotlib.pyplot as plt                         # Plotting library
import mpmath as mp                                     # Math library
import torch                                            # Main PyTorch library
import torch.optim as optim                             # Optimization algorithms
import torch.nn as nn                                   # Neural network modules
import os                                               # Operating system utilities
import glob                                             # Pattern matching
import imageio.v2 as imageio                            # Image manipulation - for creating GIFs
from IPython.display import Image, display              # For plotting GIFs
import pandas as pd                                     # Data manipulation
import rasterio                                         # Geospatial raster data

from torch.utils.data import Dataset, DataLoader        # Dataset and batch data loading
from datetime import datetime, timedelta                # Date/time utilities
from rasterio.plot import show                           # Plot raster data

import deepSSF_utils                                    # deepSSF utility functions

# Get today's date
today_date = datetime.today().strftime('%Y-%m-%d')
print(today_date)

# Set random seed for reproducibility
# seed = 42

3.12.5 | packaged by Anaconda, Inc. | (main, Sep 12 2024, 18:18:29) [MSC v.1929 64 bit (AMD64)]
2025-07-09

Set the device (accelerator - cuda for NVIDIA GPU or mps for Mac)

# Set the device to be used (GPU or CPU)
if torch.cuda.is_available():
    device = "cuda"
elif torch.backends.mps.is_available():  # For Mac M1/M2/M3
    device = "mps"
else:
    device = "cpu"
    
print(f"Using {device} device")

if torch.backends.mps.is_available():
    # Set default tensor type for PyTorch
    torch.set_default_dtype(torch.float32)
    print('Set default tensor type to float32')

Using cpu device

Select individual and create directory to save model weights and outputs

# select the id to train the model on
buffalo_id = 2005
# in our case the actual dataset will be slightly smaller due to steps being removed that were outside the extent
n_samples = 10297
# n_samples = 100

Create a directory to save the outputs

If we have already run this code today, we will add update index to create a new folder

# Count existing directories with similar pattern
pattern = f'{base_path}/Python/outputs/model_training/id{buffalo_id}_deepSSF_training_*_{today_date}'
existing_dirs = glob.glob(pattern)
dir_index = len(existing_dirs) + 1

# Create directory with index
output_dir = f'{base_path}/Python/outputs/model_training/id{buffalo_id}_deepSSF_training_{dir_index}_{today_date}'
os.makedirs(output_dir, exist_ok=True)

print(f"Created directory: {output_dir}")

# To use an existing directory for loading trained model
# output_dir = f'{base_path}/Python/outputs/model_training/id2005_2025-04-01'

Created directory: ../Python/outputs/model_training/id2005_deepSSF_training_1_2025-07-09

Set paths to data

# Specify the path to CSV file
# csv_file_path = f'{base_path}/buffalo_local_data_id/buffalo_{buffalo_id}_data_df_lag_1hr_n{n_samples}.csv'
csv_file_path = f'{base_path}/buffalo_local_data_id/buffalo_temporal_cont_{buffalo_id}_data_df_lag_1hr_n{n_samples}.csv'

# Path to your TIF file
ndvi_path = f'{base_path}/buffalo_local_layers_id/buffalo_{buffalo_id}_ndvi_cent101x101_lag_1hr_n{n_samples}.tif'
# Path to your TIF file
canopy_path = f'{base_path}/buffalo_local_layers_id/buffalo_{buffalo_id}_canopy_cent101x101_lag_1hr_n{n_samples}.tif'
# Path to your TIF file
herby_path = f'{base_path}/buffalo_local_layers_id/buffalo_{buffalo_id}_herby_cent101x101_lag_1hr_n{n_samples}.tif'
# Path to your TIF file
slope_path = f'{base_path}/buffalo_local_layers_id/buffalo_{buffalo_id}_slope_cent101x101_lag_1hr_n{n_samples}.tif'

# Path to your TIF file
# pres_path = f'{base_path}/buffalo_local_layers_id/buffalo_{buffalo_id}_pres_cent101x101_lag_1hr_n{n_samples}.tif'
pres_path = f'{base_path}/buffalo_local_layers_id/fixed_buffalo_{buffalo_id}_pres_cent101x101_lag_1hr_n{n_samples}.tif'
# pres_path = f'{base_path}/buffalo_local_layers_id/within_cell_buffalo_{buffalo_id}_pres_cent101x101_lag_1hr_n{n_samples}.tif'

Read buffalo data

# Read the CSV file into a DataFrame
buffalo_df = pd.read_csv(csv_file_path)
print(buffalo_df.shape)

# Lag the values in column 'A' by one index to get the bearing of the previous step
buffalo_df['bearing_tm1'] = buffalo_df['bearing'].shift(1)
# Pad the missing value with a specified value, e.g., 0
buffalo_df['bearing_tm1'] = buffalo_df['bearing_tm1'].fillna(0)

# Display the first few rows of the DataFrame
print(buffalo_df.head())

(10103, 43)
             x_            y_                    t_    id           x1_  \
0  41969.310875 -1.435671e+06  2018-07-25T01:04:23Z  2005  41969.310875   
1  41921.521939 -1.435654e+06  2018-07-25T02:04:39Z  2005  41921.521939   
2  41779.439594 -1.435601e+06  2018-07-25T03:04:17Z  2005  41779.439594   
3  41841.203272 -1.435635e+06  2018-07-25T04:04:39Z  2005  41841.203272   
4  41655.463332 -1.435604e+06  2018-07-25T05:04:27Z  2005  41655.463332   

            y1_           x2_           y2_     x2_cent    y2_cent  ...  \
0 -1.435671e+06  41921.521939 -1.435654e+06  -47.788936  16.857110  ...   
1 -1.435654e+06  41779.439594 -1.435601e+06 -142.082345  53.568427  ...   
2 -1.435601e+06  41841.203272 -1.435635e+06   61.763677 -34.322938  ...   
3 -1.435635e+06  41655.463332 -1.435604e+06 -185.739939  31.003534  ...   
4 -1.435604e+06  41618.651923 -1.435608e+06  -36.811409  -4.438037  ...   

    bearing  bearing_sin  bearing_cos        ta    cos_ta         x_min  \
0  2.802478     0.332652    -0.943050  1.367942  0.201466  40706.810875   
1  2.781049     0.352783    -0.935705 -0.021429  0.999770  40659.021939   
2 -0.507220    -0.485749     0.874098  2.994917 -0.989262  40516.939594   
3  2.976198     0.164641    -0.986354 -2.799767 -0.942144  40578.703272   
4 -3.021610    -0.119695    -0.992811  0.285377  0.959556  40392.963332   

          x_max         y_min         y_max  bearing_tm1  
0  43231.810875 -1.436934e+06 -1.434409e+06     0.000000  
1  43184.021939 -1.436917e+06 -1.434392e+06     2.802478  
2  43041.939594 -1.436863e+06 -1.434338e+06     2.781049  
3  43103.703272 -1.436898e+06 -1.434373e+06    -0.507220  
4  42917.963332 -1.436867e+06 -1.434342e+06     2.976198  

[5 rows x 44 columns]

Spatial data

NDVI

# Using rasterio
with rasterio.open(ndvi_path) as ndvi:
    # Read all layers/channels into a single numpy array
    # rasterio indexes channels starting from 1, hence the range is 1 to src.count + 1
    ndvi_stack = ndvi.read([i for i in range(1, ndvi.count + 1)])

print(ndvi_stack.shape)

(10103, 101, 101)

Normalise the layers

# Replace NaNs in the original array with -1, which represents water
ndvi_stack = np.nan_to_num(ndvi_stack, nan=-1.0)

# Convert the numpy array to a PyTorch tensor, which is the format required for training the model
ndvi_tens = torch.from_numpy(ndvi_stack)
print(ndvi_tens.dtype)

# Print the shape of the PyTorch tensor
print(ndvi_tens.shape)

# Print the mean, max, and min values of the NDVI tensor
ndvi_mean = torch.mean(ndvi_tens)
ndvi_max = torch.max(ndvi_tens)
ndvi_min = torch.min(ndvi_tens)
print("Mean = ", ndvi_mean)
print("Max = ", ndvi_max)
print("Min = ", ndvi_min)

# Normalizing the data
ndvi_tens_norm = (ndvi_tens - ndvi_min) / (ndvi_max - ndvi_min)
print("Mean = ", torch.mean(ndvi_tens_norm))
print("Max = ", torch.max(ndvi_tens_norm))
print("Min = ", torch.min(ndvi_tens_norm))

torch.float32
torch.Size([10103, 101, 101])
Mean =  tensor(0.3048)
Max =  tensor(0.8220)
Min =  tensor(-0.2894)
Mean =  tensor(0.5347)
Max =  tensor(1.)
Min =  tensor(0.)

Plot a single NDVI layer

for i in range(0, 1):
    plt.imshow(ndvi_tens_norm[i].numpy())
    plt.colorbar()
    plt.show()

Canopy cover

# Using rasterio
with rasterio.open(canopy_path) as canopy:
    # Read all layers/channels into a single numpy array
    # rasterio indexes channels starting from 1, hence the range is 1 to src.count + 1
    canopy_stack = canopy.read([i for i in range(1, canopy.count + 1)])

print(canopy_stack.shape)

(10103, 101, 101)

# Convert the numpy array to a PyTorch tensor, which is the format required for training the model
canopy_tens = torch.from_numpy(canopy_stack)
print(canopy_tens.shape)

# Print the mean, max, and min values of the canopy tensor
print("Mean = ", torch.mean(canopy_tens))
canopy_max = torch.max(canopy_tens)
canopy_min = torch.min(canopy_tens)
print("Max = ", canopy_max)
print("Min = ", canopy_min)

# Normalizing the data
canopy_tens_norm = (canopy_tens - canopy_min) / (canopy_max - canopy_min)
print("Mean = ", torch.mean(canopy_tens_norm))
print("Max = ", torch.max(canopy_tens_norm))
print("Min = ", torch.min(canopy_tens_norm))

torch.Size([10103, 101, 101])
Mean =  tensor(44.3548)
Max =  tensor(82.5000)
Min =  tensor(0.)
Mean =  tensor(0.5376)
Max =  tensor(1.)
Min =  tensor(0.)

for i in range(0, 1):
    plt.imshow(canopy_tens_norm[i].numpy())
    plt.colorbar()
    plt.show()

Herbaceous vegetation

# Using rasterio
with rasterio.open(herby_path) as herby:
    # Read all layers/channels into a single numpy array
    # rasterio indexes channels starting from 1, hence the range is 1 to src.count + 1
    herby_stack = herby.read([i for i in range(1, herby.count + 1)])

print(herby_stack.shape)

(10103, 101, 101)

# Convert the numpy array to a PyTorch tensor, which is the format required for training the model
herby_tens = torch.from_numpy(herby_stack)
print(herby_tens.shape)

# Print the mean, max, and min values of the herby tensor
print("Mean = ", torch.mean(herby_tens))
herby_max = torch.max(herby_tens)
herby_min = torch.min(herby_tens)
print("Max = ", herby_max)
print("Min = ", herby_min)

# Normalizing the data
herby_tens_norm = (herby_tens - herby_min) / (herby_max - herby_min)
print("Mean = ", torch.mean(herby_tens_norm))
print("Max = ", torch.max(herby_tens_norm))
print("Min = ", torch.min(herby_tens_norm))

torch.Size([10103, 101, 101])
Mean =  tensor(0.8069)
Max =  tensor(1.)
Min =  tensor(0.)
Mean =  tensor(0.8069)
Max =  tensor(1.)
Min =  tensor(0.)

for i in range(0, 1):
    plt.imshow(herby_tens_norm[i])
    plt.colorbar()
    plt.show()

Slope

# Using rasterio
with rasterio.open(slope_path) as slope:
    # Read all layers/channels into a single numpy array
    # rasterio indexes channels starting from 1, hence the range is 1 to src.count + 1
    slope_stack = slope.read([i for i in range(1, slope.count + 1)])

print(slope_stack.shape)

(10103, 101, 101)

# Convert the numpy array to a PyTorch tensor, which is the format required for training the model
slope_tens = torch.from_numpy(slope_stack)
print(slope_tens.shape)

# Print the mean, max, and min values of the slope tensor
print("Mean = ", torch.mean(slope_tens))
slope_max = torch.max(slope_tens)
slope_min = torch.min(slope_tens)
print("Max = ", slope_max)
print("Min = ", slope_min)

# Normalizing the data
slope_tens_norm = (slope_tens - slope_min) / (slope_max - slope_min)
print("Mean = ", torch.mean(slope_tens_norm))
print("Max = ", torch.max(slope_tens_norm))
print("Min = ", torch.min(slope_tens_norm))

torch.Size([10103, 101, 101])
Mean =  tensor(0.7779)
Max =  tensor(12.2981)
Min =  tensor(0.0006)
Mean =  tensor(0.0632)
Max =  tensor(1.)
Min =  tensor(0.)

for i in range(0, 1):
    plt.imshow(slope_tens_norm[i])
    plt.colorbar()
    plt.show()

Presence records - target of model

This is what the model is trying to predict, which is the location of the next step.

# Using rasterio
with rasterio.open(pres_path) as pres:
    # Read all layers/channels into a single numpy array
    # rasterio indexes channels starting from 1, hence the range is 1 to src.count + 1
    pres_stack = pres.read([i for i in range(1, pres.count + 1)])

print(pres_stack.shape)

(10103, 101, 101)

for i in range(0, 1):
    plt.imshow(pres_stack[i])
    plt.show()

Combine the spatial layers into channels

# Stack the channels along a new axis; here, 1 is commonly used for the channel axis in PyTorch
combined_stack = torch.stack([ndvi_tens_norm,
                              canopy_tens_norm,
                              herby_tens_norm,
                              slope_tens_norm],
                              dim=1)

print(combined_stack.shape)

torch.Size([10103, 4, 101, 101])

From the size we can see that there are 10103 samples (steps), 4 channels (layers = covariates) and 101x101 pixels.

Defining data sets and data loaders

Creating a dataset class

This custom PyTorch Dataset organizes all your input (spatial data, scalar covariates, bearing, and target) in a single object, allowing you to neatly manage how samples are accessed. The __init__ method prepares and stores all the data, __len__ returns the total number of samples, and __getitem__ retrieves a single sample by index—enabling straightforward batching and iteration when used with a DataLoader.

class buffalo_data(Dataset):

    def __init__(self):
        # data loading. Here we are just using the combined_stack as the spatial covariates
        self.spatial_data_x = combined_stack

        # the scalar data that will be converted to spatial data and added as channels to the spatial covariates
        self.scalar_to_grid_data = torch.from_numpy(buffalo_df[['hour_t2_sin',
                                                                'hour_t2_cos',
                                                                'yday_t2_sin',
                                                                'yday_t2_cos']].values).float()

        # the bearing data that will be added as a channel to the spatial covariates
        self.bearing_x = torch.from_numpy(buffalo_df[['bearing_tm1']].values).float()

        # the target data
        self.target = torch.tensor(pres_stack)

        # number of samples
        self.n_samples = self.spatial_data_x.shape[0]

    def __len__(self):
        # allows for the use of len() function
        return self.n_samples

    def __getitem__(self, index):
        # allows for indexing of the dataset
        return self.spatial_data_x[index], self.scalar_to_grid_data[index], self.bearing_x[index], self.target[index]

Now we can create an instance of the dataset class and check that is working as expected.

# Create an instance of our custom buffalo_data Dataset:
dataset = buffalo_data()

# Print the total number of samples loaded (determined by n_samples in the dataset):
n_samples_loaded = dataset.n_samples
print(n_samples_loaded)

# Retrieve *all* samples (using the slice dataset[:] invokes __getitem__ on all indices).
# This returns a tuple of (spatial data, scalar-to-grid data, bearing data, target labels).
features1, features2, features3, labels = dataset[:]

# Examine the dimensions of each returned tensor for verification:

# Spatial data
print(features1.shape)

# Scalar-to-grid data
print(features2.shape)

# Bearing data
print(features3.shape)

# Target labels
print(labels.shape)

10103
torch.Size([10103, 4, 101, 101])
torch.Size([10103, 4])
torch.Size([10103, 1])
torch.Size([10103, 101, 101])

Split into training, validation and test sets

Uncomment either the random split or the sequential split to use the one you prefer. A sequential split requires more of an extrapolation, as the individual could have encountered a different environmental region in the latter part of the time series, or it could be from a different time period, such as part of a season that isn’t in the training data.

For testing we will therefore use the sequential split as more of a test for the model.

training_split = 0.8 # 80% of the data will be used for training
validation_split = 0.1 # 10% of the data will be used for validation (deciding when to stop training)
test_split = 0.1 # 10% of the data will be used for testing (model evaluation)

To split the samples randomly we can use the random_split function from PyTorch

# dataset_train, dataset_val, dataset_test = torch.utils.data.random_split(dataset, [training_split, validation_split, test_split])

To split the data sequentially, we can use the Subset function from PyTorch, which allows us to select a subset of the data based on the indices we provide.

# For sequential split (need integers)
n_samples = len(dataset)
n_train = int(training_split * n_samples)
n_val = int(validation_split * n_samples)
n_test = n_samples - n_train - n_val  # Ensure they sum to total number of samples

# Get the start and end indices for each split
train_end = int(training_split * n_samples)
val_end = int((training_split + validation_split) * n_samples)

train_indices = list(range(0, train_end))
val_indices = list(range(train_end, val_end))
test_indices = list(range(val_end, len(dataset)))

# to split the samples sequentially
dataset_train = torch.utils.data.Subset(dataset, train_indices)
dataset_val = torch.utils.data.Subset(dataset, val_indices)
dataset_test = torch.utils.data.Subset(dataset, test_indices)

print("Number of training samples: ",   len(dataset_train))
print("Number of validation samples: ", len(dataset_val))
print("Number of testing samples: ",    len(dataset_test))

Number of training samples:  8082
Number of validation samples:  1010
Number of testing samples:  1011

Create dataloaders

The DataLoader in PyTorch wraps an iterable around the Dataset to enable easy access to the samples.

# Define the batch size for how many samples to process at once in each step:
batch_size = 32

# Create a DataLoader for the training dataset with a batch size of batch_size, and shuffle samples
# so that the model doesn't see data in the same order each epoch.
dataloader_train = DataLoader(dataset=dataset_train,
                              batch_size=batch_size,
                              shuffle=True)

# Create a DataLoader for the validation dataset, also with a batch size of batch_size and shuffling.
# Even though it's not always mandatory to shuffle validation data, some users keep the same setting.
dataloader_val = DataLoader(dataset=dataset_val,
                            batch_size=batch_size,
                            shuffle=True)

# Create a DataLoader for the test dataset, likewise with a batch size of batch_size and shuffling.
# As we want to index the testing data for plotting, we will not shuffle the test data.
dataloader_test = DataLoader(dataset=dataset_test,
                             batch_size=batch_size,
                             shuffle=False)

Check that the data loader is working as expected.

# Display image and label.
# next(iter(dataloader_train)) returns the next batch of the training data
features1, features2, features3, labels = next(iter(dataloader_train))
print(f"Feature 1 batch shape:  {features1.size()}")
print(f"Feature 2 batch shape:  {features2.size()}")
print(f"Feature 3 batch shape:  {features3.size()}")
print(f"Labels batch shape:     {labels.size()}")

Feature 1 batch shape:  torch.Size([32, 4, 101, 101])
Feature 2 batch shape:  torch.Size([32, 4])
Feature 3 batch shape:  torch.Size([32, 1])
Labels batch shape:     torch.Size([32, 101, 101])

Define the model

Deep learning can be considered as a sequence of blocks, each of which perform some (typically nonlinear) transformation on input data to produce some output. Providing each block has the appropriate inputs, they can be combined to build a larger network that is capable of achieving complex and abstract transformations and can be used to represent complex processes.

A block is modular component of a neural network, in our case defined as a Python class (type of object with certain functionality described by its definition) inheriting from torch.nn.Module in PyTorch. A block encapsulates a sequence of operations, including layers (such as fully connected layers or convolutional layers) and activation functions, to process input data. Each block has a forward method (i.e. instructions) that defines the data flow through the network during inference or training.

Convolutional block for the habitat selection subnetwork

This block is a convolutional layer that takes in the spatial covariates (including the layers created from the scalar values such as time), goes through a series of convolution operations and ReLU activation functions and outputs a feature map, which is the habitat selection probability surface.

class Conv2d_block_spatial(nn.Module):
    def __init__(self, params):
        super(Conv2d_block_spatial, self).__init__()

        # define the parameters
        self.batch_size = params.batch_size
        self.input_channels = params.input_channels
        self.output_channels = params.output_channels
        self.kernel_size = params.kernel_size
        self.stride = params.stride
        self.padding = params.padding
        self.image_dim = params.image_dim
        self.device = params.device

        # define the layers - nn.Sequential allows for the definition of layers in a sequential manner
        self.conv2d = nn.Sequential(

        # convolutional layer 1
        nn.Conv2d(in_channels=self.input_channels,
                  out_channels=self.output_channels,
                  kernel_size=self.kernel_size,
                  stride=self.stride,
                  padding=self.padding),
        # ReLU activation function
        nn.ReLU(),

        # convolutional layer 2
        nn.Conv2d(in_channels=self.output_channels,
                  out_channels=self.output_channels,
                  kernel_size=self.kernel_size,
                  stride=self.stride,
                  padding=self.padding),
        # ReLU activation function
        nn.ReLU(),

        # convolutional layer 3, which outputs a single layer, which is the habitat selection map
        nn.Conv2d(in_channels=self.output_channels,
                  out_channels=1,
                  kernel_size=self.kernel_size,
                  stride=self.stride,
                  padding=self.padding)
        )

    # define the forward pass of the model, i.e. how the data flows through the model
    def forward(self, x):

        # self.conv2d(x) passes the input through the convolutional layers, and the squeeze function removes the channel dimension, resulting in a 2D tensor (habitat selection map)
        # print("Shape before squeeze:", self.conv2d(x).shape) # Debugging print
        conv2d_spatial = self.conv2d(x).squeeze(dim = 1)

        # normalise to sum to 1
        # print("Shape before logsumexp:", conv2d_spatial.shape) # Debugging print
        conv2d_spatial = conv2d_spatial - torch.logsumexp(conv2d_spatial, dim = (1, 2), keepdim = True)

        # output the habitat selection map
        return conv2d_spatial

Convolutional block for the movement subnetwork

This block is also convolutional layer, with the same inputs, but this block also has max pooling layers to reduce the spatial resolution of the feature maps whilst preserving the most prominent features in the feature maps, and outputs a ‘flattened’ feature map. A flattened feature map is a 1D tensor (a vector) that can be used as input to a fully connected layer.

class Conv2d_block_toFC(nn.Module):
    def __init__(self, params):
        super(Conv2d_block_toFC, self).__init__()

        # define the parameters
        self.batch_size = params.batch_size
        self.input_channels = params.input_channels
        self.output_channels_movement = params.output_channels_movement
        self.kernel_size = params.kernel_size
        self.stride = params.stride
        self.kernel_size_mp = params.kernel_size_mp
        self.stride_mp = params.stride_mp
        self.padding = params.padding
        self.image_dim = params.image_dim
        self.device = params.device

        # define the layers - nn.Sequential allows for the definition of layers in a sequential manner
        self.conv2d = nn.Sequential(

        # convolutional layer 1
        nn.Conv2d(in_channels=self.input_channels,
                  out_channels=self.output_channels_movement,
                  kernel_size=self.kernel_size,
                  stride=self.stride,
                  padding=self.padding),
        # ReLU activation function
        nn.ReLU(),

        # max pooling layer 1 (reduces the spatial dimensions of the data whilst retaining the most important features)
        nn.MaxPool2d(kernel_size=self.kernel_size_mp,
                     stride=self.stride_mp),

        # convolutional layer 2
        nn.Conv2d(in_channels=self.output_channels_movement,
                  out_channels=self.output_channels_movement,
                  kernel_size=self.kernel_size,
                  stride=self.stride,
                  padding=self.padding),
        # ReLU activation function
        nn.ReLU(),

        # max pooling layer 2
        nn.MaxPool2d(kernel_size=self.kernel_size_mp,
                     stride=self.stride_mp),

        # # to add a third convolutional layer, uncomment the following lines
        # # convolutional layer 3
        # nn.Conv2d(in_channels=self.output_channels_movement,
        #           out_channels=self.output_channels_movement,
        #           kernel_size=self.kernel_size,
        #           stride=self.stride,
        #           padding=self.padding),
        # # ReLU activation function
        # nn.ReLU(),

        # # max pooling layer 3
        # nn.MaxPool2d(kernel_size=self.kernel_size_mp,
        #              stride=self.stride_mp),

        # flatten the data to pass through the fully connected layer
        nn.Flatten())

    def forward(self, x):

        # self.conv2d(x) passes the input through the convolutional layers, and outputs a 1D tensor
        return self.conv2d(x)

Fully connected block for the movement subnetwork

This block takes in the flattened feature map from the previous block, passes through several fully connected layers, which extracts information from the spatial covariates that is relevant for movement, and outputs the parameters that define the movement kernel.

class FCN_block_all_movement(nn.Module):
    def __init__(self, params):
        super(FCN_block_all_movement, self).__init__()

        # define the parameters
        self.batch_size = params.batch_size
        self.dense_dim_in_all = params.dense_dim_in_all
        self.dense_dim_hidden = params.dense_dim_hidden
        self.image_dim = params.image_dim
        self.device = params.device
        self.num_movement_params = params.num_movement_params
        self.dropout = params.dropout

        # define the layers - nn.Sequential allows for the definition of layers in a sequential manner
        self.ffn = nn.Sequential(

            # fully connected layer 1 (the dense_dim_in_all is the number of input features,
            # and should match the output of the Conv2d_block_toFC block).
            # the dense_dim_hidden is the number of neurons in the hidden layer, and doesn't need to be the same as the input features
            nn.Linear(self.dense_dim_in_all, self.dense_dim_hidden),
            # dropout layer (helps to reduce overfitting)
            nn.Dropout(self.dropout),
            # ReLU activation function
            nn.ReLU(),

            # fully connected layer 2
            # the number of input neurons should match the output from the previous layer
            nn.Linear(self.dense_dim_hidden, self.dense_dim_hidden),
            # dropout layer
            nn.Dropout(self.dropout),
            # ReLU activation function
            nn.ReLU(),

            # fully connected layer 3
            # the number of input neurons should match the output from the previous layer,
            # and the number of output neurons should match the number of movement parameters
            nn.Linear(self.dense_dim_hidden, self.num_movement_params)

        )

    def forward(self, x):

        # self.ffn(x) passes the input through the fully connected layers, and outputs a 1D tensor (vector of movement parameters)
        return self.ffn(x)

Block to convert the movement parameters to a probability distribution

What the block does

This block is a bit longer and more involved, but there are no parameters in here that need to be learned (estimated). It is just a series of operations that are applied to the movement parameters to convert them to a probability distribution.

This block takes in the movement parameters and converts them to a probability distribution. This essentially just applies the appropriate density functions using the parameter values predicted by the movement blocks, which in our case is a finite mixture of Gamma distributions and a finite mixture of von Mises distributions.

The formulation of predicting parameters and converting them to a movement kernel ensures that the movement kernel is very flexible, and can be any combination of distributions, which need not all be the same (e.g., a step length distribution may be combination of a Gamma and a log-normal distribution).

Constraints

One constraint to ensure that we can perform backpropagation is that the entire forward pass, including the block below that produces the density functions, must be differentiable with respect to the parameters of the model. PyTorch’s torch.distributions module and its special functions (e.g., torch.special) provide differentiable implementations for many common distributions. Examples are the

• Gamma function for the (log) Gamma distribution, torch.lgamma()
• The modified Bessel function of the first kind of order 0 for the von Mises distribution, torch.special.i0()

Some of the movement parameters, such as the shape and scale of the Gamma distribution, must be positive. We therefore exponentiate them in this block to ensure that they are positive. This means that the model is actually learning the log of the shape and scale parameters. For the von Mises mu parameters however, they can be any value, so we do not need to exponentiate them. We could constrain them to be between -pi and pi, but this is not necessary as the von Mises distribution is periodic, so any value will be equivalent to another value that is within the range -pi to pi.

Notes

To help with identifiability, it is possible to fix certain parameter values, such as the mu parameters in the mixture of von Mises distributions to pi and -pi for instance (one would then reduce the number of predicted parameters by the previous block, as these no longer need to be predicted).

We can also transform certain parameters such that they are being estimated in a similar range (analagous to standardising variables in linear regression). In our case we know that the scale parameter of one of the Gamma distributions is around 500. What we can then do after exponentiating is multiply the scale parameter by 500, so the model is learning the log of the scale parameter divided by 500. This will ensure that this parameter is in a similar range to the other parameters, and can help with convergence. To do this we:

Pull out the relevant parameters from the input tensor (output of previous block) - gamma_scale2 = torch.exp(x[:, 4]).unsqueeze(0).unsqueeze(0)

Multiply the scale parameter by 500, so the model is learning the log of the scale parameter divided by 500 - gamma_scale2 = gamma_scale2 * 500

Consideration for the centre cell

As the centre cell has a distance of exactly 0 (using the distance layer we created), this can cause numerical issues for the gamma distribution, as when the shape parameter is less than one the mode approaches infinity. To avoid this, we can add a small value to the central cell of the distance layer, so that the distance is never exactly 0.

To get the value to use in the central cell, we will calculate the average distance from the very centre to any point in the cell (assuming that the distance within the cell is continuous). This comes out to be:

\(\int_{-0.5}^{0.5} \int_{-0.5}^{0.5} \sqrt{x^2 + y^2} \, dx \, dy\)

We calculate a constant numerically below:

def integrand(x, y):
    return mp.sqrt(x**2 + y**2)

val = mp.quad(lambda Y:
              mp.quad(lambda X: integrand(X, Y),
                      [-0.5, 0.5]),
                      [-0.5, 0.5])

print(val)

0.382597668656132

class Params_to_Grid_Block(nn.Module):
    def __init__(self, params):
        super(Params_to_Grid_Block, self).__init__()

        # define the parameters
        self.batch_size = params.batch_size
        self.image_dim = params.image_dim
        self.pixel_size = params.pixel_size

        # create distance and bearing layers
        # determine the distance of each pixel from the centre of the image
        self.center = self.image_dim // 2
        y, x = np.indices((self.image_dim, self.image_dim))
        self.distance_layer = torch.from_numpy(np.sqrt((self.pixel_size*(x - self.center))**2 +
                                                       (self.pixel_size*(y - self.center))**2)).float()
        # change the centre cell to the average distance from the centre to the edge of the pixel

        # average distance from the centre to any point within the pixel
        # calculated as a double integral of sqrt(x^2 + y^2) dx dy over the area of the pixel
        self.distance_layer[self.center, self.center] = 0.3826*self.pixel_size

        # determine the bearing of each pixel from the centre of the image
        self.bearing_layer = torch.from_numpy(np.arctan2(self.center - y,
                                                         x - self.center)).float()
        self.device = params.device


    # Gamma densities (on the log-scale) for the mixture distribution
    def gamma_density(self, x, shape, scale):
        # Ensure all tensors are on the same device as x
        shape = shape.to(x.device)
        scale = scale.to(x.device)
        # return -1*torch.lgamma(shape) -shape*torch.log(scale) + (shape - 1)*torch.log(x) - x/scale

        # to account for change of variables
        return (-1*torch.lgamma(shape) -shape*torch.log(scale) + (shape - 1)*torch.log(x) - x/scale) - torch.log(x)

    # log von Mises densities (on the log-scale) for the mixture distribution
    def vonmises_density(self, x, kappa, vm_mu):
        # Ensure all tensors are on the same device as x
        kappa = kappa.to(x.device)
        vm_mu = vm_mu.to(x.device)
        return kappa*torch.cos(x - vm_mu) - 1*(np.log(2*torch.pi) + torch.log(torch.special.i0(kappa)))


    def forward(self, x, bearing):

        # parameters of the first mixture distribution
        # x are the outputs from the fully connected layers (vector of movement parameters)
        # we therefore need to extract the appropriate parameters
        # the locations are not specific to any specific parameters, as long as any aren't extracted more than once

        # Gamma distributions

        # pull out the parameters of the first gamma distribution and exponentiate them to ensure they are positive
        # the unsqueeze function adds a new dimension to the tensor
        # we do this twice to match the dimensions of the distance_layer,
        # and then repeat the parameter value across a grid, such that the density can be calculated at every cell/pixel
        gamma_shape1 = torch.exp(x[:, 0]).unsqueeze(0).unsqueeze(0)
        gamma_shape1 = gamma_shape1.repeat(self.image_dim, self.image_dim, 1)
        # this just changes the order of the dimensions to match the distance_layer
        gamma_shape1 = gamma_shape1.permute(2, 0, 1)

        gamma_scale1 = torch.exp(x[:, 1]).unsqueeze(0).unsqueeze(0)
        gamma_scale1 = gamma_scale1.repeat(self.image_dim, self.image_dim, 1)
        gamma_scale1 = gamma_scale1.permute(2, 0, 1)

        # gamma_weight1 = torch.exp(x[:, 2]).unsqueeze(0).unsqueeze(0)
        gamma_weight1 = x[:, 2].unsqueeze(0).unsqueeze(0)
        gamma_weight1 = gamma_weight1.repeat(self.image_dim, self.image_dim, 1)
        gamma_weight1 = gamma_weight1.permute(2, 0, 1)

        # parameters of the second mixture distribution
        gamma_shape2 = torch.exp(x[:, 3]).unsqueeze(0).unsqueeze(0)
        gamma_shape2 = gamma_shape2.repeat(self.image_dim, self.image_dim, 1)
        gamma_shape2 = gamma_shape2.permute(2, 0, 1)

        gamma_scale2 = torch.exp(x[:, 4]).unsqueeze(0).unsqueeze(0)
        gamma_scale2 = gamma_scale2 * 500 ### transform the scale parameter so it can be estimated near the same range as the other parameters
        gamma_scale2 = gamma_scale2.repeat(self.image_dim, self.image_dim, 1)
        gamma_scale2 = gamma_scale2.permute(2, 0, 1)

        # gamma_weight2 = torch.exp(x[:, 5]).unsqueeze(0).unsqueeze(0)
        gamma_weight2 = x[:, 5].unsqueeze(0).unsqueeze(0)
        gamma_weight2 = gamma_weight2.repeat(self.image_dim, self.image_dim, 1)
        gamma_weight2 = gamma_weight2.permute(2, 0, 1)

        # Apply softmax to the mixture weights to ensure they sum to 1
        gamma_weights = torch.stack([gamma_weight1, gamma_weight2], dim=0)
        gamma_weights = torch.nn.functional.softmax(gamma_weights, dim=0)
        gamma_weight1 = gamma_weights[0]
        gamma_weight2 = gamma_weights[1]

        # calculation of Gamma densities
        gamma_density_layer1 = self.gamma_density(self.distance_layer,
                                                  gamma_shape1,
                                                  gamma_scale1).to(device)

        gamma_density_layer2 = self.gamma_density(self.distance_layer,
                                                  gamma_shape2,
                                                  gamma_scale2).to(device)

        # combining both densities to create a mixture distribution using logsumexp
        logsumexp_gamma_corr = torch.max(gamma_density_layer1, gamma_density_layer2)
        gamma_density_layer = logsumexp_gamma_corr + torch.log(gamma_weight1 * torch.exp(gamma_density_layer1 - logsumexp_gamma_corr) +
                                                               gamma_weight2 * torch.exp(gamma_density_layer2 - logsumexp_gamma_corr))
        # print(torch.sum(gamma_density_layer))
        # print(torch.sum(torch.exp(gamma_density_layer)))


        ## Von Mises Distributions

        # calculate the new bearing from the turning angle
        # takes in the bearing from the previous step and adds the turning angle, which is estimated by the model
        # we do not exponentiate the von Mises mu parameters as we want to allow them to be negative
        bearing_new1 = x[:, 6] + bearing[:, 0]

        # the new bearing becomes the mean of the von Mises distribution
        vonmises_mu1 = bearing_new1.unsqueeze(0).unsqueeze(0)
        vonmises_mu1 = vonmises_mu1.repeat(self.image_dim, self.image_dim, 1)
        vonmises_mu1 = vonmises_mu1.permute(2, 0, 1)

        # parameters of the first von Mises distribution
        vonmises_kappa1 = torch.exp(x[:, 7]).unsqueeze(0).unsqueeze(0)
        vonmises_kappa1 = vonmises_kappa1.repeat(self.image_dim, self.image_dim, 1)
        vonmises_kappa1 = vonmises_kappa1.permute(2, 0, 1)

        # vonmises_weight1 = torch.exp(x[:, 8]).unsqueeze(0).unsqueeze(0)
        vonmises_weight1 = x[:, 8].unsqueeze(0).unsqueeze(0)
        vonmises_weight1 = vonmises_weight1.repeat(self.image_dim, self.image_dim, 1)
        vonmises_weight1 = vonmises_weight1.permute(2, 0, 1)

        # vm_mu and weight for the second von Mises distribution
        bearing_new2 = x[:, 9] + bearing[:, 0]

        vonmises_mu2 = bearing_new2.unsqueeze(0).unsqueeze(0)
        vonmises_mu2 = vonmises_mu2.repeat(self.image_dim, self.image_dim, 1)
        vonmises_mu2 = vonmises_mu2.permute(2, 0, 1)

        # parameters of the second von Mises distribution
        vonmises_kappa2 = torch.exp(x[:, 10]).unsqueeze(0).unsqueeze(0)
        vonmises_kappa2 = vonmises_kappa2.repeat(self.image_dim, self.image_dim, 1)
        vonmises_kappa2 = vonmises_kappa2.permute(2, 0, 1)

        # vonmises_weight2 = torch.exp(x[:, 11]).unsqueeze(0).unsqueeze(0)
        vonmises_weight2 = x[:, 11].unsqueeze(0).unsqueeze(0)
        vonmises_weight2 = vonmises_weight2.repeat(self.image_dim, self.image_dim, 1)
        vonmises_weight2 = vonmises_weight2.permute(2, 0, 1)

        # Apply softmax to the weights
        vonmises_weights = torch.stack([vonmises_weight1, vonmises_weight2], dim=0)
        vonmises_weights = torch.nn.functional.softmax(vonmises_weights, dim=0)
        vonmises_weight1 = vonmises_weights[0]
        vonmises_weight2 = vonmises_weights[1]

        # calculation of von Mises densities
        vonmises_density_layer1 = self.vonmises_density(self.bearing_layer,
                                                        vonmises_kappa1,
                                                        vonmises_mu1).to(device)

        vonmises_density_layer2 = self.vonmises_density(self.bearing_layer,
                                                        vonmises_kappa2,
                                                        vonmises_mu2).to(device)

        # combining both densities to create a mixture distribution using the logsumexp trick
        logsumexp_vm_corr = torch.max(vonmises_density_layer1, vonmises_density_layer2)
        vonmises_density_layer = logsumexp_vm_corr + torch.log(vonmises_weight1 * torch.exp(vonmises_density_layer1 - logsumexp_vm_corr) +
                                                               vonmises_weight2 * torch.exp(vonmises_density_layer2 - logsumexp_vm_corr))
        # print(torch.sum(vonmises_density_layer))
        # print(torch.sum(torch.exp(vonmises_density_layer)))

        # combining the two distributions
        movement_grid = gamma_density_layer + vonmises_density_layer # Gamma and von Mises densities are on the log-scale

        # normalise (on the log-scale using the log-sum-exp trick) before combining with the habitat predictions
        movement_grid = movement_grid - torch.logsumexp(movement_grid, dim = (1, 2), keepdim = True)
        # print('Movement grid norm ', torch.sum(movement_grid))
        # print(torch.sum(torch.exp(movement_grid)))

        return movement_grid

Scalar to grid block

This block takes any scalar value (e.g., time of day, day of year) and converts it to a 2D image, with the same values for all pixels.

This is so that the scalar values can be used as input to the convolutional layers.

class Scalar_to_Grid_Block(nn.Module):
    def __init__(self, params):
        super(Scalar_to_Grid_Block, self).__init__()

        # define the parameters
        self.batch_size = params.batch_size
        self.image_dim = params.image_dim
        self.device = params.device

    def forward(self, x):

        # how many scalar values are being passed in
        num_scalars = x.shape[1]
        # expand the scalar values to the spatial dimensions of the image
        scalar_map = x.view(x.shape[0], num_scalars, 1, 1).expand(x.shape[0],
                                                                  num_scalars,
                                                                  self.image_dim,
                                                                  self.image_dim)

        # return the scalar maps
        return scalar_map

Combine the blocks into the deepSSF model

Here is where we combine the blocks into a model. Similarly to the previous blocks, the model is a Python class that inherits from torch.nn.Module, which combines other torch.nn.Module modules.

For example, we can instantiate the habitat selection convolution block using self.conv_habitat = Conv2d_block_spatial(params) in the __init__ method (the ‘constructor’ for a class). We can now access that block using self.conv_habitat in the forward method.

In the forward method, we pass the input data through the habitat selection convolution block using output_habitat = self.conv_habitat(all_spatial), where all_spatial is the input data, which is a combination of the spatial covariates and the scalar values converted to 2D images.

First we instantiate the blocks, and then define the forward method, which defines the data flow through the network during inference or training.

class ConvJointModel(nn.Module):
    def __init__(self, params):
        """
        ConvJointModel:
        - Initializes blocks for scalar-to-grid transformation,
          habitat convolution, movement convolution + movement fully connected, and final parameter-to-grid transformation.
        - Accepts parameters from the params object, which we will define later.
        """
        super(ConvJointModel, self).__init__()

        # Block to convert scalar features into grid-like (spatial) features
        self.scalar_grid_output = Scalar_to_Grid_Block(params)

        # Convolutional block for habitat selection
        self.conv_habitat = Conv2d_block_spatial(params)

        # Convolutional block for movement extraction (output fed into fully connected layers)
        self.conv_movement = Conv2d_block_toFC(params)

        # Fully connected block for movement
        self.fcn_movement_all = FCN_block_all_movement(params)

        # Converts movement distribution parameters into a grid (the 2D movement kernel)
        self.movement_grid_output = Params_to_Grid_Block(params)

        # Device information from params (e.g., CPU or GPU)
        self.device = params.device

    def forward(self, x):
        """
        Forward pass:
        1. Extract scalar data and convert to grid features.
        2. Concatenate the newly created scalar-based grids with spatial data.
        3. Pass this combined input through separate sub-networks for habitat and movement.
        4. Convert movement parameters to a grid, then stack the habitat and movement outputs.
        """
        # x contains:
        # - spatial_data_x (image-like layers)
        # - scalars_to_grid (scalar features needing conversion)
        # - bearing_x (the bearing from the previous time step, the turning angle is estimated as the deviation from this)
        spatial_data_x = x[0]
        scalars_to_grid = x[1]
        bearing_x = x[2]

        # Convert scalar data to spatial (grid) form
        scalar_grids = self.scalar_grid_output(scalars_to_grid)

        # Combine the original spatial data with the newly generated scalar grids
        all_spatial = torch.cat([spatial_data_x, scalar_grids], dim=1)

        # HABITAT SUBNETWORK
        # Convolutional feature extraction for habitat selection
        output_habitat = self.conv_habitat(all_spatial)

        # MOVEMENT SUBNETWORK
        # Convolutional feature extraction (different architecture for movement)
        conv_movement = self.conv_movement(all_spatial)

        # Fully connected layers for movement (processing both spatial features and any extras)
        output_movement = self.fcn_movement_all(conv_movement)

        # Transform the movement parameters into a grid, using bearing information
        output_movement = self.movement_grid_output(output_movement, bearing_x)

        # Combine (stack) habitat and movement outputs without merging them
        output = torch.stack((output_habitat, output_movement), dim=-1)

        return output

Set the parameters for the model which will be specified in a dictionary

This Python class serves as a simple parameter container for a model that involves both spatial (e.g., convolutional layers) and non-spatial inputs. It captures all relevant hyperparameters and settings—such as image dimensions, kernel sizes, and fully connected layer dimensions—along with the target device (CPU or GPU). This structure allows easy configuration of the model without scattering parameters throughout the code.

class ModelParams():
    def __init__(self, dict_params):
        self.batch_size = dict_params["batch_size"]
        self.image_dim = dict_params["image_dim"]
        self.pixel_size = dict_params["pixel_size"]
        self.dim_in_nonspatial_to_grid = dict_params["dim_in_nonspatial_to_grid"]
        self.dense_dim_in_nonspatial = dict_params["dense_dim_in_nonspatial"]
        self.dense_dim_hidden = dict_params["dense_dim_hidden"]
        self.dense_dim_in_all = dict_params["dense_dim_in_all"]
        self.input_channels = dict_params["input_channels"]
        self.output_channels = dict_params["output_channels"]
        self.output_channels_movement = dict_params["output_channels_movement"]
        self.kernel_size = dict_params["kernel_size"]
        self.stride = dict_params["stride"]
        self.kernel_size_mp = dict_params["kernel_size_mp"]
        self.stride_mp = dict_params["stride_mp"]
        self.padding = dict_params["padding"]
        self.image_dim = dict_params["image_dim"]
        self.num_movement_params = dict_params["num_movement_params"]
        self.dropout = dict_params["dropout"]
        self.device = dict_params["device"]

Define the parameters for the model

Here we enter the specific parameter values and hyperparameters for the model. These are the values that will be used to instantiate the model.

n_max_pool_layers = 2 # used to determine the number of inputs entering the fully connected block - needs to be manually changed if the number of max pooling layers is changed

params_dict = {"batch_size": batch_size, #number of samples in each batch
               "image_dim": 101, #number of pixels along the edge of each local patch/image
               "pixel_size": 25, #number of metres along the edge of a pixel
               "input_channels": 4 + 4, #number of spatial layers in each image + number of scalar layers that are converted to a grid
               "dim_in_nonspatial_to_grid": 4, #the number of scalar predictors that are converted to a grid and appended to the spatial features
               "dense_dim_in_nonspatial": 4, #change this to however many other scalar predictors you have (bearing, velocity etc)
               "kernel_size": 3, #the size of the 2D moving windows / kernels that are being learned
               "stride": 1, #the stride used when applying the kernel.  This reduces the dimension of the output if set to greater than 1
               "kernel_size_mp": 2, #the size of the kernel that is used in max pooling operations
               "stride_mp": 2, #the stride that is used in max pooling operations
               "padding": 1, #the amount of padding to apply to images prior to applying the 2D convolution
               "num_movement_params": 12, #number of parameters used to parameterise the movement kernel
               "dropout": 0.1, #the proportion of nodes that are dropped out in the dropout layers

               # hyperparameters that change the model architecture
               "output_channels": 4, #number of convolution filters to learn
               "output_channels_movement": 2, #number of convolution filters to learn for the movement kernel
               "dense_dim_hidden": 128, #number of nodes in the hidden layers

               # this will be updated below
               "dense_dim_in_all": 1, #number of inputs entering the fully connected block once the nonspatial features have been concatenated to the spatial features
               "device": device
               }

# Now update the dictionary with calculated values
params_dict["dense_dim_in_all"] = int(((params_dict["image_dim"] - (params_dict["image_dim"] % 2))**2) * (params_dict["output_channels_movement"] / (4**n_max_pool_layers)))

Note about the model

In future scripts (such as when simulating from the deepSSF model or training the model on Sentinel-2 data) we want to load the same model, so to prevent copying and pasting, we will save the model definition to a Python file and import it into future scripts.

We do this by copying the model definition above to a Python file named deepSSF_model.py, which can be imported into future scripts using import deepSSF_model.

Ideally we would just use that file to define the model in this script, but we include it here as it’s helpful to test components of the model and see how it all works.

Just remember that if you make changes to the model in this script, you will have to copy them across the deepSSF_model.py file, or just change the model definition in the deepSSF_model.py file and call that directly from here to train it.

To call it you would uncomment the lines in the next cell.

# # Import the functions in the deepSSF_model.py file
# import deepSSF_model

# # Create an instance of the ModelParams class using the params_dict
# params = deepSSF_model.ModelParams(deepSSF_model.params_dict)

# # Create an instance of the ConvJointModel class using the params
# model = deepSSF_model.ConvJointModel(params).to(device)

# # Print the model architecture to check that it worked
# print(model)

Instantiate the model

Here we instantiate the model using the parameters defined above.

# Initialize the parameter container using the parameters defined in 'params_dict'
params = ModelParams(params_dict)

# Create an instance of the ConvJointModel using the parameters,
# and move the model to the specified device (e.g., CPU or GPU)
model = ConvJointModel(params).to(device)

# Print the model architecture
print(model)

ConvJointModel(
  (scalar_grid_output): Scalar_to_Grid_Block()
  (conv_habitat): Conv2d_block_spatial(
    (conv2d): Sequential(
      (0): Conv2d(8, 4, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
      (1): ReLU()
      (2): Conv2d(4, 4, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
      (3): ReLU()
      (4): Conv2d(4, 1, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
    )
  )
  (conv_movement): Conv2d_block_toFC(
    (conv2d): Sequential(
      (0): Conv2d(8, 2, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
      (1): ReLU()
      (2): MaxPool2d(kernel_size=2, stride=2, padding=0, dilation=1, ceil_mode=False)
      (3): Conv2d(2, 2, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
      (4): ReLU()
      (5): MaxPool2d(kernel_size=2, stride=2, padding=0, dilation=1, ceil_mode=False)
      (6): Flatten(start_dim=1, end_dim=-1)
    )
  )
  (fcn_movement_all): FCN_block_all_movement(
    (ffn): Sequential(
      (0): Linear(in_features=1250, out_features=128, bias=True)
      (1): Dropout(p=0.1, inplace=False)
      (2): ReLU()
      (3): Linear(in_features=128, out_features=128, bias=True)
      (4): Dropout(p=0.1, inplace=False)
      (5): ReLU()
      (6): Linear(in_features=128, out_features=12, bias=True)
    )
  )
  (movement_grid_output): Params_to_Grid_Block()
)

Testing model components

As we’ve defined the model, we can now test the components to ensure that they are working as expected.

We can do this block by block, or we can test the entire model.

We’ll start by testing a few of the blocks.

Testing the movement parameter to probability distribution block

Change any of the values to try out different movement kernels.

# Create a bearing tensor (e.g., 0 radians) on the desired device
test_bearing = torch.tensor([[0.0]], device=device)

# Instantiate the Params_to_Grid_Block using the given parameters
test_block = Params_to_Grid_Block(params)

# Define the parameters for the movement density
# These are the parameters that the model will learn to predict

# First Gamma distribution
gamma_shape1 = 0.75
gamma_scale1 = 100
gamma_weight1 = 0.05

# Second Gamma distribution
gamma_shape2 = 3
gamma_scale2 = 500
gamma_scale2 = gamma_scale2 / 500 # divide this by 500 (as this is what the model predicts)
gamma_weight2 = 0.95

# First von Mises distribution
vonmises_mu1 = 0.0
vonmises_kappa1 = 0.1
vonmises_weight1 = 0.75

# Second von Mises distribution
vonmises_mu2 = -np.pi
vonmises_kappa2 = 0.01
vonmises_weight2 = 0.25

# Provide parameters on a log-scale (since the block exponentiates them internally)
# Here, each group of three values represents a Gamma distribution's shape, scale, and weight, respectively.
movement_density = test_block(
    torch.tensor(
        [[
            # Gamma 1
            np.log(gamma_shape1),   np.log(gamma_scale1),   gamma_weight1, # Gamma 1 shape, scale, and weight
            # Gamma 2
            np.log(gamma_shape2),   np.log(gamma_scale2),   gamma_weight2, # Gamma 2 shape, scale, and weight
            # von Mises 1
            vonmises_mu1,           np.log(vonmises_kappa1), vonmises_weight1, # von Mises 1 mu, kappa, and weight
            # von Mises 2
            vonmises_mu2,           np.log(vonmises_kappa2), vonmises_weight2 # von Mises 2 mu, kappa, and weight
        ]],
        device=device,
        dtype=torch.float32
    ),
    test_bearing
)

# Alternatively, if you had direct (non-log) values as the model sees them:
# movement_density = test_block(torch.tensor([[-.5, -.5, -.5, -.5]], device=device))

# print(movement_density)
print(movement_density.shape)

# Plot the resulting movement density as an image
plt.imshow(np.exp(movement_density.detach().cpu().numpy()[0]))
plt.colorbar()
plt.show()

print(np.sum(np.exp(movement_density.detach().cpu().numpy()[0])))

torch.Size([1, 101, 101])

0.99999976

Pull out some testing data

To test the other blocks, and the full model, we will need some data. We can pull that out from the training set.

# Number of samples in the train dataset
print("Number of samples in the train dataset: ", len(dataloader_train.dataset))
print('\n')

# Select an index from the test dataset to retrieve a sample, between 0 and number of samples
# We picked this fairly arbitrarily, but with some interesting environmental features to illustrate the model's predictions
iteration_index = 2700

# 2. Retrieve a single sample (features and label) from the test dataset at the specified index

# sample_spatial_covs is a sample of the spatial covariates for a single step
# sample_temporal_covs is a sample of the temporal covariates for a single step
# sample_prev_bearing is a sample bearing of the previous step
# sample_next_step is the target label (what we are trying to predict) for the next step

# We set these here and will also use them later in the script to check how the model's predictions look,
# and when we extract feature maps from the convolutional layers
sample_spatial_covs, sample_temporal_covs, sample_prev_bearing, sample_next_step = dataloader_train.dataset[iteration_index]

# 3. Reshape data tensors to add a batch dimension (since the model expects batches)
sample_spatial_covs = sample_spatial_covs.unsqueeze(0).to(device)
sample_temporal_covs = sample_temporal_covs.unsqueeze(0).to(device)
sample_prev_bearing = sample_prev_bearing.unsqueeze(0).to(device)
sample_next_step = sample_next_step.unsqueeze(0).to(device)

print(f'Shape of the sample spatial covariates:  {sample_spatial_covs.shape}')
print(f'Shape of the sample temporal covariates: {sample_temporal_covs.shape}')
print(f'Shape of the sample previous bearing:    {sample_prev_bearing.shape}')
print(f'Shape of the sample next step:           {sample_next_step.shape}')

Number of samples in the train dataset:  8082


Shape of the sample spatial covariates:  torch.Size([1, 4, 101, 101])
Shape of the sample temporal covariates: torch.Size([1, 4])
Shape of the sample previous bearing:    torch.Size([1, 1])
Shape of the sample next step:           torch.Size([1, 101, 101])

For visualisation, we can return the scale of the covariates to their original values.

# 1. NDVI (Normalized Difference Vegetation Index)
ndvi_norm = sample_spatial_covs.detach().cpu()[0, 0, :, :]
ndvi_natural = (ndvi_norm * (ndvi_max - ndvi_min)) + ndvi_min

# 2. Canopy cover
canopy_norm = sample_spatial_covs.detach().cpu()[0, 1, :, :]
canopy_natural = (canopy_norm * (canopy_max - canopy_min)) + canopy_min

# 3. Herbaceous vegetation
herby_norm = sample_spatial_covs.detach().cpu()[0, 2, :, :]
herby_natural = (herby_norm * (herby_max - herby_min)) + herby_min

# 4. Slope
slope_norm = sample_spatial_covs.detach().cpu()[0, 3, :, :]
slope_natural = (slope_norm * (slope_max - slope_min)) + slope_min

Pull out the scalar values

# Convert the PyTorch tensor x2 to a NumPy array:
#   1) Detach from the computation graph so no gradients are tracked.
#   2) Move to CPU memory.
#   3) Convert to NumPy.
# Then extract the first sample (index 0) and its respective channel for each variable:
hour_t2_sin = sample_temporal_covs.detach().cpu().numpy()[0, 0]
hour_t2_cos = sample_temporal_covs.detach().cpu().numpy()[0, 1]
yday_t2_sin = sample_temporal_covs.detach().cpu().numpy()[0, 2]
yday_t2_cos = sample_temporal_covs.detach().cpu().numpy()[0, 3]

# Convert x3 similarly and extract the bearing from the first sample and channel:
bearing = sample_prev_bearing.detach().cpu().numpy()[0, 0]

Helper functions

To return the hour and day of the year to their original values, we can use the following functions.

def recover_hour(sin_term, cos_term):
    # Calculate the angle theta
    theta = np.arctan2(sin_term, cos_term)
    # Calculate hour_t2
    hour = (12 * theta) / np.pi % 24
    return hour

def recover_yday(sin_term, cos_term):
    # Calculate the angle theta
    theta = np.arctan2(sin_term, cos_term)
    # Calculate hour_t2
    yday = (365 * theta) / (2 * np.pi)  % 365
    return yday

Calculate the hour, day of year and previous bearing of the test sample

hour_t2 = recover_hour(hour_t2_sin, hour_t2_cos)
hour_t2_integer = int(hour_t2)  # Convert to integer
print(f'Hour:               {hour_t2_integer}')

yday_t2 = recover_yday(yday_t2_sin, yday_t2_cos)
yday_t2_integer = int(yday_t2)  # Convert to integer
print(f'Day of the year:    {yday_t2_integer}')

bearing_degrees = np.degrees(bearing) % 360
bearing_degrees = round(bearing_degrees, 1)  # Round to 2 decimal places
bearing_degrees = int(bearing_degrees)  # Convert to integer
print(f'Bearing (radians):  {bearing}')
print(f'Bearing (degrees):  {bearing_degrees}')

Hour:               17
Day of the year:    324
Bearing (radians):  2.7272613048553467
Bearing (degrees):  156

Grab the row and column of the observed next step (label or target)

# Find the coordinates of the element that is 1
target = sample_next_step.detach().cpu().numpy()[0,:,:]
coordinates = np.where(target == 1)
# Extract the coordinates
row, column = coordinates[0][0], coordinates[1][0]
print(f"The location of the next step is (row, column): ({row}, {column})")

The location of the next step is (row, column): (46, 44)

Plot the sample covariates

# Plot the covariates
fig, axs = plt.subplots(2, 2, figsize=(9, 7.5))

# Plot NDVI
im1 = axs[0, 0].imshow(ndvi_natural.numpy(), cmap='viridis')
axs[0, 0].set_title('NDVI')
fig.colorbar(im1, ax=axs[0, 0])

# Plot Canopy cover
im2 = axs[0, 1].imshow(canopy_natural.numpy(), cmap='viridis')
axs[0, 1].set_title('Canopy cover')
fig.colorbar(im2, ax=axs[0, 1])

# Plot Herbaceous vegetation
im3 = axs[1, 0].imshow(herby_natural.numpy(), cmap='viridis')
axs[1, 0].set_title('Herbaceous vegetation')
fig.colorbar(im3, ax=axs[1, 0])

# Plot Slope
im4 = axs[1, 1].imshow(slope_natural.numpy(), cmap='viridis')
axs[1, 1].set_title('Slope')
fig.colorbar(im4, ax=axs[1, 1])

filename_covs = f'{output_dir}/covs_id{buffalo_id}_yday{yday_t2_integer}_hour{hour_t2_integer}_bearing{bearing_degrees}_next_r{row}_c{column}.png'
plt.tight_layout()
plt.savefig(filename_covs, dpi=300, bbox_inches='tight') # if we want to save the figure
plt.show()
plt.close()  # Close the figure to free memory

Plot the target (observed location of the next step)

The model is trying to maximise the probability at the location of the next step, which is the target.

filename_target = f'{output_dir}/target_id{buffalo_id}_yday{yday_t2_integer}_hour{hour_t2_integer}_bearing{bearing_degrees}_next_r{row}_c{column}.png'

plt.imshow(target)
plt.tight_layout()
plt.savefig(filename_target, dpi=300, bbox_inches='tight') # if we want to save the figure
plt.colorbar()
plt.show()

Testing the scalar to grid function

This should just create a grid with the same value for all pixels.

# x2 contains the scalar inputs
print(sample_temporal_covs.shape)  # Check the shape of the scalar input
print(sample_temporal_covs[0, :])  # Print out the first set of scalars

# Create an instance of the scalar-to-grid block using model parameters
test_block = Scalar_to_Grid_Block(params)

# Convert scalars into spatial grid representation
scalar_maps = test_block(sample_temporal_covs).detach().cpu()
# print(scalar_maps)  # Optionally, to inspect raw output
print(scalar_maps.shape)  # Check the shape of the generated spatial maps

# Visualize one channel of the first example's scalar map
# (Values are should be repeated across the grid for each scalar)
scalar_index = 2  # Change this index to visualize other scalar maps
plt.imshow(scalar_maps[0, scalar_index]) # change the second index to see the other scalar maps
plt.colorbar()
plt.clim(-1, 1) # Set the color limits to match the range of the scalar values (sine and cosine of temporal parameters)
plt.text(scalar_maps.shape[2] // 2, scalar_maps.shape[3] // 2,
         f'Value: {round(sample_temporal_covs[0, scalar_index].item(), 2)}',
         ha='center', va='center', color='white', fontsize=12)
plt.show()

torch.Size([1, 4])
tensor([-0.9673, -0.2538, -0.6422,  0.7665])
torch.Size([1, 4, 101, 101])

Test the full model

The model is initialised with random weights and hasn’t been trained yet, so the output will not be meaningful. However, we can check that the model runs without errors and that the output is the correct shape.

# Put the model in evaluation mode (affects layers like dropout, batch norm, etc.)
model.eval()

# Pass the data through the model
test = model((sample_spatial_covs, sample_temporal_covs, sample_prev_bearing))

# Print the shape of the output
print(test.shape)

torch.Size([1, 101, 101, 2])

Habitat predictions

plt.imshow(test.detach().cpu().numpy()[0,:,:,0])
plt.colorbar()
plt.show()

Movement predictions

# Print input bearing
print(f'Bearing (radians): {bearing}')
print(f'Bearing (degrees): {bearing_degrees}')

# Plot the movement density on the log-scale
plt.imshow(test.detach().cpu().numpy()[0,:,:,1])
plt.colorbar()
plt.show()

# Plot the movement density on the natural scale
plt.imshow(np.exp(test.detach().cpu().numpy()[0,:,:,1]))
plt.colorbar()
plt.show()

Bearing (radians): 2.7272613048553467
Bearing (degrees): 156

Next-step probability distribution

# Combine the habitat selection and the movement probabilities (unnormalised)
next_step = (test[:, :, :, 0] + test[:, :, :, 1])

# Plot the combined output
plt.imshow(next_step.detach().cpu().numpy()[0,:,:])
plt.colorbar()
plt.show()

Prepare for training

Loss function

We use a custom negative log likelihood loss function. Essentially what this does is extracts the next-step log-probability at the location of the observed next step, and then takes the negative of this value. This is the loss that we want to minimise, as we want to maximise the probability of the observed next step.

We will also save this loss function in a script called deepSSF_loss.py so that we can import it into future scripts.

class negativeLogLikeLoss(nn.Module):
    """
    Custom negative log-likelihood loss that operates on a 4D prediction tensor
    (batch, height, width, channels). The forward pass:
    1. Sums across channel 3 (two log-densities, habitat selection and movement predictions) to obtain a combined log-density.
    2. Multiplies this log-density by the target, which is 0 everywhere except for at the location of the next step, effectively extracting that value,
    then multiplies by -1 such that the function can be minimised (and the probabilities maximised).
    3. Applies the user-specified reduction (mean, sum, or none).
    """

    def __init__(self, reduction='mean'):
        """
        Args:
            reduction (str): Specifies the reduction to apply to the output:
                             'mean', 'sum', or 'none'.
        """
        super(negativeLogLikeLoss, self).__init__()
        assert reduction in ['mean', 'sum', 'none'], \
            "reduction should be 'mean', 'sum', or 'none'"
        self.reduction = reduction

    def forward(self, predict, target):
        """
        Forward pass of the negative log-likelihood loss.

        Args:
            predict (Tensor): A tensor of shape (B, H, W, 2) with log-densities
                              across two channels to be summed.
            target  (Tensor): A tensor of the same spatial dimensions (B, H, W)
                              indicating where the log-densities should be evaluated.

        Returns:
            Tensor: The computed negative log-likelihood loss. Shape depends on
                    the reduction method.
        """

        habitat_probability_surface = predict[:, :, :, 0] 
        movement_probability_surface = predict[:, :, :, 1] 

        # Sum the log-densities from the two channels
        predict_prod = habitat_probability_surface + movement_probability_surface

        # Check for NaNs in the combined predictions
        if torch.isnan(predict_prod).any():
            print("NaNs detected in predict_prod")
            print("predict_prod:", predict_prod)
            raise ValueError("NaNs detected in predict_prod")

        # Normalise the next-step log-densities using the log-sum-exp trick
        # predict_prod = predict_prod - torch.logsumexp(predict_prod, dim = (1, 2), keepdim = True)

        # Compute negative log-likelihood by multiplying log-densities with target
        # and then flipping the sign
        negLogLike = -1 * (predict_prod * target)

        # Check for NaNs after computing negative log-likelihood
        if torch.isnan(negLogLike).any():
            print("NaNs detected in negLogLike")
            print("negLogLike:", negLogLike)
            raise ValueError("NaNs detected in negLogLike")

        # Just extract the value at the next step
        negLogLike = negLogLike.sum(dim=(1, 2))

        # Calculate the loss on the habitat selection surface
        habitat_loss = -1 * (habitat_probability_surface * target)
        habitat_loss = habitat_loss.sum(dim=(1, 2))

        # Calculate the loss on the movement surface
        movement_loss = -1 * (movement_probability_surface * target)
        movement_loss = movement_loss.sum(dim=(1, 2))

        # Apply the specified reduction
        if self.reduction == 'mean':
            return torch.mean(negLogLike), torch.mean(habitat_loss), torch.mean(movement_loss)
        elif self.reduction == 'sum':
            return torch.sum(negLogLike), torch.sum(habitat_loss), torch.sum(movement_loss)
        elif self.reduction == 'none':
            return negLogLike, habitat_loss, movement_loss

        # Default return (though it should never reach here without hitting an if)
        return negLogLike, habitat_loss, movement_loss

Test the loss function

# Define the negative log-likelihood loss function with mean reduction
loss_fn = negativeLogLikeLoss(reduction='mean')

# Calculate the loss using the model outputs and the targets
total_loss, habitat_loss, movement_loss = loss_fn(model((sample_spatial_covs, sample_temporal_covs, sample_prev_bearing)), sample_next_step)
print(f'Total loss:     {total_loss}')
print(f'Habitat loss:   {habitat_loss}')
print(f'Movement loss:  {movement_loss}')

Total loss:     15.513991355895996
Habitat loss:   9.28042984008789
Movement loss:  6.2335615158081055

Early stopping code

This code will be used to stop training if the validation loss does not improve after a certain number of epochs.

When the loss of the validation data (which is held out from the training data) decreases (i.e. the model improves), the model weights are saved. Each time the validation loss does not decrease, a counter is incremented. If the counter reaches the patience value, the training loop will break and the model will stop training. The ‘final’ model is then the model that had the lowest validation loss.

We have saved this code in a script called deepSSF_early_stopping.py so that we can import it into future scripts.

class EarlyStopping:
    def __init__(self, patience=5, verbose=False, delta=0, path='checkpoint.pt', trace_func=print):
        """
        Args:
            patience (int): How long to wait after last time validation loss improved.
                            Default: 5
            verbose (bool): If True, prints a message for each validation loss improvement.
                            Default: False
            delta (float): Minimum change in the monitored quantity to qualify as an improvement.
                            Default: 0
            path (str): Path for the checkpoint to be saved to.
                            Default: 'checkpoint.pt'
            trace_func (function): trace print function.
                            Default: print
        """

        self.patience = patience
        self.verbose = verbose
        self.counter = 0
        self.best_score = None
        self.early_stop = False
        self.val_loss_min = float('inf')
        self.delta = delta
        self.path = path
        self.trace_func = trace_func

    def __call__(self, val_loss, model):

        # takes the validation loss and the model as inputs
        score = -val_loss

        # save the model's weights if the validation loss decreases
        if self.best_score is None:
            self.best_score = score
            self.save_checkpoint(val_loss, model)

        # if the validation loss does not decrease, increment the counter
        elif score < self.best_score + self.delta:
            self.counter += 1
            self.trace_func(f'EarlyStopping counter: {self.counter} out of {self.patience}')
            if self.counter >= self.patience:
                self.early_stop = True
        else:
            self.best_score = score
            self.save_checkpoint(val_loss, model)
            self.counter = 0

    def save_checkpoint(self, val_loss, model):
        '''Saves model when validation loss decrease.'''
        if self.verbose:
            self.trace_func(f'Validation loss decreased ({self.val_loss_min:.6f} --> {val_loss:.6f}).  Saving model ...')
        torch.save(model.state_dict(), self.path)
        self.val_loss_min = val_loss

Path to save the model weights

path_save_weights = f'{output_dir}/checkpoint_deepSSF_buffalo{buffalo_id}.pt'
print(path_save_weights)

../Python/outputs/model_training/id2005_deepSSF_training_1_2025-07-09/checkpoint_deepSSF_buffalo2005.pt

Set the learning rate

The learning rate is a hyperparameter that controls how much we are adjusting the weights of our network with respect to the loss gradient. A larger learning rate means that the optimiser will take larger steps, but it may overshoot the minimum. A smaller learning rate means that the optimiser will take smaller steps, but it may take a long time to converge.

We can therefore use an adaptive learning rate, which will adjust the learning rate during training. If the loss does not decrease after a certain number of epochs (also called the patience), the learning rate will be reduced by a factor of 10.

The patience of the learning rate should be less than the patience of the early stopping code, as we want to reduce the learning rate before we stop training.

# Set the initial learning rates for each process
initial_learning_rate_movement = 1e-5
initial_learning_rate_habitat = 1e-4

Instantiate the loss function, optimiser, learning rate scheduler and early stopping code

In this chunk, we set up all the components needed for training and validating a neural network model:

1. Loss Function: Uses a negative log-likelihood loss (negativeLogLikeLoss) with mean reduction.
2. Optimizer: Implements the Adam optimization algorithm, updating model parameters based on the computed gradients and a specified learning rate.
3. Scheduler: Automatically reduces the learning rate when the monitored metric (e.g., validation loss) stops improving.
4. Early Stopping: Monitors validation performance and stops training if the metric fails to improve after a certain number of epochs (patience).

# Define the negative log-likelihood loss function with mean reduction
loss_fn = negativeLogLikeLoss(reduction='mean')

# Create a combined optimiser for all movement-related parameters
movement_params = list(model.conv_movement.parameters()) + list(model.fcn_movement_all.parameters())

# Define separate optimizers for each component
optimiser_movement = optim.Adam(movement_params, lr=initial_learning_rate_movement)
optimiser_habitat = optim.Adam(model.conv_habitat.parameters(), lr=initial_learning_rate_habitat)

# Put optimisers into a tuple to call in the training loop
optimisers = (optimiser_movement, optimiser_habitat)

# Create separate schedulers for each optimizer
scheduler_movement = torch.optim.lr_scheduler.ReduceLROnPlateau(
    optimiser_movement, 'min', factor=0.1, patience=5)
scheduler_habitat = torch.optim.lr_scheduler.ReduceLROnPlateau(
    optimiser_habitat, 'min', factor=0.1, patience=5)

# EarlyStopping stops training after 'patience=10' epochs with no improvement,
#    optionally saving the best model weights
early_stopping = EarlyStopping(patience=25, verbose=True, path=path_save_weights)

Training loop

This code defines the main training loop for a single epoch. It iterates over batches from the training dataloader, moves the data to the correct device (e.g., CPU or GPU), calculates the loss, and performs backpropagation to update the model parameters. It also prints periodic updates of the current loss.

def train_loop(dataloader_train, model, loss_fn, optimisers, skip_epoch0_training=False):
    """
    Runs the training process for one epoch using the given dataloader, model,
    loss function, and optimizer. Prints progress updates every few batches.
    """

    # Unpack optimisers
    optimiser_movement, optimiser_habitat = optimisers

    # 1. Total number of training examples
    num_train_batches = len(dataloader_train)
    size = len(dataloader_train.dataset)

    # 2. Put model in training mode (affects layers like dropout, batchnorm)
    model.train()

    # 3. Variable to accumulate the total loss over the epoch
    epoch_loss = 0.0

    # 4. Loop over batches in the training dataloader
    for batch, (x1, x2, x3, y) in enumerate(dataloader_train):

        # Move the batch of data to the specified device (CPU/GPU)
        x1 = x1.to(device)
        x2 = x2.to(device)
        x3 = x3.to(device)
        y = y.to(device)

        # Forward pass: compute the model output and loss
        with torch.set_grad_enabled(not skip_epoch0_training):
            outputs = model((x1, x2, x3))
            total_loss, habitat_loss, movement_loss = loss_fn(outputs, y)

        epoch_loss += total_loss.detach()  # Use detach to prevent memory leaks

        # Only perform optimization if not skipping training
        if not skip_epoch0_training:
            # Backpropagation: compute gradients and update parameters
            # Reset gradients before the next iteration

            # Zero all gradients
            optimiser_movement.zero_grad()
            optimiser_habitat.zero_grad()

            # Single backward pass on total loss
            total_loss.backward()

            # For movement optimizer: save habitat gradients, then zero them out
            habitat_grads = []
            for param in model.conv_habitat.parameters():
                # Save the gradient
                if param.grad is not None:
                    habitat_grads.append(param.grad.clone())
                else:
                    habitat_grads.append(None)
                # Zero out habitat gradient for movement update
                param.grad = None

            # Update movement parameters
            optimiser_movement.step()

            # For habitat optimizer: restore habitat gradients and zero movement gradients
            for param in model.conv_movement.parameters():
                param.grad = None
            for param in model.fcn_movement_all.parameters():
                param.grad = None

            # Restore habitat gradients
            for i, param in enumerate(model.conv_habitat.parameters()):
                param.grad = habitat_grads[i]

            # Update habitat parameters
            optimiser_habitat.step()

        # Print an update every 5 batches to keep track of training progress
        if batch % 20 == 0:
            loss_val = total_loss.item()
            current = batch * batch_size + len(x1)
            if skip_epoch0_training:
                print(f"[Observation only] loss: {loss_val:>15f}  [{current:>5d}/{size:>5d}]")
            else:
                print(f"loss: {loss_val:>15f}  [{current:>5d}/{size:>5d}]")

        torch.cuda.empty_cache()

    # Compute the average training loss and print it
    epoch_loss /= num_train_batches
    if skip_epoch0_training:
        print(f"\nAvg training loss (observation only): {epoch_loss:>15f}")
    else:
        print(f"\nAvg training loss: {epoch_loss:>15f}")
    train_losses.append(epoch_loss.item())

Test loop

The test loop is similar to the training loop, but it does not perform backpropagation. It calculates the loss on the test set and returns the average loss.

def test_loop(dataloader_test, model, loss_fn):
    """
    Evaluates the model on the provided test dataset by computing
    the average loss over all batches.
    No gradients are computed during this process (torch.no_grad()).
    """

    # 1. Set the model to evaluation mode (affects layers like dropout, batchnorm).
    model.eval()

    size = len(dataloader_test.dataset)
    num_batches = len(dataloader_test)

    test_loss = 0

    # 2. Disable gradient computation to speed up evaluation and reduce memory usage
    with torch.no_grad():
        # 3. Loop through each batch in the test dataloader
        for x1, x2, x3, y in dataloader_test:

            # Move the batch of data to the appropriate device (CPU/GPU)
            # x1, x2, x3 are the spatial covariates, temporal covariates, and bearing, respectively
            # y is the label (observed location of the next step)
            x1 = x1.to(device)
            x2 = x2.to(device)
            x3 = x3.to(device)
            y = y.to(device)

            # Compute the loss on the test set (no backward pass needed)
            total_loss, habitat_loss, movement_loss = loss_fn(model((x1, x2, x3)), y)
            test_loss += total_loss.detach()

    # 4. Compute average test loss over all batches
    test_loss /= num_batches

    torch.cuda.empty_cache()

    # Print the average test loss
    print(f"Avg test loss:    {test_loss:>15f} \n")

Train the model

Here we have the main training process that loops over multiple epochs. Each epoch involves:

1. Training the model on a training dataset.
2. Validating the model on a validation dataset to monitor its performance and adjust the learning rate (via scheduler).
3. Checking for early stopping conditions. If triggered, the best model weights are restored, and a test evaluation is performed.

Additionally, commented-out code at the end shows how you might visualise and save intermediate training results (such as predicted probability surfaces) for diagnostic or research purposes. The saved images can then be combined into an animation.

epochs = 120
train_losses = []  # Track training losses across epochs
val_losses = []   # Track validation losses across epochs
val_habitat_losses = []  # Track validation habitat losses across epochs
val_movement_losses = []  # Track validation movement losses across epochs

# Initialize the parameter container using the parameters defined in 'params_dict'
params = ModelParams(params_dict)
# Create an instance of the ConvJointModel using the parameters,
# and move the model to the specified device (e.g., CPU or GPU)
model = ConvJointModel(params).to(device)
# Print the model architecture
print(model)

# Define the negative log-likelihood loss function with mean reduction
loss_fn = negativeLogLikeLoss(reduction='mean') #, alpha=0.5

# Set the initial learning rates for each process
initial_learning_rate_movement = 1e-5
initial_learning_rate_habitat = 1e-4

# Create a combined optimiser for all movement-related parameters
movement_params = list(model.conv_movement.parameters()) + list(model.fcn_movement_all.parameters())

# Define separate optimizers for each component
optimiser_movement = optim.Adam(movement_params, lr=initial_learning_rate_movement)
optimiser_habitat = optim.Adam(model.conv_habitat.parameters(), lr=initial_learning_rate_habitat)

# Put optimisers into a tuple to call in the training loop
optimisers = (optimiser_movement, optimiser_habitat)

# Create separate schedulers for each optimizer
scheduler_movement = torch.optim.lr_scheduler.ReduceLROnPlateau(
    optimiser_movement, 'min', factor=0.1, patience=5)
scheduler_habitat = torch.optim.lr_scheduler.ReduceLROnPlateau(
    optimiser_habitat, 'min', factor=0.1, patience=5)

# Reinitialise early stopping
early_stopping = EarlyStopping(patience=10, verbose=True, path=path_save_weights)

# Create directory for saving training images
os.makedirs(f'{output_dir}/training_images', exist_ok=True)

for t in range(epochs):

    # Initialise variables to store during training
    train_loss = 0.0
    num_train_batches = len(dataloader_train)

    val_loss = 0.0
    val_loss_habitat = 0.0
    val_loss_movement = 0.0
    num_batches = len(dataloader_val)

    print(f"Epoch {t+1}\n-------------------------------")

    # Skip training in the first epoch, but still calculate losses
    skip_training = (t == 0)

    # 1. Run the training loop for one epoch using the training dataloader
    train_loop(dataloader_train, model, loss_fn, optimisers, skip_epoch0_training=skip_training)

    # 2. Evaluate model performance on the validation dataset
    model.eval()  # Switch to evaluation mode for proper layer behavior
    with torch.no_grad():

        # Loop through each batch in the validation dataloader
        for x1, x2, x3, y in dataloader_val:
            # Move data to the chosen device (CPU/GPU)
            x1 = x1.to(device)
            x2 = x2.to(device)
            x3 = x3.to(device)
            y = y.to(device)

            # Accumulate validation loss
            total_loss, habitat_loss, movement_loss = loss_fn(model((x1, x2, x3)), y)
            val_loss += total_loss.detach()
            val_loss_habitat += habitat_loss.detach()
            val_loss_movement += movement_loss.detach()

    # # 3. Step the scheduler based on the validation loss (adjusts learning rate if needed)
    # scheduler.step(val_loss)
    scheduler_movement.step(val_loss_movement)
    scheduler_habitat.step(val_loss_habitat)

    # 4. Compute the average validation loss and print it, along with the current learning rate
    val_loss /= num_batches
    val_loss_habitat /= num_batches
    val_loss_movement /= num_batches

    print(f"Avg validation loss:            {val_loss:>15f}")
    print(f"Avg validation habitat loss:    {val_loss_habitat:>15f}")
    print(f"Avg validation movement loss:   {val_loss_movement:>15f}")
    print(f"Movement learning rate:         {scheduler_movement.get_last_lr()}")
    print(f"Habitat learning rate:          {scheduler_habitat.get_last_lr()}")

    # 5. Track the validation loss for plotting or monitoring
    val_losses.append(val_loss.item())
    val_habitat_losses.append(val_loss_habitat.item())
    val_movement_losses.append(val_loss_movement.item())

    # Memory management - add after validation but before early stopping check
    # torch.cuda.empty_cache()
    # gc.collect()

    # 6. Early stopping: if no improvement in validation loss for a set patience, stop training
    early_stopping(val_loss, model)
    if early_stopping.early_stop:
        print("Early stopping")
        # Restore the best model weights saved by EarlyStopping
        model.load_state_dict(torch.load(path_save_weights, weights_only=True, map_location=device))
        test_loop(dataloader_test, model, loss_fn)  # Evaluate on test set once training stops
        break
    else:
        model.eval()
        print("\n")

    torch.cuda.empty_cache()


    # ----------------------------------------------------
    # The following code demonstrates how
    # to optionally visualize or save intermediate results
    # (e.g., habitat probability surface, movement probability,
    # and next-step probability surfaces).

    # uncomment the code all in one go to run it (it should be inside the training loop)
    # ----------------------------------------------------

    # Extract training and validation losses for plotting

    # Convert the list of tensors to a single tensor
    train_losses_np = torch.tensor(train_losses).detach().cpu().numpy()
    val_losses_np = torch.tensor(val_losses).detach().cpu().numpy()
    val_habitat_losses_np = torch.tensor(val_habitat_losses).detach().cpu().numpy()
    val_movement_losses_np = torch.tensor(val_movement_losses).detach().cpu().numpy()

    # Number of epochs
    n_epochs = len(val_losses)

    # -----------------------------------------------------------
    # 1. Retrieve a single test example (covariates and labels)
    #    at the specified 'iteration_index' from the test dataset
    # -----------------------------------------------------------
    x1, x2, x3, labels = dataloader_train.dataset[iteration_index]

    # -----------------------------------------------------------
    # 2. Add a batch dimension and move tensors to the device
    #    for model inference
    # -----------------------------------------------------------
    x1 = x1.unsqueeze(0).to(device)
    x2 = x2.unsqueeze(0).to(device)
    x3 = x3.unsqueeze(0).to(device)

    # -----------------------------------------------------------
    # 3. Run the model on the single test example
    # -----------------------------------------------------------
    test = model((x1, x2, x3))

    # -----------------------------------------------------------
    # 4. Extract habitat and movement outputs;
    #    convert them to NumPy arrays for visualization
    # -----------------------------------------------------------
    hab_density = test.detach().cpu().numpy()[0, :, :, 0]
    movement_density = test.detach().cpu().numpy()[0, :, :, 1]

    # -----------------------------------------------------------
    # 5. Generate masks to exclude certain border cells for
    #    color scale reasons (setting them to -inf).
    # -----------------------------------------------------------
    x_mask = np.ones_like(hab_density)
    y_mask = np.ones_like(hab_density)

    # Mask out a few columns (0-2 and 98-end) and rows (0-2 and 98-end)
    x_mask[:, :3] = -np.inf
    x_mask[:, 98:] = -np.inf
    y_mask[:3, :] = -np.inf
    y_mask[98:, :] = -np.inf

    # Apply the masks to habitat density
    hab_density_mask = hab_density * x_mask * y_mask

    # Combine habitat and movement densities to represent
    # next-step probability
    step_density = hab_density + movement_density
    step_density_mask = step_density * x_mask * y_mask

    # Plot the covariates
    fig, axs = plt.subplots(2, 2, figsize=(9, 7.5))

    # # Plot NDVI
    # im1 = axs[0, 0].imshow(ndvi_natural.numpy(), cmap='viridis')
    # axs[0, 0].set_title('NDVI')
    # fig.colorbar(im1, ax=axs[0, 0])

    # Plot Training and Validation Loss
    axs[0, 0].plot(range(n_epochs), train_losses_np, label='Training Loss', color='blue')
    axs[0, 0].plot(range(n_epochs), val_losses_np, label='Validation Loss', color='red')
    # axs[0, 0].plot(range(n_epochs), val_habitat_losses_np, label='Validation Habitat Loss', color='green')
    # axs[0, 0].plot(range(n_epochs), val_movement_losses_np, label='Validation Movement Loss', color='orange')
    axs[0, 0].set_xlim(0, 120)
    axs[0, 0].set_title('Training and validation loss')
    axs[0, 0].set_xlabel('Epoch')
    axs[0, 0].set_ylabel('Loss')
    axs[0, 0].legend()

    # Plot habitat selection log-probability
    im2 = axs[0, 1].imshow(hab_density_mask, cmap='viridis')
    axs[0, 1].set_title('Habitat log-probability')
    fig.colorbar(im2, ax=axs[0, 1])

    # Plot movement log-probability
    im3 = axs[1, 0].imshow(movement_density, cmap='viridis')
    axs[1, 0].set_title('Movement log-probability')
    fig.colorbar(im3, ax=axs[1, 0])

    # Plot next-step log-probability
    im4 = axs[1, 1].imshow(step_density_mask, cmap='viridis')
    axs[1, 1].set_title('Next-step log-probability')
    fig.colorbar(im4, ax=axs[1, 1])

    filename_covs = f'{output_dir}/training_images/training_epoch{t}_id{buffalo_id}_yday{yday_t2_integer}_hour{hour_t2_integer}_bearing{bearing_degrees}_next_r{row}_c{column}.png'
    plt.tight_layout()
    plt.savefig(filename_covs, dpi=300) # creates inconsistent image sizes >>> , bbox_inches='tight'
    # plt.show()
    plt.close()  # Close the figure to free memory

print("Done!")

ConvJointModel(
  (scalar_grid_output): Scalar_to_Grid_Block()
  (conv_habitat): Conv2d_block_spatial(
    (conv2d): Sequential(
      (0): Conv2d(8, 4, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
      (1): ReLU()
      (2): Conv2d(4, 4, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
      (3): ReLU()
      (4): Conv2d(4, 1, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
    )
  )
  (conv_movement): Conv2d_block_toFC(
    (conv2d): Sequential(
      (0): Conv2d(8, 2, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
      (1): ReLU()
      (2): MaxPool2d(kernel_size=2, stride=2, padding=0, dilation=1, ceil_mode=False)
      (3): Conv2d(2, 2, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
      (4): ReLU()
      (5): MaxPool2d(kernel_size=2, stride=2, padding=0, dilation=1, ceil_mode=False)
      (6): Flatten(start_dim=1, end_dim=-1)
    )
  )
  (fcn_movement_all): FCN_block_all_movement(
    (ffn): Sequential(
      (0): Linear(in_features=1250, out_features=128, bias=True)
      (1): Dropout(p=0.1, inplace=False)
      (2): ReLU()
      (3): Linear(in_features=128, out_features=128, bias=True)
      (4): Dropout(p=0.1, inplace=False)
      (5): ReLU()
      (6): Linear(in_features=128, out_features=12, bias=True)
    )
  )
  (movement_grid_output): Params_to_Grid_Block()
)
Epoch 1
-------------------------------
[Observation only] loss:       15.487193  [   32/ 8082]
[Observation only] loss:       15.850222  [  672/ 8082]
[Observation only] loss:       15.332795  [ 1312/ 8082]
[Observation only] loss:       15.690580  [ 1952/ 8082]
[Observation only] loss:       15.394622  [ 2592/ 8082]
[Observation only] loss:       15.047852  [ 3232/ 8082]
[Observation only] loss:       15.726194  [ 3872/ 8082]
[Observation only] loss:       15.706273  [ 4512/ 8082]
[Observation only] loss:       15.779350  [ 5152/ 8082]
[Observation only] loss:       15.651419  [ 5792/ 8082]
[Observation only] loss:       15.897856  [ 6432/ 8082]
[Observation only] loss:       15.443916  [ 7072/ 8082]
[Observation only] loss:       15.240913  [ 7712/ 8082]

Avg training loss (observation only):       15.575754
Avg validation loss:                  15.368755
Avg validation habitat loss:           9.229955
Avg validation movement loss:          6.138800
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (inf --> 15.368755).  Saving model ...


Epoch 2
-------------------------------
loss:       15.966509  [   32/ 8082]
loss:       16.146564  [  672/ 8082]
loss:       15.247894  [ 1312/ 8082]
loss:       15.961136  [ 1952/ 8082]
loss:       14.851363  [ 2592/ 8082]
loss:       15.272834  [ 3232/ 8082]
loss:       15.011624  [ 3872/ 8082]
loss:       15.370444  [ 4512/ 8082]
loss:       15.039049  [ 5152/ 8082]
loss:       15.737649  [ 5792/ 8082]
loss:       15.330084  [ 6432/ 8082]
loss:       15.872482  [ 7072/ 8082]
loss:       15.787772  [ 7712/ 8082]

Avg training loss:       15.514331
Avg validation loss:                  15.255189
Avg validation habitat loss:           9.227181
Avg validation movement loss:          6.028007
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (15.368755 --> 15.255189).  Saving model ...


Epoch 3
-------------------------------
loss:       15.674539  [   32/ 8082]
loss:       15.841510  [  672/ 8082]
loss:       15.836590  [ 1312/ 8082]
loss:       14.921206  [ 1952/ 8082]
loss:       15.138053  [ 2592/ 8082]
loss:       15.132509  [ 3232/ 8082]
loss:       15.444994  [ 3872/ 8082]
loss:       15.366295  [ 4512/ 8082]
loss:       15.505117  [ 5152/ 8082]
loss:       15.739033  [ 5792/ 8082]
loss:       15.300654  [ 6432/ 8082]
loss:       14.954218  [ 7072/ 8082]
loss:       15.804819  [ 7712/ 8082]

Avg training loss:       15.329888
Avg validation loss:                  15.095748
Avg validation habitat loss:           9.212420
Avg validation movement loss:          5.883327
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (15.255189 --> 15.095748).  Saving model ...


Epoch 4
-------------------------------
loss:       15.247113  [   32/ 8082]
loss:       14.966749  [  672/ 8082]
loss:       15.749109  [ 1312/ 8082]
loss:       15.845623  [ 1952/ 8082]
loss:       14.982266  [ 2592/ 8082]
loss:       15.168871  [ 3232/ 8082]
loss:       15.247325  [ 3872/ 8082]
loss:       15.815178  [ 4512/ 8082]
loss:       16.185966  [ 5152/ 8082]
loss:       15.230312  [ 5792/ 8082]
loss:       14.579098  [ 6432/ 8082]
loss:       15.381664  [ 7072/ 8082]
loss:       14.566433  [ 7712/ 8082]

Avg training loss:       15.165933
Avg validation loss:                  14.955804
Avg validation habitat loss:           9.155921
Avg validation movement loss:          5.799882
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (15.095748 --> 14.955804).  Saving model ...


Epoch 5
-------------------------------
loss:       14.800121  [   32/ 8082]
loss:       16.026505  [  672/ 8082]
loss:       14.700174  [ 1312/ 8082]
loss:       15.216954  [ 1952/ 8082]
loss:       15.160725  [ 2592/ 8082]
loss:       14.885765  [ 3232/ 8082]
loss:       15.173277  [ 3872/ 8082]
loss:       15.175369  [ 4512/ 8082]
loss:       15.083509  [ 5152/ 8082]
loss:       15.064681  [ 5792/ 8082]
loss:       14.418967  [ 6432/ 8082]
loss:       15.269823  [ 7072/ 8082]
loss:       15.114481  [ 7712/ 8082]

Avg training loss:       15.027689
Avg validation loss:                  14.828360
Avg validation habitat loss:           9.091925
Avg validation movement loss:          5.736434
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.955804 --> 14.828360).  Saving model ...


Epoch 6
-------------------------------
loss:       15.649390  [   32/ 8082]
loss:       14.294285  [  672/ 8082]
loss:       14.011728  [ 1312/ 8082]
loss:       15.519311  [ 1952/ 8082]
loss:       14.934125  [ 2592/ 8082]
loss:       14.906005  [ 3232/ 8082]
loss:       15.855635  [ 3872/ 8082]
loss:       14.978543  [ 4512/ 8082]
loss:       15.219209  [ 5152/ 8082]
loss:       15.395997  [ 5792/ 8082]
loss:       15.371474  [ 6432/ 8082]
loss:       15.567981  [ 7072/ 8082]
loss:       14.756718  [ 7712/ 8082]

Avg training loss:       14.955096
Avg validation loss:                  14.773286
Avg validation habitat loss:           9.063852
Avg validation movement loss:          5.709431
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.828360 --> 14.773286).  Saving model ...


Epoch 7
-------------------------------
loss:       14.679993  [   32/ 8082]
loss:       15.302557  [  672/ 8082]
loss:       14.944973  [ 1312/ 8082]
loss:       14.925349  [ 1952/ 8082]
loss:       15.247684  [ 2592/ 8082]
loss:       14.393382  [ 3232/ 8082]
loss:       15.381978  [ 3872/ 8082]
loss:       15.087513  [ 4512/ 8082]
loss:       14.903479  [ 5152/ 8082]
loss:       14.965706  [ 5792/ 8082]
loss:       14.814209  [ 6432/ 8082]
loss:       14.330323  [ 7072/ 8082]
loss:       14.276510  [ 7712/ 8082]

Avg training loss:       14.911941
Avg validation loss:                  14.698968
Avg validation habitat loss:           9.052184
Avg validation movement loss:          5.646785
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.773286 --> 14.698968).  Saving model ...


Epoch 8
-------------------------------
loss:       15.793880  [   32/ 8082]
loss:       14.910190  [  672/ 8082]
loss:       14.595140  [ 1312/ 8082]
loss:       14.614725  [ 1952/ 8082]
loss:       15.888414  [ 2592/ 8082]
loss:       13.564251  [ 3232/ 8082]
loss:       15.210400  [ 3872/ 8082]
loss:       15.971617  [ 4512/ 8082]
loss:       14.066001  [ 5152/ 8082]
loss:       15.758875  [ 5792/ 8082]
loss:       14.869189  [ 6432/ 8082]
loss:       14.751146  [ 7072/ 8082]
loss:       15.293427  [ 7712/ 8082]

Avg training loss:       14.880715
Avg validation loss:                  14.694047
Avg validation habitat loss:           9.046191
Avg validation movement loss:          5.647855
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.698968 --> 14.694047).  Saving model ...


Epoch 9
-------------------------------
loss:       14.884699  [   32/ 8082]
loss:       15.309649  [  672/ 8082]
loss:       15.117187  [ 1312/ 8082]
loss:       14.595469  [ 1952/ 8082]
loss:       15.653646  [ 2592/ 8082]
loss:       14.216175  [ 3232/ 8082]
loss:       14.956153  [ 3872/ 8082]
loss:       14.409956  [ 4512/ 8082]
loss:       14.535713  [ 5152/ 8082]
loss:       15.719199  [ 5792/ 8082]
loss:       14.425685  [ 6432/ 8082]
loss:       14.519854  [ 7072/ 8082]
loss:       14.875132  [ 7712/ 8082]

Avg training loss:       14.854919
Avg validation loss:                  14.660014
Avg validation habitat loss:           9.045353
Avg validation movement loss:          5.614663
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.694047 --> 14.660014).  Saving model ...


Epoch 10
-------------------------------
loss:       15.084428  [   32/ 8082]
loss:       14.334600  [  672/ 8082]
loss:       14.364571  [ 1312/ 8082]
loss:       14.219510  [ 1952/ 8082]
loss:       16.039957  [ 2592/ 8082]
loss:       15.185510  [ 3232/ 8082]
loss:       14.642065  [ 3872/ 8082]
loss:       14.612602  [ 4512/ 8082]
loss:       13.912735  [ 5152/ 8082]
loss:       14.630155  [ 5792/ 8082]
loss:       14.162275  [ 6432/ 8082]
loss:       14.451645  [ 7072/ 8082]
loss:       15.631488  [ 7712/ 8082]

Avg training loss:       14.838990
Avg validation loss:                  14.627380
Avg validation habitat loss:           9.039882
Avg validation movement loss:          5.587497
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.660014 --> 14.627380).  Saving model ...


Epoch 11
-------------------------------
loss:       16.076319  [   32/ 8082]
loss:       15.393507  [  672/ 8082]
loss:       14.155714  [ 1312/ 8082]
loss:       13.874891  [ 1952/ 8082]
loss:       15.566840  [ 2592/ 8082]
loss:       15.137169  [ 3232/ 8082]
loss:       15.217298  [ 3872/ 8082]
loss:       14.505178  [ 4512/ 8082]
loss:       14.603851  [ 5152/ 8082]
loss:       14.553517  [ 5792/ 8082]
loss:       14.246157  [ 6432/ 8082]
loss:       14.389214  [ 7072/ 8082]
loss:       14.694841  [ 7712/ 8082]

Avg training loss:       14.823997
Avg validation loss:                  14.592532
Avg validation habitat loss:           9.036983
Avg validation movement loss:          5.555549
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.627380 --> 14.592532).  Saving model ...


Epoch 12
-------------------------------
loss:       14.919466  [   32/ 8082]
loss:       14.745094  [  672/ 8082]
loss:       15.635286  [ 1312/ 8082]
loss:       14.690367  [ 1952/ 8082]
loss:       14.325845  [ 2592/ 8082]
loss:       14.733474  [ 3232/ 8082]
loss:       14.673908  [ 3872/ 8082]
loss:       15.327178  [ 4512/ 8082]
loss:       13.800572  [ 5152/ 8082]
loss:       15.478816  [ 5792/ 8082]
loss:       14.400972  [ 6432/ 8082]
loss:       14.855070  [ 7072/ 8082]
loss:       15.354678  [ 7712/ 8082]

Avg training loss:       14.813732
Avg validation loss:                  14.574661
Avg validation habitat loss:           9.035049
Avg validation movement loss:          5.539614
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.592532 --> 14.574661).  Saving model ...


Epoch 13
-------------------------------
loss:       14.599110  [   32/ 8082]
loss:       14.963902  [  672/ 8082]
loss:       13.945543  [ 1312/ 8082]
loss:       14.805696  [ 1952/ 8082]
loss:       13.888733  [ 2592/ 8082]
loss:       15.481693  [ 3232/ 8082]
loss:       14.108137  [ 3872/ 8082]
loss:       14.177662  [ 4512/ 8082]
loss:       14.585087  [ 5152/ 8082]
loss:       15.082737  [ 5792/ 8082]
loss:       15.406599  [ 6432/ 8082]
loss:       14.775091  [ 7072/ 8082]
loss:       15.089314  [ 7712/ 8082]

Avg training loss:       14.805429
Avg validation loss:                  14.569777
Avg validation habitat loss:           9.031487
Avg validation movement loss:          5.538291
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.574661 --> 14.569777).  Saving model ...


Epoch 14
-------------------------------
loss:       14.343162  [   32/ 8082]
loss:       14.303407  [  672/ 8082]
loss:       14.553021  [ 1312/ 8082]
loss:       14.179476  [ 1952/ 8082]
loss:       15.772667  [ 2592/ 8082]
loss:       14.487302  [ 3232/ 8082]
loss:       14.998327  [ 3872/ 8082]
loss:       15.898334  [ 4512/ 8082]
loss:       14.597420  [ 5152/ 8082]
loss:       15.214395  [ 5792/ 8082]
loss:       15.202866  [ 6432/ 8082]
loss:       14.289448  [ 7072/ 8082]
loss:       14.730522  [ 7712/ 8082]

Avg training loss:       14.794465
Avg validation loss:                  14.558420
Avg validation habitat loss:           9.030169
Avg validation movement loss:          5.528248
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.569777 --> 14.558420).  Saving model ...


Epoch 15
-------------------------------
loss:       14.339008  [   32/ 8082]
loss:       14.714615  [  672/ 8082]
loss:       15.160459  [ 1312/ 8082]
loss:       14.501738  [ 1952/ 8082]
loss:       14.491870  [ 2592/ 8082]
loss:       15.139768  [ 3232/ 8082]
loss:       14.355402  [ 3872/ 8082]
loss:       14.810845  [ 4512/ 8082]
loss:       15.244129  [ 5152/ 8082]
loss:       13.955277  [ 5792/ 8082]
loss:       14.514570  [ 6432/ 8082]
loss:       14.785206  [ 7072/ 8082]
loss:       14.801582  [ 7712/ 8082]

Avg training loss:       14.790505
Avg validation loss:                  14.545255
Avg validation habitat loss:           9.030570
Avg validation movement loss:          5.514685
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.558420 --> 14.545255).  Saving model ...


Epoch 16
-------------------------------
loss:       15.066435  [   32/ 8082]
loss:       14.762068  [  672/ 8082]
loss:       15.077748  [ 1312/ 8082]
loss:       13.875414  [ 1952/ 8082]
loss:       14.386683  [ 2592/ 8082]
loss:       15.112674  [ 3232/ 8082]
loss:       14.299816  [ 3872/ 8082]
loss:       15.671088  [ 4512/ 8082]
loss:       14.702100  [ 5152/ 8082]
loss:       14.842596  [ 5792/ 8082]
loss:       14.896635  [ 6432/ 8082]
loss:       15.484968  [ 7072/ 8082]
loss:       14.762105  [ 7712/ 8082]

Avg training loss:       14.783350
Avg validation loss:                  14.537501
Avg validation habitat loss:           9.026664
Avg validation movement loss:          5.510835
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.545255 --> 14.537501).  Saving model ...


Epoch 17
-------------------------------
loss:       14.367132  [   32/ 8082]
loss:       14.218849  [  672/ 8082]
loss:       15.487064  [ 1312/ 8082]
loss:       15.100477  [ 1952/ 8082]
loss:       14.079762  [ 2592/ 8082]
loss:       15.423223  [ 3232/ 8082]
loss:       14.470475  [ 3872/ 8082]
loss:       14.977905  [ 4512/ 8082]
loss:       15.139036  [ 5152/ 8082]
loss:       15.647018  [ 5792/ 8082]
loss:       14.377527  [ 6432/ 8082]
loss:       14.794130  [ 7072/ 8082]
loss:       15.221723  [ 7712/ 8082]

Avg training loss:       14.776567
Avg validation loss:                  14.522481
Avg validation habitat loss:           9.027071
Avg validation movement loss:          5.495412
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.537501 --> 14.522481).  Saving model ...


Epoch 18
-------------------------------
loss:       14.945076  [   32/ 8082]
loss:       14.613690  [  672/ 8082]
loss:       15.086289  [ 1312/ 8082]
loss:       13.730463  [ 1952/ 8082]
loss:       13.229572  [ 2592/ 8082]
loss:       14.766150  [ 3232/ 8082]
loss:       14.102248  [ 3872/ 8082]
loss:       15.005037  [ 4512/ 8082]
loss:       15.221640  [ 5152/ 8082]
loss:       15.635056  [ 5792/ 8082]
loss:       14.590537  [ 6432/ 8082]
loss:       14.890905  [ 7072/ 8082]
loss:       14.766616  [ 7712/ 8082]

Avg training loss:       14.776791
Avg validation loss:                  14.514667
Avg validation habitat loss:           9.023906
Avg validation movement loss:          5.490763
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.522481 --> 14.514667).  Saving model ...


Epoch 19
-------------------------------
loss:       14.531549  [   32/ 8082]
loss:       15.593445  [  672/ 8082]
loss:       14.410435  [ 1312/ 8082]
loss:       14.350975  [ 1952/ 8082]
loss:       15.163882  [ 2592/ 8082]
loss:       14.009061  [ 3232/ 8082]
loss:       14.828666  [ 3872/ 8082]
loss:       14.838782  [ 4512/ 8082]
loss:       14.073042  [ 5152/ 8082]
loss:       14.711904  [ 5792/ 8082]
loss:       14.281801  [ 6432/ 8082]
loss:       15.139931  [ 7072/ 8082]
loss:       15.199179  [ 7712/ 8082]

Avg training loss:       14.768194
Avg validation loss:                  14.519581
Avg validation habitat loss:           9.021039
Avg validation movement loss:          5.498543
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 1 out of 10


Epoch 20
-------------------------------
loss:       14.414771  [   32/ 8082]
loss:       14.542617  [  672/ 8082]
loss:       14.483595  [ 1312/ 8082]
loss:       14.956654  [ 1952/ 8082]
loss:       15.102265  [ 2592/ 8082]
loss:       14.041032  [ 3232/ 8082]
loss:       14.944235  [ 3872/ 8082]
loss:       15.292775  [ 4512/ 8082]
loss:       14.952427  [ 5152/ 8082]
loss:       15.059939  [ 5792/ 8082]
loss:       15.519001  [ 6432/ 8082]
loss:       14.759561  [ 7072/ 8082]
loss:       14.893908  [ 7712/ 8082]

Avg training loss:       14.764187
Avg validation loss:                  14.518028
Avg validation habitat loss:           9.020075
Avg validation movement loss:          5.497952
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 2 out of 10


Epoch 21
-------------------------------
loss:       15.020992  [   32/ 8082]
loss:       14.867832  [  672/ 8082]
loss:       14.235175  [ 1312/ 8082]
loss:       13.979012  [ 1952/ 8082]
loss:       15.168735  [ 2592/ 8082]
loss:       14.914450  [ 3232/ 8082]
loss:       15.201292  [ 3872/ 8082]
loss:       14.494457  [ 4512/ 8082]
loss:       14.115108  [ 5152/ 8082]
loss:       14.426221  [ 5792/ 8082]
loss:       14.261929  [ 6432/ 8082]
loss:       14.608278  [ 7072/ 8082]
loss:       15.495301  [ 7712/ 8082]

Avg training loss:       14.765177
Avg validation loss:                  14.502775
Avg validation habitat loss:           9.022175
Avg validation movement loss:          5.480600
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.514667 --> 14.502775).  Saving model ...


Epoch 22
-------------------------------
loss:       15.242526  [   32/ 8082]
loss:       14.291880  [  672/ 8082]
loss:       14.496867  [ 1312/ 8082]
loss:       15.120615  [ 1952/ 8082]
loss:       14.430217  [ 2592/ 8082]
loss:       15.111511  [ 3232/ 8082]
loss:       15.196911  [ 3872/ 8082]
loss:       14.655392  [ 4512/ 8082]
loss:       14.942902  [ 5152/ 8082]
loss:       15.210142  [ 5792/ 8082]
loss:       13.989799  [ 6432/ 8082]
loss:       14.689806  [ 7072/ 8082]
loss:       14.983494  [ 7712/ 8082]

Avg training loss:       14.758531
Avg validation loss:                  14.520660
Avg validation habitat loss:           9.020647
Avg validation movement loss:          5.500016
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 1 out of 10


Epoch 23
-------------------------------
loss:       14.549969  [   32/ 8082]
loss:       15.134369  [  672/ 8082]
loss:       15.411176  [ 1312/ 8082]
loss:       14.960958  [ 1952/ 8082]
loss:       14.136788  [ 2592/ 8082]
loss:       14.820084  [ 3232/ 8082]
loss:       13.372366  [ 3872/ 8082]
loss:       16.235594  [ 4512/ 8082]
loss:       14.364975  [ 5152/ 8082]
loss:       14.612951  [ 5792/ 8082]
loss:       13.739614  [ 6432/ 8082]
loss:       15.559229  [ 7072/ 8082]
loss:       14.697345  [ 7712/ 8082]

Avg training loss:       14.754134
Avg validation loss:                  14.474580
Avg validation habitat loss:           9.020129
Avg validation movement loss:          5.454450
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.502775 --> 14.474580).  Saving model ...


Epoch 24
-------------------------------
loss:       15.041821  [   32/ 8082]
loss:       15.024039  [  672/ 8082]
loss:       14.405666  [ 1312/ 8082]
loss:       15.575415  [ 1952/ 8082]
loss:       14.561767  [ 2592/ 8082]
loss:       15.114403  [ 3232/ 8082]
loss:       14.261971  [ 3872/ 8082]
loss:       14.975683  [ 4512/ 8082]
loss:       15.282489  [ 5152/ 8082]
loss:       15.553199  [ 5792/ 8082]
loss:       14.289170  [ 6432/ 8082]
loss:       15.026646  [ 7072/ 8082]
loss:       15.162102  [ 7712/ 8082]

Avg training loss:       14.747905
Avg validation loss:                  14.503793
Avg validation habitat loss:           9.016137
Avg validation movement loss:          5.487655
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 1 out of 10


Epoch 25
-------------------------------
loss:       14.147286  [   32/ 8082]
loss:       15.219295  [  672/ 8082]
loss:       14.705502  [ 1312/ 8082]
loss:       14.668311  [ 1952/ 8082]
loss:       14.881212  [ 2592/ 8082]
loss:       14.601383  [ 3232/ 8082]
loss:       14.653970  [ 3872/ 8082]
loss:       15.455982  [ 4512/ 8082]
loss:       14.305542  [ 5152/ 8082]
loss:       13.911459  [ 5792/ 8082]
loss:       15.510782  [ 6432/ 8082]
loss:       15.258265  [ 7072/ 8082]
loss:       15.298699  [ 7712/ 8082]

Avg training loss:       14.745745
Avg validation loss:                  14.479233
Avg validation habitat loss:           9.015648
Avg validation movement loss:          5.463584
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 2 out of 10


Epoch 26
-------------------------------
loss:       14.247448  [   32/ 8082]
loss:       15.281538  [  672/ 8082]
loss:       14.364931  [ 1312/ 8082]
loss:       15.609540  [ 1952/ 8082]
loss:       14.456483  [ 2592/ 8082]
loss:       15.138524  [ 3232/ 8082]
loss:       14.576565  [ 3872/ 8082]
loss:       14.660721  [ 4512/ 8082]
loss:       15.451391  [ 5152/ 8082]
loss:       14.425240  [ 5792/ 8082]
loss:       15.281254  [ 6432/ 8082]
loss:       14.729147  [ 7072/ 8082]
loss:       13.930916  [ 7712/ 8082]

Avg training loss:       14.739621
Avg validation loss:                  14.490713
Avg validation habitat loss:           9.015299
Avg validation movement loss:          5.475416
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 3 out of 10


Epoch 27
-------------------------------
loss:       14.470182  [   32/ 8082]
loss:       14.806958  [  672/ 8082]
loss:       13.972561  [ 1312/ 8082]
loss:       15.469616  [ 1952/ 8082]
loss:       14.847608  [ 2592/ 8082]
loss:       14.695177  [ 3232/ 8082]
loss:       14.765018  [ 3872/ 8082]
loss:       15.867525  [ 4512/ 8082]
loss:       14.774474  [ 5152/ 8082]
loss:       14.578945  [ 5792/ 8082]
loss:       14.247944  [ 6432/ 8082]
loss:       15.129126  [ 7072/ 8082]
loss:       15.508838  [ 7712/ 8082]

Avg training loss:       14.737145
Avg validation loss:                  14.484796
Avg validation habitat loss:           9.016390
Avg validation movement loss:          5.468407
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 4 out of 10


Epoch 28
-------------------------------
loss:       15.209536  [   32/ 8082]
loss:       14.969676  [  672/ 8082]
loss:       14.928401  [ 1312/ 8082]
loss:       14.796721  [ 1952/ 8082]
loss:       15.221884  [ 2592/ 8082]
loss:       15.469288  [ 3232/ 8082]
loss:       13.888922  [ 3872/ 8082]
loss:       14.053439  [ 4512/ 8082]
loss:       15.272770  [ 5152/ 8082]
loss:       15.787294  [ 5792/ 8082]
loss:       14.545238  [ 6432/ 8082]
loss:       15.133281  [ 7072/ 8082]
loss:       14.361634  [ 7712/ 8082]

Avg training loss:       14.732975
Avg validation loss:                  14.449441
Avg validation habitat loss:           9.011135
Avg validation movement loss:          5.438306
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.474580 --> 14.449441).  Saving model ...


Epoch 29
-------------------------------
loss:       14.829226  [   32/ 8082]
loss:       14.419657  [  672/ 8082]
loss:       14.388544  [ 1312/ 8082]
loss:       14.849034  [ 1952/ 8082]
loss:       14.471621  [ 2592/ 8082]
loss:       14.765641  [ 3232/ 8082]
loss:       14.684801  [ 3872/ 8082]
loss:       14.300472  [ 4512/ 8082]
loss:       13.888409  [ 5152/ 8082]
loss:       14.762572  [ 5792/ 8082]
loss:       13.533930  [ 6432/ 8082]
loss:       14.697838  [ 7072/ 8082]
loss:       14.800223  [ 7712/ 8082]

Avg training loss:       14.730077
Avg validation loss:                  14.452222
Avg validation habitat loss:           9.012033
Avg validation movement loss:          5.440190
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 1 out of 10


Epoch 30
-------------------------------
loss:       13.908449  [   32/ 8082]
loss:       15.326937  [  672/ 8082]
loss:       14.532779  [ 1312/ 8082]
loss:       14.864893  [ 1952/ 8082]
loss:       14.883301  [ 2592/ 8082]
loss:       14.919694  [ 3232/ 8082]
loss:       15.290260  [ 3872/ 8082]
loss:       14.712601  [ 4512/ 8082]
loss:       15.059449  [ 5152/ 8082]
loss:       14.347562  [ 5792/ 8082]
loss:       16.158821  [ 6432/ 8082]
loss:       14.247931  [ 7072/ 8082]
loss:       14.604774  [ 7712/ 8082]

Avg training loss:       14.728552
Avg validation loss:                  14.451700
Avg validation habitat loss:           9.011653
Avg validation movement loss:          5.440050
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 2 out of 10


Epoch 31
-------------------------------
loss:       14.668903  [   32/ 8082]
loss:       14.946527  [  672/ 8082]
loss:       15.516362  [ 1312/ 8082]
loss:       14.821546  [ 1952/ 8082]
loss:       14.230898  [ 2592/ 8082]
loss:       13.805585  [ 3232/ 8082]
loss:       14.427662  [ 3872/ 8082]
loss:       15.405170  [ 4512/ 8082]
loss:       15.169492  [ 5152/ 8082]
loss:       15.809995  [ 5792/ 8082]
loss:       15.345498  [ 6432/ 8082]
loss:       14.856279  [ 7072/ 8082]
loss:       14.802890  [ 7712/ 8082]

Avg training loss:       14.725499
Avg validation loss:                  14.443869
Avg validation habitat loss:           9.011971
Avg validation movement loss:          5.431900
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.449441 --> 14.443869).  Saving model ...


Epoch 32
-------------------------------
loss:       14.548349  [   32/ 8082]
loss:       14.067467  [  672/ 8082]
loss:       14.208516  [ 1312/ 8082]
loss:       14.624399  [ 1952/ 8082]
loss:       14.405368  [ 2592/ 8082]
loss:       15.069077  [ 3232/ 8082]
loss:       15.562119  [ 3872/ 8082]
loss:       15.290662  [ 4512/ 8082]
loss:       15.286899  [ 5152/ 8082]
loss:       14.281475  [ 5792/ 8082]
loss:       14.719427  [ 6432/ 8082]
loss:       14.990799  [ 7072/ 8082]
loss:       15.227346  [ 7712/ 8082]

Avg training loss:       14.720349
Avg validation loss:                  14.437970
Avg validation habitat loss:           9.009909
Avg validation movement loss:          5.428061
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.443869 --> 14.437970).  Saving model ...


Epoch 33
-------------------------------
loss:       14.499824  [   32/ 8082]
loss:       14.116893  [  672/ 8082]
loss:       14.741138  [ 1312/ 8082]
loss:       15.272849  [ 1952/ 8082]
loss:       14.005964  [ 2592/ 8082]
loss:       14.728837  [ 3232/ 8082]
loss:       14.384792  [ 3872/ 8082]
loss:       14.849600  [ 4512/ 8082]
loss:       14.718244  [ 5152/ 8082]
loss:       14.815147  [ 5792/ 8082]
loss:       14.533769  [ 6432/ 8082]
loss:       14.507512  [ 7072/ 8082]
loss:       14.787693  [ 7712/ 8082]

Avg training loss:       14.718873
Avg validation loss:                  14.449002
Avg validation habitat loss:           9.009474
Avg validation movement loss:          5.439528
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 1 out of 10


Epoch 34
-------------------------------
loss:       14.976305  [   32/ 8082]
loss:       13.903684  [  672/ 8082]
loss:       15.016391  [ 1312/ 8082]
loss:       14.653476  [ 1952/ 8082]
loss:       15.582041  [ 2592/ 8082]
loss:       13.842564  [ 3232/ 8082]
loss:       13.825401  [ 3872/ 8082]
loss:       15.258045  [ 4512/ 8082]
loss:       15.417095  [ 5152/ 8082]
loss:       14.600894  [ 5792/ 8082]
loss:       15.001360  [ 6432/ 8082]
loss:       14.956659  [ 7072/ 8082]
loss:       15.330032  [ 7712/ 8082]

Avg training loss:       14.715694
Avg validation loss:                  14.432675
Avg validation habitat loss:           9.010555
Avg validation movement loss:          5.422120
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.437970 --> 14.432675).  Saving model ...


Epoch 35
-------------------------------
loss:       15.184846  [   32/ 8082]
loss:       14.373926  [  672/ 8082]
loss:       15.017561  [ 1312/ 8082]
loss:       14.607674  [ 1952/ 8082]
loss:       14.934904  [ 2592/ 8082]
loss:       15.325996  [ 3232/ 8082]
loss:       15.412434  [ 3872/ 8082]
loss:       15.109729  [ 4512/ 8082]
loss:       15.222236  [ 5152/ 8082]
loss:       14.526888  [ 5792/ 8082]
loss:       13.992868  [ 6432/ 8082]
loss:       13.761662  [ 7072/ 8082]
loss:       15.358542  [ 7712/ 8082]

Avg training loss:       14.713015
Avg validation loss:                  14.424857
Avg validation habitat loss:           9.010808
Avg validation movement loss:          5.414051
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.432675 --> 14.424857).  Saving model ...


Epoch 36
-------------------------------
loss:       15.107381  [   32/ 8082]
loss:       15.567760  [  672/ 8082]
loss:       14.569506  [ 1312/ 8082]
loss:       14.493985  [ 1952/ 8082]
loss:       15.275295  [ 2592/ 8082]
loss:       14.251236  [ 3232/ 8082]
loss:       14.445346  [ 3872/ 8082]
loss:       14.355404  [ 4512/ 8082]
loss:       14.564994  [ 5152/ 8082]
loss:       14.724042  [ 5792/ 8082]
loss:       14.693325  [ 6432/ 8082]
loss:       14.608762  [ 7072/ 8082]
loss:       15.422507  [ 7712/ 8082]

Avg training loss:       14.712835
Avg validation loss:                  14.433568
Avg validation habitat loss:           9.008803
Avg validation movement loss:          5.424763
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 1 out of 10


Epoch 37
-------------------------------
loss:       14.493797  [   32/ 8082]
loss:       15.038944  [  672/ 8082]
loss:       14.486239  [ 1312/ 8082]
loss:       13.584826  [ 1952/ 8082]
loss:       14.531980  [ 2592/ 8082]
loss:       14.409586  [ 3232/ 8082]
loss:       15.127234  [ 3872/ 8082]
loss:       15.023055  [ 4512/ 8082]
loss:       13.412717  [ 5152/ 8082]
loss:       13.813923  [ 5792/ 8082]
loss:       14.683863  [ 6432/ 8082]
loss:       14.655014  [ 7072/ 8082]
loss:       14.656239  [ 7712/ 8082]

Avg training loss:       14.707429
Avg validation loss:                  14.429677
Avg validation habitat loss:           9.008612
Avg validation movement loss:          5.421065
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 2 out of 10


Epoch 38
-------------------------------
loss:       15.022018  [   32/ 8082]
loss:       14.865845  [  672/ 8082]
loss:       13.732735  [ 1312/ 8082]
loss:       15.153743  [ 1952/ 8082]
loss:       14.401704  [ 2592/ 8082]
loss:       14.544643  [ 3232/ 8082]
loss:       13.837459  [ 3872/ 8082]
loss:       14.870276  [ 4512/ 8082]
loss:       14.304262  [ 5152/ 8082]
loss:       13.935591  [ 5792/ 8082]
loss:       14.454280  [ 6432/ 8082]
loss:       15.163536  [ 7072/ 8082]
loss:       14.173571  [ 7712/ 8082]

Avg training loss:       14.705484
Avg validation loss:                  14.448496
Avg validation habitat loss:           9.008964
Avg validation movement loss:          5.439531
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 3 out of 10


Epoch 39
-------------------------------
loss:       14.323171  [   32/ 8082]
loss:       14.603864  [  672/ 8082]
loss:       14.408138  [ 1312/ 8082]
loss:       14.566893  [ 1952/ 8082]
loss:       15.444710  [ 2592/ 8082]
loss:       15.103184  [ 3232/ 8082]
loss:       14.654699  [ 3872/ 8082]
loss:       14.959540  [ 4512/ 8082]
loss:       14.135837  [ 5152/ 8082]
loss:       15.108889  [ 5792/ 8082]
loss:       15.416900  [ 6432/ 8082]
loss:       14.220594  [ 7072/ 8082]
loss:       14.096067  [ 7712/ 8082]

Avg training loss:       14.705142
Avg validation loss:                  14.433189
Avg validation habitat loss:           9.009150
Avg validation movement loss:          5.424040
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 4 out of 10


Epoch 40
-------------------------------
loss:       14.016753  [   32/ 8082]
loss:       15.408770  [  672/ 8082]
loss:       14.996299  [ 1312/ 8082]
loss:       15.124666  [ 1952/ 8082]
loss:       14.884027  [ 2592/ 8082]
loss:       13.722424  [ 3232/ 8082]
loss:       14.341768  [ 3872/ 8082]
loss:       15.567417  [ 4512/ 8082]
loss:       15.436991  [ 5152/ 8082]
loss:       14.819772  [ 5792/ 8082]
loss:       14.816588  [ 6432/ 8082]
loss:       14.676475  [ 7072/ 8082]
loss:       14.241707  [ 7712/ 8082]

Avg training loss:       14.704233
Avg validation loss:                  14.434471
Avg validation habitat loss:           9.006763
Avg validation movement loss:          5.427708
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 5 out of 10


Epoch 41
-------------------------------
loss:       13.968441  [   32/ 8082]
loss:       14.405289  [  672/ 8082]
loss:       14.243867  [ 1312/ 8082]
loss:       13.721926  [ 1952/ 8082]
loss:       15.151093  [ 2592/ 8082]
loss:       15.179749  [ 3232/ 8082]
loss:       14.514939  [ 3872/ 8082]
loss:       14.724803  [ 4512/ 8082]
loss:       14.318686  [ 5152/ 8082]
loss:       15.513855  [ 5792/ 8082]
loss:       14.318115  [ 6432/ 8082]
loss:       15.300887  [ 7072/ 8082]
loss:       14.495792  [ 7712/ 8082]

Avg training loss:       14.700992
Avg validation loss:                  14.409428
Avg validation habitat loss:           9.006707
Avg validation movement loss:          5.402720
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.424857 --> 14.409428).  Saving model ...


Epoch 42
-------------------------------
loss:       14.046747  [   32/ 8082]
loss:       14.819637  [  672/ 8082]
loss:       13.873738  [ 1312/ 8082]
loss:       15.069615  [ 1952/ 8082]
loss:       14.390927  [ 2592/ 8082]
loss:       16.240936  [ 3232/ 8082]
loss:       14.927973  [ 3872/ 8082]
loss:       14.946957  [ 4512/ 8082]
loss:       13.982300  [ 5152/ 8082]
loss:       15.302126  [ 5792/ 8082]
loss:       15.332224  [ 6432/ 8082]
loss:       14.887763  [ 7072/ 8082]
loss:       14.849896  [ 7712/ 8082]

Avg training loss:       14.700878
Avg validation loss:                  14.405620
Avg validation habitat loss:           9.005770
Avg validation movement loss:          5.399851
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.409428 --> 14.405620).  Saving model ...


Epoch 43
-------------------------------
loss:       15.070433  [   32/ 8082]
loss:       14.142093  [  672/ 8082]
loss:       14.900124  [ 1312/ 8082]
loss:       15.968881  [ 1952/ 8082]
loss:       15.219810  [ 2592/ 8082]
loss:       14.619028  [ 3232/ 8082]
loss:       15.164410  [ 3872/ 8082]
loss:       14.023293  [ 4512/ 8082]
loss:       14.562815  [ 5152/ 8082]
loss:       14.083849  [ 5792/ 8082]
loss:       14.875497  [ 6432/ 8082]
loss:       15.584556  [ 7072/ 8082]
loss:       15.643020  [ 7712/ 8082]

Avg training loss:       14.694345
Avg validation loss:                  14.405678
Avg validation habitat loss:           9.005078
Avg validation movement loss:          5.400600
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 1 out of 10


Epoch 44
-------------------------------
loss:       15.094213  [   32/ 8082]
loss:       13.260601  [  672/ 8082]
loss:       14.425416  [ 1312/ 8082]
loss:       14.198008  [ 1952/ 8082]
loss:       14.927270  [ 2592/ 8082]
loss:       14.514199  [ 3232/ 8082]
loss:       15.044137  [ 3872/ 8082]
loss:       15.080734  [ 4512/ 8082]
loss:       13.810625  [ 5152/ 8082]
loss:       14.707187  [ 5792/ 8082]
loss:       15.045238  [ 6432/ 8082]
loss:       13.402869  [ 7072/ 8082]
loss:       14.630686  [ 7712/ 8082]

Avg training loss:       14.691822
Avg validation loss:                  14.392460
Avg validation habitat loss:           9.007193
Avg validation movement loss:          5.385266
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.405620 --> 14.392460).  Saving model ...


Epoch 45
-------------------------------
loss:       14.630781  [   32/ 8082]
loss:       14.565153  [  672/ 8082]
loss:       14.777569  [ 1312/ 8082]
loss:       14.681658  [ 1952/ 8082]
loss:       15.340382  [ 2592/ 8082]
loss:       15.352612  [ 3232/ 8082]
loss:       14.845429  [ 3872/ 8082]
loss:       14.920072  [ 4512/ 8082]
loss:       14.709377  [ 5152/ 8082]
loss:       15.150167  [ 5792/ 8082]
loss:       14.820426  [ 6432/ 8082]
loss:       14.042624  [ 7072/ 8082]
loss:       15.076562  [ 7712/ 8082]

Avg training loss:       14.694293
Avg validation loss:                  14.397444
Avg validation habitat loss:           9.003758
Avg validation movement loss:          5.393686
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 1 out of 10


Epoch 46
-------------------------------
loss:       14.692982  [   32/ 8082]
loss:       14.622877  [  672/ 8082]
loss:       14.966135  [ 1312/ 8082]
loss:       14.052504  [ 1952/ 8082]
loss:       15.286769  [ 2592/ 8082]
loss:       14.486897  [ 3232/ 8082]
loss:       15.183427  [ 3872/ 8082]
loss:       13.568325  [ 4512/ 8082]
loss:       13.975044  [ 5152/ 8082]
loss:       14.616806  [ 5792/ 8082]
loss:       14.213511  [ 6432/ 8082]
loss:       15.171218  [ 7072/ 8082]
loss:       14.711535  [ 7712/ 8082]

Avg training loss:       14.688455
Avg validation loss:                  14.423701
Avg validation habitat loss:           9.007339
Avg validation movement loss:          5.416360
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 2 out of 10


Epoch 47
-------------------------------
loss:       14.116798  [   32/ 8082]
loss:       14.421505  [  672/ 8082]
loss:       15.508658  [ 1312/ 8082]
loss:       15.143512  [ 1952/ 8082]
loss:       14.312105  [ 2592/ 8082]
loss:       14.464258  [ 3232/ 8082]
loss:       14.191757  [ 3872/ 8082]
loss:       15.111117  [ 4512/ 8082]
loss:       14.239563  [ 5152/ 8082]
loss:       14.343143  [ 5792/ 8082]
loss:       14.053215  [ 6432/ 8082]
loss:       15.614652  [ 7072/ 8082]
loss:       14.306411  [ 7712/ 8082]

Avg training loss:       14.688252
Avg validation loss:                  14.420854
Avg validation habitat loss:           9.008953
Avg validation movement loss:          5.411902
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 3 out of 10


Epoch 48
-------------------------------
loss:       14.614623  [   32/ 8082]
loss:       14.554979  [  672/ 8082]
loss:       15.073971  [ 1312/ 8082]
loss:       15.290283  [ 1952/ 8082]
loss:       14.254403  [ 2592/ 8082]
loss:       14.011925  [ 3232/ 8082]
loss:       14.790536  [ 3872/ 8082]
loss:       13.779184  [ 4512/ 8082]
loss:       13.933920  [ 5152/ 8082]
loss:       14.944463  [ 5792/ 8082]
loss:       14.982913  [ 6432/ 8082]
loss:       14.457654  [ 7072/ 8082]
loss:       14.047795  [ 7712/ 8082]

Avg training loss:       14.681406
Avg validation loss:                  14.391467
Avg validation habitat loss:           9.005760
Avg validation movement loss:          5.385708
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
Validation loss decreased (14.392460 --> 14.391467).  Saving model ...


Epoch 49
-------------------------------
loss:       15.563956  [   32/ 8082]
loss:       14.409714  [  672/ 8082]
loss:       15.059926  [ 1312/ 8082]
loss:       14.641969  [ 1952/ 8082]
loss:       14.362056  [ 2592/ 8082]
loss:       15.436552  [ 3232/ 8082]
loss:       14.979259  [ 3872/ 8082]
loss:       14.065618  [ 4512/ 8082]
loss:       12.887695  [ 5152/ 8082]
loss:       15.610971  [ 5792/ 8082]
loss:       13.715605  [ 6432/ 8082]
loss:       14.608070  [ 7072/ 8082]
loss:       14.326908  [ 7712/ 8082]

Avg training loss:       14.680515
Avg validation loss:                  14.424132
Avg validation habitat loss:           9.008206
Avg validation movement loss:          5.415929
Movement learning rate:         [1e-05]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 1 out of 10


Epoch 50
-------------------------------
loss:       14.479509  [   32/ 8082]
loss:       14.407170  [  672/ 8082]
loss:       14.838150  [ 1312/ 8082]
loss:       15.157505  [ 1952/ 8082]
loss:       15.777323  [ 2592/ 8082]
loss:       14.642736  [ 3232/ 8082]
loss:       14.021012  [ 3872/ 8082]
loss:       15.064647  [ 4512/ 8082]
loss:       14.878405  [ 5152/ 8082]
loss:       13.949443  [ 5792/ 8082]
loss:       14.795333  [ 6432/ 8082]
loss:       14.874767  [ 7072/ 8082]
loss:       15.023170  [ 7712/ 8082]

Avg training loss:       14.677996
Avg validation loss:                  14.418661
Avg validation habitat loss:           9.006645
Avg validation movement loss:          5.412018
Movement learning rate:         [1.0000000000000002e-06]
Habitat learning rate:          [0.0001]
EarlyStopping counter: 2 out of 10


Epoch 51
-------------------------------
loss:       13.860001  [   32/ 8082]
loss:       15.233062  [  672/ 8082]
loss:       15.113272  [ 1312/ 8082]
loss:       14.674922  [ 1952/ 8082]
loss:       15.725531  [ 2592/ 8082]
loss:       14.455893  [ 3232/ 8082]
loss:       14.780814  [ 3872/ 8082]
loss:       14.623874  [ 4512/ 8082]
loss:       15.325014  [ 5152/ 8082]
loss:       14.545890  [ 5792/ 8082]
loss:       14.104012  [ 6432/ 8082]
loss:       15.049157  [ 7072/ 8082]
loss:       15.398517  [ 7712/ 8082]

Avg training loss:       14.678690
Avg validation loss:                  14.403462
Avg validation habitat loss:           9.018920
Avg validation movement loss:          5.384542
Movement learning rate:         [1.0000000000000002e-06]
Habitat learning rate:          [1e-05]
EarlyStopping counter: 3 out of 10


Epoch 52
-------------------------------
loss:       14.788683  [   32/ 8082]
loss:       13.488342  [  672/ 8082]
loss:       14.435366  [ 1312/ 8082]
loss:       13.664868  [ 1952/ 8082]
loss:       15.421284  [ 2592/ 8082]
loss:       14.140110  [ 3232/ 8082]
loss:       14.733646  [ 3872/ 8082]
loss:       14.113314  [ 4512/ 8082]
loss:       14.357676  [ 5152/ 8082]
loss:       14.553676  [ 5792/ 8082]
loss:       13.481390  [ 6432/ 8082]
loss:       15.013284  [ 7072/ 8082]
loss:       14.898931  [ 7712/ 8082]

Avg training loss:       14.673674
Avg validation loss:                  14.406239
Avg validation habitat loss:           9.012079
Avg validation movement loss:          5.394159
Movement learning rate:         [1.0000000000000002e-06]
Habitat learning rate:          [1e-05]
EarlyStopping counter: 4 out of 10


Epoch 53
-------------------------------
loss:       13.768931  [   32/ 8082]
loss:       14.353657  [  672/ 8082]
loss:       14.123418  [ 1312/ 8082]
loss:       15.040705  [ 1952/ 8082]
loss:       14.722520  [ 2592/ 8082]
loss:       14.156158  [ 3232/ 8082]
loss:       13.541374  [ 3872/ 8082]
loss:       14.467157  [ 4512/ 8082]
loss:       14.682388  [ 5152/ 8082]
loss:       15.036253  [ 5792/ 8082]
loss:       14.503506  [ 6432/ 8082]
loss:       14.747510  [ 7072/ 8082]
loss:       14.135642  [ 7712/ 8082]

Avg training loss:       14.677609
Avg validation loss:                  14.417871
Avg validation habitat loss:           9.018906
Avg validation movement loss:          5.398966
Movement learning rate:         [1.0000000000000002e-06]
Habitat learning rate:          [1e-05]
EarlyStopping counter: 5 out of 10


Epoch 54
-------------------------------
loss:       15.290676  [   32/ 8082]
loss:       15.095904  [  672/ 8082]
loss:       14.388249  [ 1312/ 8082]
loss:       15.103960  [ 1952/ 8082]
loss:       14.331125  [ 2592/ 8082]
loss:       14.340302  [ 3232/ 8082]
loss:       15.032572  [ 3872/ 8082]
loss:       13.641793  [ 4512/ 8082]
loss:       14.513950  [ 5152/ 8082]
loss:       14.876392  [ 5792/ 8082]
loss:       15.003579  [ 6432/ 8082]
loss:       14.640949  [ 7072/ 8082]
loss:       15.066162  [ 7712/ 8082]

Avg training loss:       14.676331
Avg validation loss:                  14.404132
Avg validation habitat loss:           9.008471
Avg validation movement loss:          5.395663
Movement learning rate:         [1.0000000000000002e-06]
Habitat learning rate:          [1e-05]
EarlyStopping counter: 6 out of 10


Epoch 55
-------------------------------
loss:       13.942184  [   32/ 8082]
loss:       13.688111  [  672/ 8082]
loss:       16.257336  [ 1312/ 8082]
loss:       15.309095  [ 1952/ 8082]
loss:       16.012276  [ 2592/ 8082]
loss:       14.815039  [ 3232/ 8082]
loss:       14.440355  [ 3872/ 8082]
loss:       14.220662  [ 4512/ 8082]
loss:       13.997902  [ 5152/ 8082]
loss:       14.712907  [ 5792/ 8082]
loss:       14.384615  [ 6432/ 8082]
loss:       14.257649  [ 7072/ 8082]
loss:       15.522534  [ 7712/ 8082]

Avg training loss:       14.674029
Avg validation loss:                  14.407248
Avg validation habitat loss:           9.008042
Avg validation movement loss:          5.399207
Movement learning rate:         [1.0000000000000002e-06]
Habitat learning rate:          [1e-05]
EarlyStopping counter: 7 out of 10


Epoch 56
-------------------------------
loss:       14.967100  [   32/ 8082]
loss:       15.162800  [  672/ 8082]
loss:       14.741832  [ 1312/ 8082]
loss:       15.290051  [ 1952/ 8082]
loss:       14.549616  [ 2592/ 8082]
loss:       14.567600  [ 3232/ 8082]
loss:       15.023174  [ 3872/ 8082]
loss:       14.756496  [ 4512/ 8082]
loss:       14.786263  [ 5152/ 8082]
loss:       14.864202  [ 5792/ 8082]
loss:       14.881171  [ 6432/ 8082]
loss:       15.186590  [ 7072/ 8082]
loss:       14.744934  [ 7712/ 8082]

Avg training loss:       14.674669
Avg validation loss:                  14.410643
Avg validation habitat loss:           9.009926
Avg validation movement loss:          5.400716
Movement learning rate:         [1.0000000000000002e-06]
Habitat learning rate:          [1e-05]
EarlyStopping counter: 8 out of 10


Epoch 57
-------------------------------
loss:       14.177410  [   32/ 8082]
loss:       14.846877  [  672/ 8082]
loss:       15.258842  [ 1312/ 8082]
loss:       15.355652  [ 1952/ 8082]
loss:       14.920850  [ 2592/ 8082]
loss:       15.257000  [ 3232/ 8082]
loss:       14.099021  [ 3872/ 8082]
loss:       14.435457  [ 4512/ 8082]
loss:       14.708059  [ 5152/ 8082]
loss:       15.555657  [ 5792/ 8082]
loss:       14.424987  [ 6432/ 8082]
loss:       14.010809  [ 7072/ 8082]
loss:       14.267629  [ 7712/ 8082]

Avg training loss:       14.677167
Avg validation loss:                  14.401154
Avg validation habitat loss:           9.008599
Avg validation movement loss:          5.392555
Movement learning rate:         [1.0000000000000002e-07]
Habitat learning rate:          [1.0000000000000002e-06]
EarlyStopping counter: 9 out of 10


Epoch 58
-------------------------------
loss:       15.318058  [   32/ 8082]
loss:       14.540492  [  672/ 8082]
loss:       15.066505  [ 1312/ 8082]
loss:       14.143305  [ 1952/ 8082]
loss:       14.668678  [ 2592/ 8082]
loss:       14.069659  [ 3232/ 8082]
loss:       15.123915  [ 3872/ 8082]
loss:       15.042867  [ 4512/ 8082]
loss:       13.843367  [ 5152/ 8082]
loss:       15.329596  [ 5792/ 8082]
loss:       13.813165  [ 6432/ 8082]
loss:       14.756277  [ 7072/ 8082]
loss:       14.579089  [ 7712/ 8082]

Avg training loss:       14.676712
Avg validation loss:                  14.400331
Avg validation habitat loss:           9.008992
Avg validation movement loss:          5.391338
Movement learning rate:         [1.0000000000000002e-07]
Habitat learning rate:          [1.0000000000000002e-06]
EarlyStopping counter: 10 out of 10
Early stopping
Avg test loss:          14.327008 

Done!

Make a GIF of the training images

First, here’s a function to call to make a gif from a given directory.

# Example sorting by the epoch number
def extract_epoch(filename):
    # Extract the epoch number from the filename
    # Adjust the extraction based on your naming pattern
    import re
    match = re.search(r'training_epoch(\d+)_', filename)
    if match:
        return int(match.group(1))
    return 0

def create_gif(image_folder, output_filename, fps=10):
    """
    Creates a GIF from a sequence of images in a folder.

    Parameters:
    - image_folder: Path to the folder containing images
    - output_filename: Name of the output GIF file
    - duration: Duration of each frame in seconds
    """
    # Get all png files in the specified folder, sorted by name
    images = sorted(glob.glob(os.path.join(image_folder, '*.png')), key=extract_epoch)

    # Check if any images were found
    if not images:
        print(f"No images found in {image_folder}")
        return

    # Read all images
    frames = [imageio.imread(image) for image in images]

    # Save as GIF
    imageio.mimsave(output_filename, frames, fps=fps, loop=0)

    display(Image(filename=output_filename))

    print(f"GIF created successfully: {output_filename}")

# Path to your images
image_folder =  f'{output_dir}/training_images'

# Output GIF filename
output_filename = f'{output_dir}/training_gif_id{buffalo_id}_yday{yday_t2_integer}_hour{hour_t2_integer}_bearing{bearing_degrees}_next_r{row}_c{column}.gif'

# Create the GIF
create_gif(image_folder, output_filename, fps=10)

<IPython.core.display.Image object>

GIF created successfully: ../Python/outputs/model_training/id2005_deepSSF_training_1_2025-07-09/training_gif_id2005_yday324_hour17_bearing156_next_r46_c44.gif

# to look at the parameters (weights and biases) of the model
# print(model.state_dict())

Loading in previous models

As we’ve trained the model, the model parameters are already stored in the model object. But as we were training the model, we were saving it to file, and that, and other trained models can be loaded.

The model parameters that are being loaded must match the model object that has been defined above. If the model object has changed, the model parameters will not be able to be loaded.

path_save_weights

'../Python/outputs/model_training/id2005_deepSSF_training_1_2025-07-09/checkpoint_deepSSF_buffalo2005.pt'

If loading a previously trained model

# to load previously saved weights
# path_save_weights = f'{output_dir}/checkpoint_deepSSF_buffalo2005_2025-04-01.pt'

model.load_state_dict(torch.load(path_save_weights,
                                 weights_only=True,
                                 map_location=torch.device('cpu')))

<All keys matched successfully>

View model outputs

Create a directory to save model outputs

Save the validation loss as a dataframe

# Directory for saving the loss dataframe
filename_loss_csv = f'{output_dir}/deepSSF_val_loss_buffalo{buffalo_id}.csv'

# Check if val_losses is defined (which means a model has been trained in this session)
try:

    # Convert the list of tensors to a single tensor
    val_losses_tensor = torch.tensor(val_losses)

    print("val_losses has been defined - storing as csv\n")

    # Number of epochs
    n_epochs = len(val_losses)
    print(f'Number of epochs: {n_epochs}')

    val_losses_df = pd.DataFrame({
        "epoch": range(1, n_epochs + 1),
        "val_losses": val_losses_tensor.detach().cpu().numpy()
    })

    print(val_losses_df.head())

    # Save the validation losses to a CSV file
    val_losses_df.to_csv(filename_loss_csv, index=False)

# if val_losses hasn't been defined (for if you are loading model weights from a saved object)
except NameError:

    # This code runs if val_losses is not defined
    print("val_losses has not been defined - loading from saved csv\n")
    # Initialize it with a default value

    # Read the val_losses csv file
    val_losses_df = pd.read_csv(filename_loss_csv)
    print(val_losses_df.head())

    # Number of epochs
    n_epochs = len(val_losses_df)
    print(f'\nNumber of epochs: {n_epochs}')

val_losses has been defined - storing as csv

Number of epochs: 58
   epoch  val_losses
0      1   15.368755
1      2   15.255189
2      3   15.095748
3      4   14.955804
4      5   14.828360

# Directory for saving the loss dataframe
filename_train_loss_csv = f'{output_dir}/deepSSF_train_loss_buffalo{buffalo_id}.csv'

# Check if train_losses is defined (which means a model has been trained in this session)
try:

    # Convert the list of tensors to a single tensor
    train_losses_tensor = torch.tensor(train_losses)

    print("train_losses has been defined - storing as csv\n")

    train_losses_df = pd.DataFrame({
        "epoch": np.linspace(1, n_epochs, len(train_losses)),
        "train_losses": train_losses_tensor.detach().cpu().numpy()
    })

    print(train_losses_df.head)

    # Save the train losses to a CSV file
    train_losses_df.to_csv(filename_train_loss_csv, index=False)

# if train_losses hasn't been defined (for if you are loading model weights from a saved object)
except NameError:

    # This code runs if train_losses is not defined
    print("train_losses has not been defined - loading from saved csv\n")
    # Initialize it with a default value

    # Read the train_losses csv file
    train_losses_df = pd.read_csv(filename_train_loss_csv)
    print(train_losses_df.head())

train_losses has been defined - storing as csv

<bound method NDFrame.head of     epoch  train_losses
0     1.0     15.575754
1     2.0     15.514331
2     3.0     15.329888
3     4.0     15.165933
4     5.0     15.027689
5     6.0     14.955096
6     7.0     14.911941
7     8.0     14.880715
8     9.0     14.854919
9    10.0     14.838990
10   11.0     14.823997
11   12.0     14.813732
12   13.0     14.805429
13   14.0     14.794465
14   15.0     14.790505
15   16.0     14.783350
16   17.0     14.776567
17   18.0     14.776791
18   19.0     14.768194
19   20.0     14.764187
20   21.0     14.765177
21   22.0     14.758531
22   23.0     14.754134
23   24.0     14.747905
24   25.0     14.745745
25   26.0     14.739621
26   27.0     14.737145
27   28.0     14.732975
28   29.0     14.730077
29   30.0     14.728552
30   31.0     14.725499
31   32.0     14.720349
32   33.0     14.718873
33   34.0     14.715694
34   35.0     14.713015
35   36.0     14.712835
36   37.0     14.707429
37   38.0     14.705484
38   39.0     14.705142
39   40.0     14.704233
40   41.0     14.700992
41   42.0     14.700878
42   43.0     14.694345
43   44.0     14.691822
44   45.0     14.694293
45   46.0     14.688455
46   47.0     14.688252
47   48.0     14.681406
48   49.0     14.680515
49   50.0     14.677996
50   51.0     14.678690
51   52.0     14.673674
52   53.0     14.677609
53   54.0     14.676331
54   55.0     14.674029
55   56.0     14.674669
56   57.0     14.677167
57   58.0     14.676712>

Plot the validation loss

# Directory for saving the loss plots
filename_loss = f'{output_dir}/val_loss_buffalo{buffalo_id}.png'

# Plot the validation losses
plt.plot(train_losses_df['epoch'], train_losses_df['train_losses'], label='Training Loss', color='blue')  # Plot training loss in blue
plt.plot(val_losses_df['epoch'], val_losses_df['val_losses'], label='Validation Loss', color='red')  # Plot validation loss in red
plt.title('Validation Losses')
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.legend()  # Show legend to distinguish lines
plt.savefig(filename_loss, dpi=300, bbox_inches='tight')
plt.show()

Test the model on sample covariates

Use the same sample covariates as above (selected before training the model) from the test dataset, and then run the model on these covariates. To select a different sample, change the iteration_index in the code above and re-run those chunks.

Print the scalar values and plot the sample covariates

# Print relevant information about the current prediction context
# such as time of day, day of year, and bearing angles in both radians and degrees.
print(f'Hour:               {hour_t2_integer}')
print(f'Day of the year:    {yday_t2_integer}')
print(f'Bearing (radians):  {bearing}')
print(f'Bearing (degrees):  {bearing_degrees}')

# Plot the covariates
fig, axs = plt.subplots(2, 2, figsize=(9, 7.5))

# Plot NDVI
im1 = axs[0, 0].imshow(ndvi_natural.numpy(), cmap='viridis')
axs[0, 0].set_title('NDVI')
fig.colorbar(im1, ax=axs[0, 0])

# Plot Canopy cover
im2 = axs[0, 1].imshow(canopy_natural.numpy(), cmap='viridis')
axs[0, 1].set_title('Canopy cover')
fig.colorbar(im2, ax=axs[0, 1])

# Plot Herbaceous vegetation
im3 = axs[1, 0].imshow(herby_natural.numpy(), cmap='viridis')
axs[1, 0].set_title('Herbaceous vegetation')
fig.colorbar(im3, ax=axs[1, 0])

# Plot Slope
im4 = axs[1, 1].imshow(slope_natural.numpy(), cmap='viridis')
axs[1, 1].set_title('Slope')
fig.colorbar(im4, ax=axs[1, 1])

filename_covs = f'{output_dir}/id{buffalo_id}_sample{iteration_index}_yday{yday_t2_integer}_hour{hour_t2_integer}_bearing{bearing_degrees}_next_r{row}_c{column}.png'
plt.tight_layout()
plt.savefig(filename_covs, dpi=300, bbox_inches='tight') # if we want to save the figure
plt.show()
plt.close()  # Close the figure to free memory

Hour:               17
Day of the year:    324
Bearing (radians):  2.7272613048553467
Bearing (degrees):  156

Run the model on the sample covariates

# -------------------------------------------------------------------------
# Switch the model to evaluation mode (e.g., disables dropout, etc.)
# -------------------------------------------------------------------------
model.eval()

# -------------------------------------------------------------------------
# Pass the inputs through the model; 'test' will have shape [batch, H, W, 2]
# -------------------------------------------------------------------------
test = model((sample_spatial_covs, sample_temporal_covs, sample_prev_bearing))
# test = model((x1, x2, x3))
print(test.shape)

# -------------------------------------------------------------------------
# Extract and exponentiate the habitat density channel
#    (at index 0 in the last dimension)
# -------------------------------------------------------------------------
hab_density = test.detach().cpu().numpy()[0, :, :, 0]
hab_density_exp = np.exp(hab_density)
# print(np.sum(hab_density_exp))  # Debug: check the sum of exponentiated values

# -------------------------------------------------------------------------
# Create masks to remove unwanted edge cells from visualization
#    (setting them to -∞ affects the color scale in plots)
# -------------------------------------------------------------------------
x_mask = np.ones_like(hab_density)
y_mask = np.ones_like(hab_density)

x_mask[:, :3] = -np.inf
x_mask[:, 98:] = -np.inf
y_mask[:3, :] = -np.inf
y_mask[98:, :] = -np.inf

# -------------------------------------------------------------------------
# Apply the masks to the habitat density (log scale) and exponentiated version
# -------------------------------------------------------------------------
hab_density_mask = hab_density * x_mask * y_mask
hab_density_exp_mask = hab_density_exp * x_mask * y_mask

# -------------------------------------------------------------------------
# Plot and save the habitat density in log scale
# -------------------------------------------------------------------------
plt.imshow(hab_density_mask)
plt.colorbar()
plt.title('Habitat selection probability (log)')
plt.savefig(f'{output_dir}/hab_log_prob_id{buffalo_id}.png', dpi=300, bbox_inches='tight')
plt.show()
plt.close()

# -------------------------------------------------------------------------
# Plot and save the habitat density in probability (exponentiated) scale
# -------------------------------------------------------------------------
plt.imshow(hab_density_exp_mask)
plt.colorbar()
plt.title('Habitat selection probability')
plt.savefig(f'{output_dir}/hab_prob_id{buffalo_id}.png', dpi=300, bbox_inches='tight')
plt.show()
plt.close()

# -------------------------------------------------------------------------
# Extract and exponentiate the movement density channel
#    (at index 1 in the last dimension)
# -------------------------------------------------------------------------
move_density = test.detach().cpu().numpy()[0, :, :, 1]
move_density_exp = np.exp(move_density)
# print(np.sum(move_density_exp))  # Debug: check the sum of exponentiated values

# -------------------------------------------------------------------------
# Apply the same masking strategy to movement densities
# -------------------------------------------------------------------------
move_density_mask = move_density * x_mask * y_mask
move_density_exp_mask = move_density_exp * x_mask * y_mask

# -------------------------------------------------------------------------
# Plot and save the movement density in log scale
# -------------------------------------------------------------------------
plt.imshow(move_density_mask)
plt.colorbar()
plt.title('Movement probability (log)')
plt.savefig(f'{output_dir}/move_log_prob_id{buffalo_id}.png', dpi=300, bbox_inches='tight')
plt.show()
plt.close()

# -------------------------------------------------------------------------
# Plot and save the movement density in probability (exponentiated) scale
# -------------------------------------------------------------------------
plt.imshow(move_density_exp_mask)
plt.colorbar()
plt.title('Movement probability')
plt.savefig(f'{output_dir}/move_prob_id{buffalo_id}.png', dpi=300, bbox_inches='tight')
plt.show()
plt.close()

# -------------------------------------------------------------------------
# Compute the next-step density by adding habitat + movement (log-space)
#     Then exponentiate and normalize
# -------------------------------------------------------------------------
step_density = test[0, :, :, 0] + test[0, :, :, 1]
step_density = step_density.detach().cpu().numpy()
step_density_exp = np.exp(step_density)
# print('Sum of step density exp = ', np.sum(step_density_exp))  # Debug

step_density_exp_norm = step_density_exp / np.sum(step_density_exp)
# print('Sum of step density exp norm = ', np.sum(step_density_exp_norm))  # Debug

# -------------------------------------------------------------------------
# Apply masks to the step densities (log and exponentiated + normalized)
# -------------------------------------------------------------------------
step_density_mask = step_density * x_mask * y_mask
step_density_exp_norm_mask = step_density_exp_norm * x_mask * y_mask

# -------------------------------------------------------------------------
# Plot and save the combined next-step probability surface in log scale
# -------------------------------------------------------------------------
plt.imshow(step_density_mask)
plt.colorbar()
plt.title('Next-step probability (log)')
plt.savefig(f'{output_dir}/step_log_prob_id{buffalo_id}.png', dpi=300, bbox_inches='tight')
plt.show()
plt.close()

# -------------------------------------------------------------------------
# Plot and save the combined next-step probability surface in probability scale
# -------------------------------------------------------------------------
plt.imshow(step_density_exp_norm_mask)
plt.colorbar()
plt.title('Next-step probability')
plt.savefig(f'{output_dir}/step_prob_id{buffalo_id}.png', dpi=300, bbox_inches='tight')
plt.show()
plt.close()

torch.Size([1, 101, 101, 2])

Extracting convolution layer outputs

In the convolutional blocks, each convolutional layer learns a set of filters (kernels) that extract different features from the input data. In the habitat selection subnetwork, the convolution filters (and their associated bias parameters - not shown below) are the only parameters that are trained, and it is the filters that transform the set of input covariates into the habitat selection probabilities. They do this by maximising features of the inputs that correlate with observed next-steps.

For each convolutional layer, there are typically a number of filters. For the habitat selection subnetwork, we used 4 filters in the first two layers, and a single filter in the last layer. Each of these filters has a number of channels which correspond one-to-one with the input layers. The outputs of the filter channels are then combined to produce a feature map, with a single feature map produced for each filter. In successive layers, the feature maps become the input layers, and the filters operate on these layers. Because there are multiple filters in ech layer, they can ‘specialise’ in extracting different features from the input layers.

By visualizing and inspecting these filters, and the corresponding feature maps, we can:

• Gain interpretability: Understand what kind of features the network is detecting—e.g., edges, shapes, or textures.
• Debug: Check if the filters have meaningful patterns or if something went wrong (e.g., all zeros or random noise).
• Compare layers: See how early layers often learn low-level patterns while deeper layers learn more abstract features.

We will first set up some activation hooks for storing the feature maps. Activation hooks are placed at certain points within the model’s forward pass and store intermediate results. We will also extract the convolution filters (which are weights of the model and as such don’t require hooks - we can access them directly).

We will then run the sample covariates through the model and extract the feature maps from the habitat selection convolutional block, and plot them along with the covariates and convolution filters.

Note that there are also ReLU activation functions in the convolutional blocks, which are not shown below. These are applied to the feature maps, and set all negative values to zero. They are not learned parameters, but are part of the forward pass of the model.

Convolutional layer 1

Activation hook

# -----------------------------------------------------------
# Create a dictionary to store activation outputs
# -----------------------------------------------------------
activation = {}

def get_activation(name):
    """
    Returns a hook function that can be registered on a layer
    to capture its output (i.e., feature maps) after the forward pass.

    Args:
        name (str): The key under which the activation is stored in the 'activation' dict.
    """
    def hook(model, input, output):
        # Detach and save the layer's output in the dictionary
        activation[name] = output.detach()
    return hook

# -----------------------------------------------------------
# Register a forward hook on the first convolution layer
#    in the model's 'conv_habitat' block
# -----------------------------------------------------------
model.conv_habitat.conv2d[0].register_forward_hook(get_activation("hab_conv1"))

# -----------------------------------------------------------
# Perform a forward pass through the model with the desired input
#    The feature maps from the hooked layer will be stored in 'activation'
# -----------------------------------------------------------
# out = model((x1, x2, x3))  # e.g., model((spatial_data_x, scalars_to_grid, bearing_x))
test = model((sample_spatial_covs, sample_temporal_covs, sample_prev_bearing))

# -----------------------------------------------------------
# Retrieve the captured feature maps from the dictionary
#    and move them to the CPU for inspection
# -----------------------------------------------------------
feat_maps1 = activation["hab_conv1"].cpu()
print("Feature map shape:", feat_maps1.shape)
# Typically shape: (batch_size, out_channels, height, width)

# -----------------------------------------------------------
# Visualize the feature maps for the first sample in the batch
# -----------------------------------------------------------
feat_maps1_sample = feat_maps1[0]  # Shape: (out_channels, H, W)
num_maps1 = feat_maps1_sample.shape[0]
print("Number of feature maps:", num_maps1)

Feature map shape: torch.Size([1, 4, 101, 101])
Number of feature maps: 4

Stack spatial and scalar (as grid) covariates

For plotting. Also create a vector of names to index over.

covariate_stack = torch.cat([sample_spatial_covs.detach().cpu(), scalar_maps], dim=1)
print(covariate_stack.shape)

covariate_names = ['NDVI', 'Canopy cover', 'Herbaceous vegetation', 'Slope',
                   'Hour sin', 'Hour cos', 'yday sin', 'yday cos']

torch.Size([1, 8, 101, 101])

Extract filters and plot

# -------------------------------------------------------------------------
# Check or print the convolution layer in conv_habitat (for debugging)
# -------------------------------------------------------------------------
print(model.conv_habitat.conv2d)

# -------------------------------------------------------------------------
# Set the model to evaluation mode (disables dropout, etc.)
# -------------------------------------------------------------------------
model.eval()

# -------------------------------------------------------------------------
# Extract the weights (filters) from the first convolution layer in conv_habitat
# -------------------------------------------------------------------------
filters_c1 = model.conv_habitat.conv2d[0].weight.data.clone().cpu()
print("Filters shape:", filters_c1.shape)
# Typically (out_channels, in_channels, kernel_height, kernel_width)

# -------------------------------------------------------------------------
# Visualize each filter’s first channel in a grid of subplots
# -------------------------------------------------------------------------
num_filters_c1 = filters_c1.shape[1]
print(num_filters_c1)

for z in range(num_maps1):

    fig, axes = plt.subplots(2, num_filters_c1, figsize=(2*num_filters_c1, 4))
    for i in range(num_filters_c1):

        # Add the covariates as the first row of subplots
        axes[0,i].imshow(covariate_stack[0, i].detach().cpu().numpy(), cmap='viridis')
        axes[0,i].axis('off')
        axes[0,i].set_title(f'{covariate_names[i]}')
        if i > 3:
            im1 = axes[0,i].imshow(covariate_stack[0, i].detach().cpu().numpy(), cmap='viridis')
            im1.set_clim(-1, 1)
            axes[0,i].text(scalar_maps.shape[2] // 2, scalar_maps.shape[3] // 2,
                f'Value: {round(sample_temporal_covs[0, i-4].item(), 2)}',
                ha='center', va='center', color='white', fontsize=12)

        kernel = filters_c1[z, i, :, :]  # Show the first input channel
        im = axes[1,i].imshow(kernel, cmap='viridis')
        axes[1,i].axis('off')
        axes[1,i].set_title(f'Layer 1, Filter {z+1}')
        # Annotate each cell with the numeric value
        for (j, k), val in np.ndenumerate(kernel):
            axes[1,i].text(k, j, f'{val:.2f}', ha='center', va='center', color='white')

    plt.tight_layout()
    plt.savefig(f'{output_dir}/id{buffalo_id}_conv_layer1_filters{z}.png', dpi=300, bbox_inches='tight')
    plt.show()


    # -----------------------------------------------------------
    # Loop over each feature map channel and save them as images.
    #    Multiply by x_mask * y_mask if you need to mask out edges.
    # -----------------------------------------------------------

    plt.figure()
    plt.imshow(feat_maps1_sample[z].numpy() * x_mask * y_mask, cmap='viridis')
    plt.title(f"Layer 1, Feature Map {z+1}")
    # Hide axis if you prefer: plt.axis('off')
    plt.savefig(f'{output_dir}/id{buffalo_id}_conv_layer1_feature_map{z}.png', dpi=300, bbox_inches='tight')
    plt.show()

Sequential(
  (0): Conv2d(8, 4, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
  (1): ReLU()
  (2): Conv2d(4, 4, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
  (3): ReLU()
  (4): Conv2d(4, 1, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
)
Filters shape: torch.Size([4, 8, 3, 3])
8

Convolutional layer 2

Activation hook

# -----------------------------------------------------------
# Register a forward hook on the second convolution layer
#    in the model's 'conv_habitat' block
# -----------------------------------------------------------
model.conv_habitat.conv2d[2].register_forward_hook(get_activation("hab_conv2"))

# -----------------------------------------------------------
# Perform a forward pass through the model with the desired input
#    The feature maps from the hooked layer will be stored in 'activation'
# -----------------------------------------------------------
# out = model((x1, x2, x3))  # e.g., model((spatial_data_x, scalars_to_grid, bearing_x))
test = model((sample_spatial_covs, sample_temporal_covs, sample_prev_bearing))

# -----------------------------------------------------------
# Retrieve the captured feature maps from the dictionary
#    and move them to the CPU for inspection
# -----------------------------------------------------------
feat_maps2 = activation["hab_conv2"].cpu()
print("Feature map shape:", feat_maps2.shape)
# Typically shape: (batch_size, out_channels, height, width)

# -----------------------------------------------------------
# Visualize the feature maps for the first sample in the batch
# -----------------------------------------------------------
feat_maps2_sample = feat_maps2[0]  # Shape: (out_channels, H, W)
num_maps2 = feat_maps2_sample.shape[0]
print("Number of feature maps:", num_maps2)

Feature map shape: torch.Size([1, 4, 101, 101])
Number of feature maps: 4

Extract filters and plot

# -------------------------------------------------------------------------
# Extract the weights (filters) from the second convolution layer in conv_habitat
# -------------------------------------------------------------------------
filters_c2 = model.conv_habitat.conv2d[2].weight.data.clone().cpu()
print("Filters shape:", filters_c2.shape)
# Typically (out_channels, in_channels, kernel_height, kernel_width)

# -------------------------------------------------------------------------
# Visualize each filter’s first channel in a grid of subplots
# -------------------------------------------------------------------------
num_filters_c2 = filters_c2.shape[1]
print(num_filters_c2)

for z in range(num_maps2):

    fig, axes = plt.subplots(2, num_filters_c2, figsize=(2*num_filters_c2, 4))
    for i in range(num_filters_c2):

        # Add the covariates as the first row of subplots
        axes[0,i].imshow(feat_maps1_sample[i].numpy() * x_mask * y_mask, cmap='viridis')
        axes[0,i].axis('off')
        axes[0,i].set_title(f"Layer 1, Map {z+1}")

        # if i > 3:
        #     im1 = axes[0,i].imshow(covariate_stack[0, i].detach().cpu().numpy(), cmap='viridis')
        #     im1.set_clim(-1, 1)
        #     axes[0,i].text(scalar_maps.shape[2] // 2, scalar_maps.shape[3] // 2,
        #         f'Value: {round(x2[0, i-4].item(), 2)}',
        #         ha='center', va='center', color='white', fontsize=12)

        kernel = filters_c2[z, i, :, :]  # Show the first input channel
        im = axes[1,i].imshow(kernel, cmap='viridis')
        axes[1,i].axis('off')
        axes[1,i].set_title(f'Layer 2, Filter {z+1}')
        # Annotate each cell with the numeric value
        for (j, k), val in np.ndenumerate(kernel):
            axes[1,i].text(k, j, f'{val:.2f}', ha='center', va='center', color='white')

    plt.tight_layout()
    plt.savefig(f'{output_dir}/id{buffalo_id}_conv_layer2_filters{z}.png', dpi=300, bbox_inches='tight')
    plt.show()


    # -----------------------------------------------------------
    # 6. Loop over each feature map channel and save them as images.
    #    Multiply by x_mask * y_mask if you need to mask out edges.
    # -----------------------------------------------------------

    plt.figure()
    plt.imshow(feat_maps2_sample[z].numpy() * x_mask * y_mask, cmap='viridis')
    plt.title(f"Layer 2, Feature Map {z+1}")
    # Hide axis if you prefer: plt.axis('off')
    plt.savefig(f'{output_dir}/id{buffalo_id}_conv_layer2_feature_map{z}.png', dpi=300, bbox_inches='tight')
    plt.show()

Filters shape: torch.Size([4, 4, 3, 3])
4

Convolutional layer 3

Activation hook

# -----------------------------------------------------------
# Register a forward hook on the third convolution layer
#    in the model's 'conv_habitat' block
# -----------------------------------------------------------
model.conv_habitat.conv2d[4].register_forward_hook(get_activation("hab_conv3"))

# -----------------------------------------------------------
# Perform a forward pass through the model with the desired input
#    The feature maps from the hooked layer will be stored in 'activation'
# -----------------------------------------------------------
# out = model((x1, x2, x3))  # e.g., model((spatial_data_x, scalars_to_grid, bearing_x))
test = model((sample_spatial_covs, sample_temporal_covs, sample_prev_bearing))

# -----------------------------------------------------------
# Retrieve the captured feature maps from the dictionary
#    and move them to the CPU for inspection
# -----------------------------------------------------------
feat_maps3 = activation["hab_conv3"].cpu()
print("Feature map shape:", feat_maps3.shape)
# Typically shape: (batch_size, out_channels, height, width)

# -----------------------------------------------------------
# Visualize the feature maps for the first sample in the batch
# -----------------------------------------------------------
feat_maps3_sample = feat_maps3[0]  # Shape: (out_channels, H, W)
num_maps3 = feat_maps3_sample.shape[0]
print("Number of feature maps:", num_maps3)

Feature map shape: torch.Size([1, 1, 101, 101])
Number of feature maps: 1

Extract filters and plot

# -------------------------------------------------------------------------
# Extract the weights (filters) from the second convolution layer in conv_habitat
# -------------------------------------------------------------------------
filters_c3 = model.conv_habitat.conv2d[4].weight.data.clone().cpu()
print("Filters shape:", filters_c3.shape)
# Typically (out_channels, in_channels, kernel_height, kernel_width)

# -------------------------------------------------------------------------
# Visualize each filter’s first channel in a grid of subplots
# -------------------------------------------------------------------------
num_filters_c3 = filters_c3.shape[1]
print(num_filters_c3)

for z in range(num_maps3):

    fig, axes = plt.subplots(2, num_filters_c3, figsize=(2*num_filters_c3, 4))
    for i in range(num_filters_c3):

        # Add the covariates as the first row of subplots
        axes[0,i].imshow(feat_maps2_sample[i].numpy() * x_mask * y_mask, cmap='viridis')
        axes[0,i].axis('off')
        axes[0,i].set_title(f"Layer 2, Map {z+1}")


        kernel = filters_c3[z, i, :, :]  # Show the first input channel
        im = axes[1,i].imshow(kernel, cmap='viridis')
        axes[1,i].axis('off')
        axes[1,i].set_title(f'Layer 3, Filter {z+1}')
        # Annotate each cell with the numeric value
        for (j, k), val in np.ndenumerate(kernel):
            axes[1,i].text(k, j, f'{val:.2f}', ha='center', va='center', color='white')

    plt.tight_layout()
    plt.savefig(f'{output_dir}/id{buffalo_id}_conv_layer3_filters{z}.png', dpi=300, bbox_inches='tight')
    plt.show()


    # -----------------------------------------------------------
    # 6. Loop over each feature map channel and save them as images.
    #    Multiply by x_mask * y_mask if you need to mask out edges.
    # -----------------------------------------------------------

    plt.figure()
    plt.imshow(feat_maps3_sample[z].numpy() * x_mask * y_mask, cmap='viridis')
    plt.title(f"Habitat selection log probability")
    # Hide axis if you prefer: plt.axis('off')
    plt.savefig(f'{output_dir}/id{buffalo_id}_conv_layer3_feature_map{z}.png', dpi=300, bbox_inches='tight')
    plt.show()

Filters shape: torch.Size([1, 4, 3, 3])
4

Checking estimated movement parameters

Similarly to the convolutional layers, we can set hooks to extract the predicted movement parameters from the model, and assess how variable that is across samples.

torch.cuda.empty_cache()

# -------------------------------------------------------------------------
# Create a list to store the intermediate output from the fully connected
#    movement sub-network (fcn_movement_all)
# -------------------------------------------------------------------------
intermediate_output = []

def hook(module, input, output):
    """
    Hook function that captures the output of the specified layer
    (fcn_movement_all) during the forward pass.
    """
    intermediate_output.append(output)

# -------------------------------------------------------------------------
# Register the forward hook on 'fcn_movement_all', so its outputs
#    are recorded every time the model does a forward pass.
# -------------------------------------------------------------------------
hook_handle = model.fcn_movement_all.register_forward_hook(hook)

# -------------------------------------------------------------------------
# Perform a forward pass with the model in evaluation mode,
#    disabling gradient computation.
# -------------------------------------------------------------------------
model.eval()
with torch.no_grad():
    final_output = model((sample_spatial_covs, sample_temporal_covs, sample_prev_bearing))

# -------------------------------------------------------------------------
# Inspect the captured intermediate output
#    'intermediate_output[0]' corresponds to the first (and only) forward pass.
# -------------------------------------------------------------------------
print("Intermediate output shape:", intermediate_output[0].shape)
print("Intermediate output values:", intermediate_output[0][0])

# -------------------------------------------------------------------------
# Remove the hook to avoid repeated capturing in subsequent passes
# -------------------------------------------------------------------------
hook_handle.remove()

# -------------------------------------------------------------------------
# Unpack the parameters from the FCN output (assumes a specific ordering)
# -------------------------------------------------------------------------
gamma_shape1, gamma_scale1, gamma_weight1, \
gamma_shape2, gamma_scale2, gamma_weight2, \
vonmises_mu1, vonmises_kappa1, vonmises_weight1, \
vonmises_mu2, vonmises_kappa2, vonmises_weight2 = intermediate_output[0][0]

# -------------------------------------------------------------------------
# Convert parameters from log-space (if applicable) and print them
#    Gamma and von Mises parameters
# -------------------------------------------------------------------------
# --- Gamma #1 ---
print("Gamma shape 1:", torch.exp(gamma_shape1))
print("Gamma scale 1:", torch.exp(gamma_scale1))
print("Gamma weight 1:",
      torch.exp(gamma_weight1) / (torch.exp(gamma_weight1) + torch.exp(gamma_weight2)))

# --- Gamma #2 ---
print("Gamma shape 2:", torch.exp(gamma_shape2))
print("Gamma scale 2:", torch.exp(gamma_scale2) * 500)  # scale factor 500
print("Gamma weight 2:",
      torch.exp(gamma_weight2) / (torch.exp(gamma_weight1) + torch.exp(gamma_weight2)))

# --- von Mises #1 ---
# % (2*np.pi) ensures the mu (angle) is wrapped within [0, 2π)
print("Von Mises mu 1:", vonmises_mu1 % (2*np.pi))
print("Von Mises kappa 1:", torch.exp(vonmises_kappa1))
print("Von Mises weight 1:",
      torch.exp(vonmises_weight1) / (torch.exp(vonmises_weight1) + torch.exp(vonmises_weight2)))

# --- von Mises #2 ---
print("Von Mises mu 2:", vonmises_mu2 % (2*np.pi))
print("Von Mises kappa 2:", torch.exp(vonmises_kappa2))
print("Von Mises weight 2:",
      torch.exp(vonmises_weight2) / (torch.exp(vonmises_weight1) + torch.exp(vonmises_weight2)))

Intermediate output shape: torch.Size([1, 12])
Intermediate output values: tensor([ 0.2050,  0.4531, -0.0247, -0.4483, -0.1241,  0.0866, -0.1144, -0.5838,
        -0.1406,  0.1756, -0.6430,  0.2865])
Gamma shape 1: tensor(1.2275)
Gamma scale 1: tensor(1.5732)
Gamma weight 1: tensor(0.4722)
Gamma shape 2: tensor(0.6387)
Gamma scale 2: tensor(441.6361)
Gamma weight 2: tensor(0.5278)
Von Mises mu 1: tensor(6.1687)
Von Mises kappa 1: tensor(0.5578)
Von Mises weight 1: tensor(0.3948)
Von Mises mu 2: tensor(0.1756)
Von Mises kappa 2: tensor(0.5257)
Von Mises weight 2: tensor(0.6052)

Plot the movement distributions

We can use the movement parameters to plot the step length and turning angle distributions for the sample covariates.

# -------------------------------------------------------------------------
# Define helper functions for calculating Gamma and von Mises log-densities
# -------------------------------------------------------------------------
def gamma_density(x, shape, scale):
    """
    Computes the log of the Gamma density for each value in x.

    Args:
      x (Tensor): Input values for which to compute the density.
      shape (float): Gamma shape parameter
      scale (float): Gamma scale parameter

    Returns:
      Tensor: The log of the Gamma probability density at each x.
    """
    return -1*torch.lgamma(shape) - shape*torch.log(scale) \
           + (shape - 1)*torch.log(x) - x/scale

def vonmises_density(x, kappa, vm_mu):
    """
    Computes the log of the von Mises density for each value in x.

    Args:
      x (Tensor): Input angles in radians.
      kappa (float): Concentration parameter (kappa)
      vm_mu (float): Mean direction parameter (mu)

    Returns:
      Tensor: The log of the von Mises probability density at each x.
    """
    return kappa*torch.cos(x - vm_mu) - 1*(np.log(2*torch.pi) + torch.log(torch.special.i0(kappa)))


# -------------------------------------------------------------------------
# Round and display the mixture weights for the Gamma distributions
# -------------------------------------------------------------------------
gamma_weight1_recovered = torch.exp(gamma_weight1)/(torch.exp(gamma_weight1) + torch.exp(gamma_weight2))
rounded_gamma_weight1 = round(gamma_weight1_recovered.item(), 2)

gamma_weight2_recovered = torch.exp(gamma_weight2)/(torch.exp(gamma_weight1) + torch.exp(gamma_weight2))
rounded_gamma_weight2 = round(gamma_weight2_recovered.item(), 2)

# -------------------------------------------------------------------------
# Round and display the mixture weights for the von Mises distributions
# -------------------------------------------------------------------------
vonmises_weight1_recovered = torch.exp(vonmises_weight1)/(torch.exp(vonmises_weight1) + torch.exp(vonmises_weight2))
rounded_vm_weight1 = round(vonmises_weight1_recovered.item(), 2)

vonmises_weight2_recovered = torch.exp(vonmises_weight2)/(torch.exp(vonmises_weight1) + torch.exp(vonmises_weight2))
rounded_vm_weight2 = round(vonmises_weight2_recovered.item(), 2)


# -------------------------------------------------------------------------
# 1. Plotting the Gamma mixture distribution
#    a) Generate x values
#    b) Compute individual Gamma log densities
#    c) Exponentiate and combine using recovered weights
# -------------------------------------------------------------------------
x_values = torch.linspace(1, 101, 1000).to(device)
gamma1_density = gamma_density(x_values, torch.exp(gamma_shape1), torch.exp(gamma_scale1))
gamma2_density = gamma_density(x_values, torch.exp(gamma_shape2), torch.exp(gamma_scale2)*500)
gamma_mixture_density = gamma_weight1_recovered*torch.exp(gamma1_density) \
                        + gamma_weight2_recovered*torch.exp(gamma2_density)

# Move results to CPU and convert to NumPy for plotting
x_values_np = x_values.cpu().numpy()
gamma1_density_np = np.exp(gamma1_density.cpu().numpy())
gamma2_density_np = np.exp(gamma2_density.cpu().numpy())
gamma_mixture_density_np = gamma_mixture_density.cpu().numpy()

# -------------------------------------------------------------------------
# 2. Plot the Gamma distributions and their mixture
# -------------------------------------------------------------------------
filename_gamma_distributions = f'{output_dir}/id{buffalo_id}_gamma_distributions.png'

plt.plot(x_values_np, gamma1_density_np, label=f'Gamma 1 Density: weight = {rounded_gamma_weight1}')
plt.plot(x_values_np, gamma2_density_np, label=f'Gamma 2 Density: weight = {rounded_gamma_weight2}')
plt.plot(x_values_np, gamma_mixture_density_np, label='Gamma Mixture Density')
plt.xlabel('x')
plt.ylabel('Density')
plt.title('Gamma Density Function')
plt.legend()
plt.savefig(filename_gamma_distributions, dpi=300, bbox_inches='tight')
plt.show()


# -------------------------------------------------------------------------
# 3. Plotting the von Mises mixture distribution
#    a) Generate x values from -π to π
#    b) Compute individual von Mises log densities
#    c) Exponentiate and combine using recovered weights
# -------------------------------------------------------------------------
x_values = torch.linspace(-np.pi, np.pi, 1000).to(device)
vonmises1_density = vonmises_density(x_values, torch.exp(vonmises_kappa1), vonmises_mu1)
vonmises2_density = vonmises_density(x_values, torch.exp(vonmises_kappa2), vonmises_mu2)
vonmises_mixture_density = vonmises_weight1_recovered*torch.exp(vonmises1_density) \
                           + vonmises_weight2_recovered*torch.exp(vonmises2_density)

# Move results to CPU and convert to NumPy for plotting
x_values_np = x_values.cpu().numpy()
vonmises1_density_np = np.exp(vonmises1_density.cpu().numpy())
vonmises2_density_np = np.exp(vonmises2_density.cpu().numpy())
vonmises_mixture_density_np = vonmises_mixture_density.cpu().numpy()

# -------------------------------------------------------------------------
# 4. Plot the von Mises distributions and their mixture
# -------------------------------------------------------------------------
filename_vonmises_distributions = f'{output_dir}/id{buffalo_id}_vonmises_distributions.png'

plt.plot(x_values_np, vonmises1_density_np, label=f'Von Mises 1 Density: weight = {rounded_vm_weight1}')
plt.plot(x_values_np, vonmises2_density_np, label=f'Von Mises 2 Density: weight = {rounded_vm_weight2}')
plt.plot(x_values_np, vonmises_mixture_density_np, label='Von Mises Mixture Density')
plt.xlabel('x (radians)')
plt.ylabel('Density')
plt.title('Von Mises Density Function')
plt.ylim(0, 0.5)  # Set a limit for the y-axis
plt.legend()
plt.savefig(filename_vonmises_distributions, dpi=300, bbox_inches='tight')
plt.show()

Generate a distribution of movement parameters

To see how variable the movement parameters are across samples, we can generate a distribution of movement parameters from a batch of samples.

We take the code from above that we used to create the DataLoader for the test data and increase the batch size (to get more samples to create the distribution from).

As we’re not using the test dataset any more, we’ll just put all of the samples in the same batch, and generate movement parameters for all of them.

# To use all of the test samples
print(f'There are {len(dataset_test)} samples in the test dataset')
bs = len(dataset_test) # batch size
dataloader_test = DataLoader(dataset=dataset_test, batch_size=bs, shuffle=True)

torch.cuda.empty_cache()

There are 1011 samples in the test dataset

Take all of the samples from the test dataset and put them in a single batch.

# -----------------------------------------------------------
# Fetch a batch of data from the training dataloader
# -----------------------------------------------------------
x1_batch, x2_batch, x3_batch, labels = next(iter(dataloader_test))

# Move the input batches to the same device as the model
x1_batch = x1_batch.to(device)
x2_batch = x2_batch.to(device)
x3_batch = x3_batch.to(device)

# -----------------------------------------------------------
# Register a forward hook to capture the outputs
#    from 'fcn_movement_all' during the forward pass
# -----------------------------------------------------------
hook_handle = model.fcn_movement_all.register_forward_hook(hook)

# -----------------------------------------------------------
# Perform a forward pass in evaluation mode to generate
#    and capture the sub-network's outputs in 'intermediate_output'
# -----------------------------------------------------------
model.eval()  # Disables certain layers like dropout

# Pass the batch through the model
final_output = model((x1_batch, x2_batch, x3_batch))

# -----------------------------------------------------------
# Prepare lists to store the distribution parameters
#    for each sample in the batch
# -----------------------------------------------------------
gamma_shape1_list = []
gamma_scale1_list = []
gamma_weight1_list = []
gamma_shape2_list = []
gamma_scale2_list = []
gamma_weight2_list = []
vonmises_mu1_list = []
vonmises_kappa1_list = []
vonmises_weight1_list = []
vonmises_mu2_list = []
vonmises_kappa2_list = []
vonmises_weight2_list = []

# -----------------------------------------------------------
# Extract parameters from 'intermediate_output'
#    for every sample in the batch
# -----------------------------------------------------------
for batch_output in intermediate_output:
    # Each 'batch_output' corresponds to one forward pass;
    # it might contain multiple samples if the batch size > 1
    for sample_output in batch_output:
        # Unpack the 12 parameters of the Gamma and von Mises mixtures
        gamma_shape1, gamma_scale1, gamma_weight1, \
        gamma_shape2, gamma_scale2, gamma_weight2, \
        vonmises_mu1, vonmises_kappa1, vonmises_weight1, \
        vonmises_mu2, vonmises_kappa2, vonmises_weight2 = sample_output

        # Convert log-space parameters to real space, then store
        gamma_shape1_list.append(torch.exp(gamma_shape1).item())
        gamma_scale1_list.append(torch.exp(gamma_scale1).item())
        gamma_weight1_list.append(
            (torch.exp(gamma_weight1)/(torch.exp(gamma_weight1) + torch.exp(gamma_weight2))).item()
        )
        gamma_shape2_list.append(torch.exp(gamma_shape2).item())
        gamma_scale2_list.append((torch.exp(gamma_scale2)*500).item())  # scale factor 500
        gamma_weight2_list.append(
            (torch.exp(gamma_weight2)/(torch.exp(gamma_weight1) + torch.exp(gamma_weight2))).item()
        )
        vonmises_mu1_list.append((vonmises_mu1 % (2*np.pi)).item())
        vonmises_kappa1_list.append(torch.exp(vonmises_kappa1).item())
        vonmises_weight1_list.append(
            (torch.exp(vonmises_weight1)/(torch.exp(vonmises_weight1) + torch.exp(vonmises_weight2))).item()
        )
        vonmises_mu2_list.append((vonmises_mu2 % (2*np.pi)).item())
        vonmises_kappa2_list.append(torch.exp(vonmises_kappa2).item())
        vonmises_weight2_list.append(
            (torch.exp(vonmises_weight2)/(torch.exp(vonmises_weight1) + torch.exp(vonmises_weight2))).item()
        )

Plot the distribution of movement parameters

# -----------------------------------------------------------
# Define a helper function to plot histograms
#    for the collected parameters
# -----------------------------------------------------------
def plot_histogram(data, title, xlabel):
    """
    Plots a histogram of the provided data.

    Args:
        data (list): Data points to plot in a histogram.
        title (str): Title of the histogram plot.
        xlabel (str): X-axis label.
    """
    plt.figure()
    plt.hist(data, bins=30, alpha=0.75)
    plt.title(title)
    plt.xlabel(xlabel)
    plt.ylabel('Frequency')

    # Generate a filename from the title (replace spaces with underscores)
    filename = title.replace(' ', '_') + '.png'
    filepath = os.path.join(output_dir, filename)

    # Save the figure
    plt.savefig(filepath, dpi=300, bbox_inches='tight')
    plt.show()
    plt.close()  # Close the figure to free memory\

# -----------------------------------------------------------
# Plot histograms for each parameter distribution
# -----------------------------------------------------------
plot_histogram(gamma_shape1_list, 'Gamma Shape 1 Distribution', 'Shape 1')
plot_histogram(gamma_scale1_list, 'Gamma Scale 1 Distribution', 'Scale 1')
plot_histogram(gamma_weight1_list, 'Gamma Weight 1 Distribution', 'Weight 1')
plot_histogram(gamma_shape2_list, 'Gamma Shape 2 Distribution', 'Shape 2')
plot_histogram(gamma_scale2_list, 'Gamma Scale 2 Distribution', 'Scale 2')
plot_histogram(gamma_weight2_list, 'Gamma Weight 2 Distribution', 'Weight 2')
plot_histogram(vonmises_mu1_list, 'Von Mises Mu 1 Distribution', 'Mu 1')
plot_histogram(vonmises_kappa1_list, 'Von Mises Kappa 1 Distribution', 'Kappa 1')
plot_histogram(vonmises_weight1_list, 'Von Mises Weight 1 Distribution', 'Weight 1')
plot_histogram(vonmises_mu2_list, 'Von Mises Mu 2 Distribution', 'Mu 2')
plot_histogram(vonmises_kappa2_list, 'Von Mises Kappa 2 Distribution', 'Kappa 2')
plot_histogram(vonmises_weight2_list, 'Von Mises Weight 2 Distribution', 'Weight 2')

# -----------------------------------------------------------
# Remove the hook to stop capturing outputs
#    in subsequent forward passes
# -----------------------------------------------------------
hook_handle.remove()

Validate the predictions from the model

We want to assess how well the model is generating movement, habitat selection and next-step probabilities, so we will extract those probabilities at the location of each observed step.

Importing spatial data

Instead of importing the stacks of local layers (one for each step), here we want to import the spatial covariates for the extent we want to simulate over. We use an extent that covers all of the observed locations, which refer to as the ‘landscape’.

NDVI

We have monthly NDVI layers for 2018 and 2019, which we import as a stack. The files don’t import with a time component, so we will use a function further down that indexes them correctly.

# for monthly NDVI
file_path = f'{base_path}/mapping/cropped rasters/ndvi_monthly.tif'

# read the raster file
with rasterio.open(file_path) as src:
    # Read the raster band as separate variable
    ndvi_landscape = src.read([i for i in range(1, src.count + 1)])
    # Get the metadata of the raster
    ndvi_meta = src.meta
    raster_transform = src.transform

    # Print the metadata to check for time component
    print("Metadata:", ndvi_meta)

    # Check for specific time-related metadata
    if 'TIFFTAG_DATETIME' in src.tags():
        print("Time component found:", src.tags()['TIFFTAG_DATETIME'])
    else:
        print("No explicit time component found in metadata.")

# the rasters don't contain a time component, so we will use a function later to index the layers correctly

Metadata: {'driver': 'GTiff', 'dtype': 'float32', 'nodata': nan, 'width': 2400, 'height': 2280, 'count': 24, 'crs': CRS.from_wkt('LOCAL_CS["GDA94 / Geoscience Australia Lambert",UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'), 'transform': Affine(25.0, 0.0, 0.0,
       0.0, -25.0, -1406000.0)}
No explicit time component found in metadata.

Prepare the NDVI data

There are a few things we need to do to prepare the landscape layers.

First, we need to ensure that there are no NA values in the data. For NDVI we will replace any NA values with -1 (which denotes water), as in our case that is typically why they were set to NA.

Secondly, the model expects the covariates to on the same scale as the training data. We will therefore scale the NDVI data using the same max and min scaling parameters as the training data. To get these, there are some min and max print statements in the deepSSF_train.ipynb script. When we plot the NDVI data below we will see that the values will no longer range from 0 to 1, which is because there are values in the landscape layers that are outside of the range of the training data.

# Check the coordinate reference system
print("NDVI metadata:")
print(ndvi_meta)
print("\n")

# Have a look at the affine transformation parameters that are used to convert pixel
# coordinates to geographic coordinates and vice versa
print("Affine transformation parameters:")
print(raster_transform)
print("\n")

# Check the shape (layers, row, columns) of the raster
print("Shape of the raster:")
print(ndvi_landscape.shape)

# Replace NaNs in the original array with -1, which represents water
ndvi_landscape = np.nan_to_num(ndvi_landscape, nan=-1.0)

# from the stack of local layers (training data)
ndvi_max = 0.8220
ndvi_min = -0.2772

# Convert the numpy array to a PyTorch tensor
ndvi_landscape_tens = torch.from_numpy(ndvi_landscape)

# Normalizing the data
ndvi_landscape_norm = (ndvi_landscape_tens - ndvi_min) / (ndvi_max - ndvi_min)

# Show two example layers of the scaled NDVI data
layer_index = 1
plt.imshow(ndvi_landscape_norm[layer_index,:,:].numpy())
plt.colorbar()
plt.title(f'NDVI layer index {layer_index}')
plt.show()

layer_index = 8
plt.imshow(ndvi_landscape_norm[layer_index,:,:].numpy())
plt.colorbar()
plt.title(f'NDVI layer index {layer_index}')
plt.show()

NDVI metadata:
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': nan, 'width': 2400, 'height': 2280, 'count': 24, 'crs': CRS.from_wkt('LOCAL_CS["GDA94 / Geoscience Australia Lambert",UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'), 'transform': Affine(25.0, 0.0, 0.0,
       0.0, -25.0, -1406000.0)}


Affine transformation parameters:
| 25.00, 0.00, 0.00|
| 0.00,-25.00,-1406000.00|
| 0.00, 0.00, 1.00|


Shape of the raster:
(24, 2280, 2400)

Canopy cover

Canopy cover is just a single static layer.

# Path to the canopy cover raster file
file_path = f'{base_path}/mapping/cropped rasters/canopy_cover.tif'

# read the raster file
with rasterio.open(file_path) as src:
    # Read the raster band as separate variable
    canopy_landscape = src.read(1)
    # Get the metadata of the raster
    canopy_meta = src.meta

Prepare the canopy cover data

As with the NDVI data, we need to ensure that there are no NA values in the data.

As the canopy cover values in the landscape layer are within the same range as the training data, we see that the values range from 0 to 1.

# Check the canopy metadata:
print("Canopy metadata:")
print(canopy_meta)
print("\n")

# Check the shape (rows, columns) of the canopy raster:
print("Shape of canopy raster:")
print(canopy_landscape.shape)
print("\n")

# Check for NA values in the canopy raster:
print("Number of NA values in the canopy raster:")
print(np.isnan(canopy_landscape).sum())

# Define the maximum and minimum canopy values from the stack of local layers:
canopy_max = 82.5000
canopy_min = 0.0

# Convert the canopy data from a NumPy array to a PyTorch tensor:
canopy_landscape_tens = torch.from_numpy(canopy_landscape)

# Normalise the canopy data:
canopy_landscape_norm = (canopy_landscape_tens - canopy_min) / (canopy_max - canopy_min)

# Visualise the normalised canopy cover:
plt.imshow(canopy_landscape_norm.numpy())
plt.colorbar()
plt.title('Canopy Cover')
plt.show()

Canopy metadata:
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -3.3999999521443642e+38, 'width': 2400, 'height': 2280, 'count': 1, 'crs': CRS.from_wkt('LOCAL_CS["GDA94 / Geoscience Australia Lambert",UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'), 'transform': Affine(25.0, 0.0, 0.0,
       0.0, -25.0, -1406000.0)}


Shape of canopy raster:
(2280, 2400)


Number of NA values in the canopy raster:
0

Herbaceous vegetation

# Path to the herbaceous vegetation raster file
file_path = f'{base_path}/mapping/cropped rasters/veg_herby.tif'

# read the raster file
with rasterio.open(file_path) as src:
    # Read the raster band as separate variable
    herby_landscape = src.read(1)
    # Get the metadata of the raster
    herby_meta = src.meta

# Check the herbaceous metadata:
print("Herbaceous metadata:")
print(herby_meta)
print("\n")

# Check the shape (rows, columns) of the herbaceous raster:
print("Shape of herbaceous raster:")
print(herby_landscape.shape)
print("\n")

# Check for NA values in the herby raster:
print("Number of NA values in the herbaceous vegetation raster:")
print(np.isnan(herby_landscape).sum())

# Define the maximum and minimum herbaceous values from the stack of local layers:
herby_max = 1.0
herby_min = 0.0

# Convert the herbaceous data from a NumPy array to a PyTorch tensor:
herby_landscape_tens = torch.from_numpy(herby_landscape)

# Normalize the herbaceous data:
herby_landscape_norm = (herby_landscape_tens - herby_min) / (herby_max - herby_min)

# Visualize the normalised herbaceous cover:
plt.imshow(herby_landscape_norm.numpy())
plt.colorbar()
plt.show()

Herbaceous metadata:
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -3.3999999521443642e+38, 'width': 2400, 'height': 2280, 'count': 1, 'crs': CRS.from_wkt('LOCAL_CS["GDA94 / Geoscience Australia Lambert",UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'), 'transform': Affine(25.0, 0.0, 0.0,
       0.0, -25.0, -1406000.0)}


Shape of herbaceous raster:
(2280, 2400)


Number of NA values in the herbaceous vegetation raster:
0

Slope

# Path to the slope raster file
file_path = f'{base_path}/mapping/cropped rasters/slope.tif'

# read the raster file
with rasterio.open(file_path) as src:
    # Read the raster band as separate variable
    slope_landscape = src.read(1)
    # Get the metadata of the raster
    slope_meta = src.meta

# Check the slope metadata:
print("Slope metadata:")
print(slope_meta)
print("\n")

# Check the shape (rows, columns) of the slope landscape raster:
print("Shape of slope landscape raster:")
print(slope_landscape.shape)
print("\n")

# Check for NA values in the slope raster:
print("Number of NA values in the slope raster:")
print(np.isnan(slope_landscape).sum())

# Replace NaNs in the slope array with 0.0:
slope_landscape = np.nan_to_num(slope_landscape, nan=0.0)

# Define the maximum and minimum slope values from the stack of local layers:
slope_max = 12.2981
slope_min = 0.0006

# Convert the slope landscape data from a NumPy array to a PyTorch tensor:
slope_landscape_tens = torch.from_numpy(slope_landscape)

# Normalize the slope landscape data:
slope_landscape_norm = (slope_landscape_tens - slope_min) / (slope_max - slope_min)

# Visualize the slope landscape (note: displaying the original tensor, not the normalised data):
plt.imshow(slope_landscape_tens.numpy())
plt.colorbar()
plt.show()

Slope metadata:
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': nan, 'width': 2400, 'height': 2280, 'count': 1, 'crs': CRS.from_wkt('LOCAL_CS["GDA94 / Geoscience Australia Lambert",UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'), 'transform': Affine(25.0, 0.0, 0.0,
       0.0, -25.0, -1406000.0)}


Shape of slope landscape raster:
(2280, 2400)


Number of NA values in the slope raster:
9356

Subset function (with padding)

Now that we have our landscape layers imported, we need a way to crop out the local layers that can be fed into the deepSSF model as covariates.

We will use the same subset function as in the simulation script, which we stored in the deepSSF_utils.py script. This function will take the landscape layers and a set of coordinates, and return the local layers for those coordinates.

This function also has padding for if the simulated individual was to go off the edge of the landscape, which we retain here (although we won’t need that functionality).

subset_raster = deepSSF_utils.subset_raster_with_padding_torch

Testing the subset function

Use the subset function to crop out the local layers for all covariates. Try different locations using the x and y coordinates, which are in geographic coordinates (x = easting/longitude, y = northing/latitude).

# Pick a location (x, y) from the buffalo DataFrame
x = buffalo_df['x1_'].iloc[0]
y = buffalo_df['y1_'].iloc[0]

# Define the size of the window to extract
window_size = 101

# Select the NDVI layer index
which_ndvi = 1

# Extract subsets from various raster layers using the custom function.
# Each call centres the window at the specified (x, y) location and applies padding where necessary.
ndvi_subset, origin_x, origin_y = subset_raster(ndvi_landscape_norm[which_ndvi, :, :],
                                                x, y, window_size, raster_transform)
canopy_subset, origin_x, origin_y = subset_raster(canopy_landscape_norm,
                                                  x, y, window_size, raster_transform)
herby_subset, origin_x, origin_y = subset_raster(herby_landscape_norm,
                                                 x, y, window_size, raster_transform)
slope_subset, origin_x, origin_y = subset_raster(slope_landscape_norm,
                                                 x, y, window_size, raster_transform)

# Create a 2x2 grid of subplots with a fixed figure size.
fig, axs = plt.subplots(2, 2, figsize=(10, 8))

# Plot the NDVI subset.
im0 = axs[0, 0].imshow(ndvi_subset.numpy(), cmap='viridis')
fig.colorbar(im0, ax=axs[0, 0], shrink=0.8)
axs[0, 0].set_title('NDVI Subset')

# Plot the Canopy Cover subset.
im1 = axs[0, 1].imshow(canopy_subset.numpy(), cmap='viridis')
fig.colorbar(im1, ax=axs[0, 1], shrink=0.8)
axs[0, 1].set_title('Canopy Cover Subset')

# Plot the Herbaceous Vegetation subset.
im2 = axs[1, 0].imshow(herby_subset.numpy(), cmap='viridis')
fig.colorbar(im2, ax=axs[1, 0], shrink=0.8)
axs[1, 0].set_title('Herbaceous Vegetation Subset')

# Plot the Slope subset.
im3 = axs[1, 1].imshow(slope_subset.numpy(), cmap='viridis')
fig.colorbar(im3, ax=axs[1, 1], shrink=0.8)
axs[1, 1].set_title('Slope Subset')

Text(0.5, 1.0, 'Slope Subset')

Create a mask for edge cells

Due to the padding at the edges of the covariates, convolutional layers create artifacts that can affect the colour scale of the predictions when plotting. To avoid this, we will create a mask that we can apply to the predictions to remove the edge cells.

# Create a mask to remove the edge values for plotting
# (as it affects the colour scale)
x_mask = np.ones_like(ndvi_subset)
y_mask = np.ones_like(ndvi_subset)

# Mask out bordering cells
x_mask[:, :3] = -np.inf
x_mask[:, 98:] = -np.inf
y_mask[:3, :] = -np.inf
y_mask[98:, :] = -np.inf

Setup simulation parameters

To get the simulation running we need a few extra functions.

Firstly, we need to index the NDVI layers correctly, based on the time of the simulated location. We’ll do this by creating a function that takes day of the year of the simulated location and returns the correct index for the NDVI layers.

Recall that Python indexes from 0, so when the month_index is equal to 2 for instance, this will index the third layer, which is for March.

# Create a mapping from day of the year to month index
def day_to_month_index(day_of_year):
    # Calculate the year and the day within that year
    base_date = datetime(2018, 1, 1) # base date for the calculation, which is when the NDVI layers start
    date = base_date + timedelta(days=int(day_of_year) - 1)
    year_diff = date.year - base_date.year
    month_index = (date.month - 1) + (year_diff * 12)  # month index (0-based, accounting for year change)
    return month_index

yday = 70 # day of the year, which is March 11th
month_index = day_to_month_index(yday)
print(month_index)

2

Next-step probability values

We can now calculate the next-step probabilities for each observed step. As we generate habitat selection, movement and next-step probability surfaces, we can get the predicted probability values for each one, which can be compared to the respective process in the SSF.

The process for generating the next-step probabilities is as follows:

1. Get the current location of the individual
2. Crop out the local layers for the current location
3. Run the model of the local layers to get the habitat selection, movement and next-step probability surfaces
4. Get the predicted probability values at the location of the next step
5. Store the predicted probability values and export them as a csv for comparison with the SSF

First, select the data to generate prediction values for. For testing the function we can select a subset.

# To select a subset of samples to test the function
# test_data = buffalo_df.iloc[0:10]

# To select all of the data
test_data = buffalo_df

# Get the number of samples in the test data
n_samples = len(test_data)
print(f'Number of samples: {n_samples}')

# Create empty vectors to store the predicted probabilities
habitat_probs = np.repeat(0., n_samples)
move_probs = np.repeat(0., n_samples)
next_step_probs = np.repeat(0., n_samples)

Number of samples: 10103

Loop over each step

# Create directory for saving prediction images
os.makedirs(f'{output_dir}/prediction_images', exist_ok=True)

# Start at 1 so the bearing at t - 1 is available
for i in range(1, n_samples):
# for i in range(1, 4):

  sample = test_data.iloc[i]

  # Current location (x1, y1)
  x = sample['x1_']
  y = sample['y1_']

  # Convert geographic coordinates to pixel coordinates
  px, py = ~raster_transform * (x, y)
  # print('px and py are ', px, py) # Debugging

  # Next step location (x2, y2)
  x2 = sample['x2_']
  y2 = sample['y2_']

  # Convert geographic coordinates to pixel coordinates
  px2, py2 = ~raster_transform * (x2, y2)
  # print('px2 and py2 are ', px2, py2) # Debugging

  # The difference in x and y coordinates
  d_x = x2 - x
  d_y = y2 - y
  # print('d_x and d_y are ', d_x, d_y) # Debugging

  # The difference in pixel coordinates
  d_px = px2 - px
  d_py = py2 - py
  # print('d_px and d_py are ', d_px, d_py) # Debugging

  # Temporal covariates
  hour_t2_sin = sample['hour_t2_sin']
  hour_t2_cos = sample['hour_t2_cos']
  yday_t2_sin = sample['yday_t2_sin']
  yday_t2_cos = sample['yday_t2_cos']

  # Bearing of previous step (t - 1)
  bearing = sample['bearing_tm1']

  # Hour of the day (for saving the plot)
  hour_t2 = sample['hour_t2']

  # Day of the year
  yday = sample['yday_t2']

  # Convert day of the year to month index
  month_index = day_to_month_index(yday)
  # print(month_index)

  # Extract the subset of the covariates at the location of x1, y1
  # NDVI
  ndvi_subset, origin_x, origin_y = subset_raster(ndvi_landscape_norm[month_index,:,:],
                                                  x, y,
                                                  window_size,
                                                  raster_transform)

  # Canopy cover
  canopy_subset, origin_x, origin_y = subset_raster(canopy_landscape_norm,
                                                    x, y,
                                                    window_size,
                                                    raster_transform)

  # Herbaceous vegetation
  herby_subset, origin_x, origin_y = subset_raster(herby_landscape_norm,
                                                   x, y,
                                                   window_size,
                                                   raster_transform)

  # Slope
  slope_subset, origin_x, origin_y = subset_raster(slope_landscape_norm,
                                                   x, y,
                                                   window_size,
                                                   raster_transform)

  # Location of the current step in local pixel coordinates
  px_subset = px - origin_x
  py_subset = py - origin_y
  # print('px_subset and py_subset are ', px_subset, py_subset) # Debugging

  # Location of the next step in local pixel coordinates
  px2_subset = px2 - origin_x
  py2_subset = py2 - origin_y
  # print('px2_subset and py2_subset are ', px2_subset, py2_subset, '\n') # Debugging

  # print(int(py2_subset), int(px2_subset))

  # Location of the next step in local pixel coordinates
  px2_subset_corrected = (px2 - px) + (px - origin_x)
  py2_subset_corrected = (py2 - py) + (py - origin_y)
  # print('px2_subset_corrected and py2_subset_corrected are ', px2_subset_corrected, py2_subset_corrected, '\n') # Debugging

  # Extract the value of the covariates at the location of x2, y2
  # value = ndvi_subset.detach().cpu().numpy()[(int(py2_subset), int(px2_subset))]

  # Stack the channels along a new axis
  x1 = torch.stack([ndvi_subset, canopy_subset, herby_subset, slope_subset], dim=0)

  # Add a batch dimension (required to be the correct dimension for the model)
  x1 = x1.unsqueeze(0).to(device)
  # print(x1.shape)

  # Convert lists to PyTorch tensors
  hour_t2_sin_tensor = torch.tensor(hour_t2_sin).float()
  hour_t2_cos_tensor = torch.tensor(hour_t2_cos).float()
  yday_t2_sin_tensor = torch.tensor(yday_t2_sin).float()
  yday_t2_cos_tensor = torch.tensor(yday_t2_cos).float()

  # Stack tensors
  x2 = torch.stack((hour_t2_sin_tensor.unsqueeze(0),
                    hour_t2_cos_tensor.unsqueeze(0),
                    yday_t2_sin_tensor.unsqueeze(0),
                    yday_t2_cos_tensor.unsqueeze(0)),
                    dim=1)
  # print(x2)
  # print(x2.shape)
  x2 = x2.to(device)

  # Put bearing in the correct dimension (batch_size, 1)
  bearing = torch.tensor(bearing).float().unsqueeze(0).unsqueeze(0)
  # print(bearing)
  # print(bearing.shape)
  bearing = bearing.to(device)


  # -------------------------------------------------------------------------
  # Run the model
  # -------------------------------------------------------------------------
  model_output = model((x1, x2, bearing))


  # -------------------------------------------------------------------------
  # Habitat selection probability
  # -------------------------------------------------------------------------
  hab_density = model_output.detach().cpu().numpy()[0,:,:,0]
  hab_density_exp = np.exp(hab_density)

  # Normalise the probability surface to sum to 1
  hab_density_exp_norm = hab_density_exp / np.sum(hab_density_exp)
  # print(np.sum(hab_density_exp_norm))  # Should be 1

  # Store the probability of habitat selection at the location of x2, y2
  # These probabilities are normalised in the model function
  habitat_probs[i] = hab_density_exp_norm[(int(py2_subset), int(px2_subset))]
  # print('Habitat probability value = ', habitat_probs[i])


  # -------------------------------------------------------------------------
  # Movement probability
  # -------------------------------------------------------------------------
  move_density = model_output.detach().cpu().numpy()[0,:,:,1]
  move_density_exp = np.exp(move_density)

  # Normalise the probability surface to sum to 1
  move_density_exp_norm = move_density_exp / np.sum(move_density_exp)
  # print(np.sum(move_density_exp_norm))  # Should be 1

  # Store the movement probability at the location of x2, y2
  # These probabilities are normalised in the model function
  move_probs[i] = move_density_exp_norm[(int(py2_subset), int(px2_subset))]
  # print('Movement probability value = ', move_probs[i])


  # -------------------------------------------------------------------------
  # Next step probability
  # -------------------------------------------------------------------------
  step_density = hab_density + move_density
  step_density_exp = np.exp(step_density)
  # print('Sum of step density exp = ', np.sum(step_density_exp)) # Won't be 1

  step_density_exp_norm = step_density_exp / np.sum(step_density_exp)
  # print('Sum of step density exp norm = ', np.sum(step_density_exp_norm)) # Should be 1

  # Extract the value of the covariates at the location of x2, y2
  next_step_probs[i] = step_density_exp_norm[(int(py2_subset), int(px2_subset))]
  # print('Next-step probability value = ', next_step_probs[i])


  # -------------------------------------------------------------------------
  # Plot the next-step predictions
  # -------------------------------------------------------------------------

  # Plot the first few probability surfaces - change the condition to i < n_steps to plot all
  if i < 51:

    # Mask out bordering cells
    hab_density_mask = np.log(hab_density_exp_norm) * x_mask * y_mask
    move_density_mask = np.log(move_density_exp_norm) * x_mask * y_mask
    step_density_mask = np.log(step_density_exp_norm) * x_mask * y_mask

    # Create a mask for the next step
    next_step_mask = np.ones_like(hab_density)
    next_step_mask[int(py2_subset), int(px2_subset)] = -np.inf

    # Plot the outputs
    fig_out, axs_out = plt.subplots(2, 2, figsize=(10, 8))

    # Plot NDVI
    im1 = axs_out[0, 0].imshow(ndvi_subset.numpy(), cmap='viridis')
    axs_out[0, 0].set_title('NDVI')
    fig_out.colorbar(im1, ax=axs_out[0, 0], shrink=0.7)

    # Plot habitat selection log-probability
    im2 = axs_out[0, 1].imshow(hab_density_mask * next_step_mask, cmap='viridis')
    axs_out[0, 1].set_title('Habitat selection log-probability')
    fig_out.colorbar(im2, ax=axs_out[0, 1], shrink=0.7)

    # Movement density log-probability
    im3 = axs_out[1, 0].imshow(move_density_mask * next_step_mask, cmap='viridis')
    axs_out[1, 0].set_title('Movement log-probability')
    fig_out.colorbar(im3, ax=axs_out[1, 0], shrink=0.7)

    # Next-step probability
    im4 = axs_out[1, 1].imshow(step_density_mask * next_step_mask, cmap='viridis')
    axs_out[1, 1].set_title('Next-step log-probability')
    fig_out.colorbar(im4, ax=axs_out[1, 1], shrink=0.7)

    filename_covs = f'{output_dir}/prediction_images/id{buffalo_id}_step{i+1}_yday{yday}_hour{hour_t2}.png'
    plt.tight_layout()
    plt.savefig(filename_covs, dpi=300)
    # plt.show()
    plt.close()  # Close the figure to free memory

print(next_step_probs)

[0.00000000e+00 2.81957048e-03 5.07908035e-03 ... 3.43296617e-01
 3.20735824e-04 1.59874879e-04]

Create gif of next-step predictions against the observed data

Change the extract epoch to extract step for ordering the images.

# Example sorting by the epoch number
def extract_step(filename):
    # Extract the epoch number from the filename
    # Adjust the extraction based on your naming pattern
    import re
    match = re.search(r'_step(\d+)_', filename)
    if match:
        return int(match.group(1))
    return 0

def create_gif(image_folder, output_filename, fps=10):
    """
    Creates a GIF from a sequence of images in a folder.

    Parameters:
    - image_folder: Path to the folder containing images
    - output_filename: Name of the output GIF file
    - duration: Duration of each frame in seconds
    """
    # Get all png files in the specified folder, sorted by name
    images = sorted(glob.glob(os.path.join(image_folder, '*.png')), key=extract_step)

    # Check if any images were found
    if not images:
        print(f"No images found in {image_folder}")
        return

    # Read all images
    frames = [imageio.imread(image) for image in images]

    # Save as GIF
    imageio.mimsave(output_filename, frames, fps=fps, loop=0)

    display(Image(filename=output_filename))

    print(f"GIF created successfully: {output_filename}")

# Path to your images
image_folder =  f'{output_dir}/prediction_images'

# Output GIF filename
output_filename = f'{output_dir}/prediction_gif_id{buffalo_id}_yday{yday_t2_integer}_hour{hour_t2_integer}_bearing{bearing_degrees}_next_r{row}_c{column}.gif'

# Create the GIF
create_gif(image_folder, output_filename, fps=5)

<IPython.core.display.Image object>

GIF created successfully: ../Python/outputs/model_training/id2005_deepSSF_training_1_2025-07-09/prediction_gif_id2005_yday324_hour17_bearing156_next_r46_c44.gif

Calculate the null probabilities

As each cell has a probability values, we can calculate what the probability would be if the model provided no information at all, and each cell was equally likely to be the next step. This is just 1 divided by the total number of cells.

null_prob = 1 / (window_size ** 2)
print(f'Null probability: {null_prob:.3e}')

Null probability: 9.803e-05

Compute the rolling average of the probabilities

rolling_window_size = 100 # Rolling window size

# Convert to pandas Series and compute rolling mean
rolling_mean_habitat = pd.Series(habitat_probs).rolling(window=window_size, center=True).mean()
rolling_mean_movement = pd.Series(move_probs).rolling(window=window_size, center=True).mean()
rolling_mean_next_step = pd.Series(next_step_probs).rolling(window=window_size, center=True).mean()

Plot the probabilities

We can get an idea of how variable the probabilities are for the habitat selection and movement surfaces, and for the next-step probabilities, by plotting them across the trajectory

# Plot the habitat probs through time as a line graph
plt.plot(habitat_probs[habitat_probs > 0], color='blue', label='Habitat Probabilities')
plt.plot(rolling_mean_habitat[rolling_mean_habitat > 0], color='red', label='Rolling Mean')
plt.axhline(y=null_prob, color='black', linestyle='--', label='Null Probability')  # null probs
plt.xlabel('Index')
plt.ylabel('Probability')
plt.title('Habitat Probability')
plt.legend()  # Add legend to differentiate lines
plt.savefig(f'{output_dir}/id{buffalo_id}_habitat_probs.png', dpi=300, bbox_inches='tight')
plt.show()

# Plot the movement probs through time as a line graph
plt.plot(move_probs[move_probs > 0], color='blue', label='Movement Probabilities')
plt.plot(rolling_mean_movement[rolling_mean_movement > 0], color='red', label='Rolling Mean')
plt.axhline(y=null_prob, color='black', linestyle='--', label='Null Probability')  # null probs
plt.xlabel('Index')
plt.ylabel('Probability')
plt.title('Movement Probability')
plt.legend()  # Add legend to differentiate lines
plt.savefig(f'{output_dir}/id{buffalo_id}_movement_probs.png', dpi=300, bbox_inches='tight')
plt.show()

# Plot the next step probs through time as a line graph
plt.plot(next_step_probs[next_step_probs > 0], color='blue', label='Next Step Probabilities')
plt.plot(rolling_mean_next_step[rolling_mean_next_step > 0], color='red', label='Rolling Mean')
plt.axhline(y=null_prob, color='black', linestyle='--', label='Null Probability')  # null probs
plt.xlabel('Index')
plt.ylabel('Probability')
plt.title('Next Step Probability')
plt.legend()  # Add legend to differentiate lines
plt.savefig(f'{output_dir}/id{buffalo_id}_next_step_probs.png', dpi=300, bbox_inches='tight')
plt.show()

Save the probabilities

We can save the probabilities to a csv file to compare with the SSF probabilities.

# Append the probabilities to the dataframe
buffalo_df['habitat_probs'] = habitat_probs
buffalo_df['move_probs'] = move_probs
buffalo_df['next_step_probs'] = next_step_probs

csv_filename = f'{output_dir}/deepSSF_id{buffalo_id}_n{len(test_data)}_{today_date}.csv'
print(csv_filename)
buffalo_df.to_csv(csv_filename, index=True)

../Python/outputs/model_training/id2005_deepSSF_training_1_2025-07-09/deepSSF_id2005_n10103_2025-07-09.csv