Bulk Density SOC

[lazarusA]

@Meng2077 You can add images (plots) here if you want us to have a look. And in general discuss specific issues if needed.

[Meng2077]

Yesterday, I discussed with Davide to better understand the basics of NN and explore ideas for improving our current model. We agreed that before diving into further fine-tuning or design, we will need to get the model performing reasonably well. At this stage, the main issues are: 1. Overfitting, as shown by the divergence between training and validation loss curves. 2. Underperformance compared to the simpler MLP model I built last year using scikit-learn.

Loss curves from the multivariate NN (similar in the univariate NN) show a rapid drop in validation loss followed by early plateauing, while training loss keeps decreasing — indicating overfitting.

Validation accuracy density plot for SOCconc from the univariate NN (shown in original scale, but metrics calculated in log scale remain poor)

Validation accuracy density plot for SOCconc in log scale from the 2024 scikit-learn MLP

I'll work on reproduce the simple MLP model to get some ideas. Davide suggested a few directions to explore:

1. Use different weight initializer
2. Apply log1p transformation to skewed targets
3. Experiment with a deeper (or simpler?) architecture
4. Try alternative activation functions
5. Tune other model hyperparameters

[BernhardAhrens]

@Meng2077 Last week, I discovered something that may affect the use of the DropOut layers. The states from BatchNorm (for sure) but maybe also from Dropout layers were not correctly passed. Please make sure to merge main back into your branch or I would suggest to maybe not use the dropout for now.

[BernhardAhrens]

@Meng2077 So, here is the commit which gave the decent results for me:
181cee7

In addition to what I did there, make sure to use Dense(16, 1) as the last layer (after z-transforming the outputs, sigmoid in the last layer does not make sense anymore)

NN = Chain(
        Dense(nfeatures, 32, tanh),
        Dense(32, 16, tanh),
        Dense(16, 1)# 
    )

Please play around with it more, I just played around for a bit.

batchsize and learning rates and number and width of hidden layers.

[BernhardAhrens]

I am working on including main into your branch but not quite there.

[BernhardAhrens]

@Meng2077 @lazarusA I also allowed observations to be shuffled at the train validation split. Before we had only time-series in mind. For time-series, you probably want to see how the model performs for a held-out period. For soil profile IDs, it probably does not make sense to keep 70% of the first IDs for training and the rest as validation since there is probably some structure to the order of the IDs. I introduced this option for train:

function train(hybridModel, data, save_ps; nepochs=200, batchsize=10, opt=Adam(0.01), shuffleobs = false)
    # ? split training and validation data
    (x_train, y_train), (x_val, y_val) = splitobs(data; at=0.8, shuffle=shuffleobs)

[Meng2077]

@BernhardAhrens I gave up on "playing" with the parameters, trying to run a simple grid search to get some feeling first. The EasyHybrid.constructNNModel is really handy! @lazarus I saw that prepare_hidden_chain doesn’t seem to handle single hidden layers when it's passed as vectors of ints. Is it by design, or just something not needed yet?

[BernhardAhrens]

@Meng2077 I fixed prepare_hidden_chain - that was an oversight... The fix is in main. We also have a new interactive train/validation loss plot that updates on the fly

[LeRoss92]

With the following setup I get good pure NN & parameter learning results (e.g. predicting LUCAS SOC/BD test NSE; .55/.18):

1. z-score continuous predictors (or minmax in case of ReLU)
2. 1-hot-encoding + linear embedding to about a fifth of #categories of categorical predictors
3. NN:
   • with 1-7 layers
   • each layer 1.0-3.0 times as broad as input vector
   • activation functions: GELU or swish (but ReLU works also)
   • output into sigmoid layer -> already normalized output
     (4. only Parameter Learning: re-minmax to parameter range)
     (5. only Parameter Learning: parameters/forcing into mechanistic model)
     (6. only Parameter Learning: re-minmax predictions)
4. minmax scale target
5. MSE-loss

I usually train with batch size 1000 for 100-200 epochs. On NN depth and width there is hyperparameter training necessary.

Let me know if that works, or if you have any further questions :)

[Meng2077]

Thanks for the tips @LeRoss92 ! I’m running a plain grid search over roughly the parameter ranges you suggested. When you say ‘broad,’ do you mean the number of neurons? With ~360 features, 1–3× that seems high for 5,000 data points.

[LeRoss92]

You're welcome! That means in your case 1.0-3.0 times 360 neurons per layer.

[Meng2077]

According to some experiments so fat, the activation function swish works really well, thanks!

[Meng2077]

@BernhardAhrens Just a minor thing about early stopping. I fetched the current repo this morning and merged it. It looks like the iteration and loss updates now rely on best_agg_loss instead of best_loss, which solves the earlier issue. best_loss no longer seems necessary to me, as everything it did can be handled by best_agg_loss or other existing variables. It might be good to remove it? otherwise it's a bit confusing.

[BernhardAhrens]

Good catch - I updated the train function accordingly. It's in main @Meng2077

[Meng2077]

Somehow the data sparcity handling scheme in the train.jl doesn't work for me, which keeps giving NaN losses after the 1st epoch:

Will work on it now. Maybe it's some mistake from my side..

[lazarusA]

will you mind pushing a new branch (if possible) with the specific setup for this experiment. I agree that we need a better error / nan/ message /handling for these cases. I could take this opportunity to include more checks and intermediate reports. Is weird that the initial loss is fine but then the other not, maybe a single batch is compose of sample where all are nan, the current check only work when there is only one target. I need to add need one to handle the case >1.

[lazarusA]

if the issue in the train function, then the culprit will be this line:

EasyHybrid.jl/src/train.jl

Line 245 in 3ac8e3c

if length(is_no_nan)>0 # ! be careful here, multivariate needs fine tuning

[Meng2077]

The current commit in HybridBulkDensity is the one I am using, the .toml files should be updated.

I checked the input data, the input data shouldn't be the issue. At least one of the 4 targets have valid values in the process, so if length(is_no_nan)>0 or any(.!isnan.(y)) should always be True, as there is no complete empty batch.

I (roughly) looked at the calculation, aggregation and the logging of the losses, and I don't think they are the issues, either. If the first epoch works well, then loss wrappers aren't the issue. I think it's the training process itself, maybe the extremely imbalanced dataset leads to weird gradients.

To check this assumption, I left the parameter search running last night. The idea is see if at least one parameter combination could give some meaningful results. But the process killed itself due to some reasons... Anyway, I will take a further look and keep you updated!

[Meng2077]

Ah, the results in output_tmp shows that some parameter combination did give some meaningful results. So it's gradients update problem, The function train works fine!

[Meng2077]

Preliminary results for Univariate Naive NN, Multivariate Naive NN, and NN living with Sparcity issue (SOCconc, BD and CF aren't complete).

I will continue with getting familiar with Hybrid BD model.

[BernhardAhrens]

Nice! Is that with the same SOC density validation set across models?

[BernhardAhrens]

BTW, @Meng2077 predicted should be in the x-axis and observed in the y-axis. Piñeiro G, Perelman S, Guerschman JP, Paruelo JM (2008) How to evaluate models: observed vs. predicted or predicted vs. observed? Ecological Modelling, 216, 316-322.

[Meng2077]

The validation set so far is not consistent only for data sparsity scenario. Will correct it once all four test scenarios are set. As well as the plot axis. I think they will stay like this for a while.....

[BernhardAhrens]

If I remember correctly, we wrote the Hybrid BD model before I implemented the generic hybrid model definition framework. Check here for an example https://earthyscience.github.io/EasyHybrid.jl/stable/get_started

[BernhardAhrens]

Nice - and you had the right idea with parameter bounds for organic matter bulk density and mineral soil bulk density:

    mBD = out[3, :] .* (1.5f0 - 0.75f0) .+ 0.75f0 # mineral bulk density
    oBD = sigmoid(ps.oBD) .* (0.4f0 - 0.05f0) .+ 0.05f0

This is what I hard coded in the old implementation. The oBD has to be much smaller - I used a range of (0.05, 0.40). The mBD ranges I used were probably too tight - maybe something like (0.75, 2.0)

[Meng2077]

Emm, if the bounds look ok, then we have a bigger issue. The next suspect is the unit.

[BernhardAhrens]

Did you try with the new bounds?

[Meng2077]

Sorry I skipped the latter part after you had the right idea with parameter bounds 😇

But the new bounds indeed do magic! The results look way better, though the BD still look weirdly bounded.

[BernhardAhrens]

Cool! We're getting there 😃.

At the moment, we do:

neural_param_names = [:SOCconc, :oBD]
global_param_names = [:mBD]

I think having mBD as a global parameter maybe the culprit for the bound. It makes more sense to have mBD as a neural_param since we can derive logically that it should depend on texture (how dense can a soil be packed if it clay or sand).

So try:

neural_param_names = [:SOCconc, :mBD]
global_param_names = [:oBD]

That is switch mBD and oBD. For oBD - the theoretical density of organic matter - it might more sense that this does not vary so much between organic soils. However, we could even have this as a neural_param since a peat soil may have a different organic matter bulk density than a organic layer in a deciduous or evergreen needle leave forest.

So, also try neural_param_names = [:SOCconc, :mBD, :oBD].

You can even use different predictors for the different neural_params - for example for oBD it would make sense to use more vegetation derived predictors. See for example here how to use special predictors for different parameters

EasyHybrid.jl/projects/Respiration_Fluxnet/Q10estimation_new.jl

Line 89 in c1a9351

predictors = (Rb = [:SWC_shallow, :P, :WS, :sine_dayofyear, :cos_dayofyear],

.

Maybe something for later 😉