Fortran Parameter Search

When the parameter space is low-dimensional and the model evaluation is cheap, brute-force grid search offers the simplest and most transparent calibration strategy. Fortran's raw loop performance makes it ideal for exhaustive sweeps over millions of parameter combinations.

1. Brute-Force Grid Search

Consider a model with $d$ parameters, each discretised into $N$ grid points. The total number of evaluations is:

$$N_{\text{eval}} = N^d$$

For $d = 3$ parameters with $N = 100$ grid points each, this is $10^6$ evaluations—easily tractable if each evaluation takes microseconds (as in a simple CA step or analytic formula).

At each grid point $\boldsymbol{\theta}_k$, we evaluate an objective (or “figure of merit”):

$$\text{FoM}(\boldsymbol{\theta}) = -\sum_{i=1}^{T}\bigl(N_i^{\text{obs}} - N_i^{\text{sim}}(\boldsymbol{\theta})\bigr)^2$$

The negative sum of squared residuals makes the best fit the maximum of FoM. Alternatively, one can use pattern-match metrics, Kappa statistics for spatial agreement, or fuzzy set similarity for comparing categorical maps.

2. Parallelisation Strategy

Grid search is embarrassingly parallel: each grid point is independent. In Fortran, OpenMP parallelism is added with a single directive:

!$OMP PARALLEL DO PRIVATE(i,j,k,r_val,K_val,sigma_val,fom) &
!$OMP            REDUCTION(MAX:best_fom)
do i = 1, Nr
  do j = 1, NK
    do k = 1, Nsigma
      r_val = r_min + (i-1) * dr
      K_val = K_min + (j-1) * dK
      sigma_val = sigma_min + (k-1) * dsigma
      fom = evaluate_model(r_val, K_val, sigma_val)
      if (fom > best_fom) then
        best_fom = fom
        best_r = r_val; best_K = K_val; best_sigma = sigma_val
      end if
    end do
  end do
end do
!$OMP END PARALLEL DO

The speedup is nearly linear in the number of cores. On a 64-core machine, the $10^6$ evaluations complete in roughly $10^6 / 64 \approx 15{,}000$ sequential evaluations worth of time.

Memory Considerations

Storing all FoM values requires $N^d \times 8$ bytes (double precision). For $100^3 = 10^6$ points, that is only 8 MB. For $100^4$, it is 800 MB—still feasible on modern workstations, but one should consider streaming results to disk for larger searches.

3. Fortran Code: Parameter Search for Urban CA

The following Fortran program performs a 2D grid search over growth rate $r$ and carrying capacity $K$, evaluating a logistic model against synthetic observations. Results are written to a file for visualisation.

program parameter_search
  implicit none
  integer, parameter :: Nr = 80, NK = 80, Nobs = 40
  real(8) :: r_min, r_max, K_min, K_max, dr, dK
  real(8) :: r_val, K_val, N0, t, dt_obs
  real(8) :: N_pred, residual, fom, best_fom
  real(8) :: best_r, best_K
  real(8) :: obs(Nobs), t_obs(Nobs)
  real(8) :: fom_grid(Nr, NK)
  integer :: i, j, n

  ! True parameters for synthetic data
  real(8), parameter :: r_true = 0.05d0, K_true = 500000.d0
  real(8), parameter :: sigma = 15000.d0
  N0 = 50000.d0

  ! Generate synthetic observations
  dt_obs = 2.d0
  do n = 1, Nobs
    t_obs(n) = (n-1) * dt_obs
    obs(n) = K_true / (1.d0 + (K_true/N0 - 1.d0)*exp(-r_true*t_obs(n)))
    ! Add noise (simple LCG for reproducibility)
    obs(n) = obs(n) + sigma * 0.1d0 * sin(dble(n)*13.7d0)
  end do

  ! Grid bounds
  r_min = 0.01d0;  r_max = 0.10d0
  K_min = 200000.d0; K_max = 800000.d0
  dr = (r_max - r_min) / (Nr - 1)
  dK = (K_max - K_min) / (NK - 1)

  best_fom = -1.d30

  ! Nested loop search
  do i = 1, Nr
    r_val = r_min + (i-1)*dr
    do j = 1, NK
      K_val = K_min + (j-1)*dK
      ! Evaluate figure of merit (neg. sum of squared residuals)
      fom = 0.d0
      do n = 1, Nobs
        N_pred = K_val / (1.d0 + (K_val/N0 - 1.d0)*exp(-r_val*t_obs(n)))
        residual = obs(n) - N_pred
        fom = fom - residual*residual
      end do
      fom_grid(i,j) = fom
      if (fom > best_fom) then
        best_fom = fom
        best_r = r_val
        best_K = K_val
      end if
    end do
  end do

  ! Output results
  open(unit=10, file='param_search.dat', status='replace')
  do i = 1, Nr
    do j = 1, NK
      write(10, '(3E16.8)') r_min+(i-1)*dr, K_min+(j-1)*dK, fom_grid(i,j)
    end do
    write(10, *)  ! blank line for gnuplot
  end do
  close(10)

  write(*,'(A,F8.5,A,F12.1)') 'Best: r = ', best_r, '  K = ', best_K
  write(*,'(A,E14.6)') 'Best FoM = ', best_fom
end program parameter_search

4. Python Wrapper: Sensitivity Surface

We replicate the Fortran grid search in Python (for portability) and visualise the objective function landscape. The surface reveals parameter identifiability: a sharp peak means the parameters are well-determined; an elongated ridge signals correlation.

Parameter Sensitivity Surface for Urban Growth Calibration

Python

script.py77 lines

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# --- Synthetic data (same as Fortran version) ---
r_true, K_true, sigma = 0.05, 500_000, 15_000
N0 = 50_000
Nobs = 40
t_obs = np.arange(Nobs) * 2.0

def logistic(t, r, K):
    return K / (1.0 + (K/N0 - 1.0) * np.exp(-r * t))

obs = logistic(t_obs, r_true, K_true) + np.random.normal(0, sigma, Nobs)

# --- Grid search ---
Nr, NK = 80, 80
r_vals = np.linspace(0.01, 0.10, Nr)
K_vals = np.linspace(200_000, 800_000, NK)
fom_grid = np.zeros((Nr, NK))

for i, r in enumerate(r_vals):
    for j, K in enumerate(K_vals):
        pred = logistic(t_obs, r, K)
        fom_grid[i, j] = -np.sum((obs - pred)**2)

# Best fit
best_idx = np.unravel_index(fom_grid.argmax(), fom_grid.shape)
best_r = r_vals[best_idx[0]]
best_K = K_vals[best_idx[1]]

# --- Visualisation ---
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# (a) FoM surface (heatmap)
ax = axes[0]
R, KK = np.meshgrid(r_vals, K_vals, indexing='ij')
# Normalise for visibility
fom_norm = fom_grid - fom_grid.max()
im = ax.contourf(R, KK/1e6, fom_norm, levels=30, cmap='YlGn')
ax.plot(r_true, K_true/1e6, 'r*', markersize=15, label='True')
ax.plot(best_r, best_K/1e6, 'ko', markersize=8, label='Best grid')
ax.set_xlabel('r (growth rate)')
ax.set_ylabel('K (carrying capacity, millions)')
ax.set_title('FoM Surface')
ax.legend(fontsize=9)
plt.colorbar(im, ax=ax, label='FoM (relative)')

# (b) Slice through r at best K
ax = axes[1]
j_best = best_idx[1]
ax.plot(r_vals, fom_grid[:, j_best], color='teal', linewidth=2)
ax.axvline(r_true, color='crimson', linestyle='--', label='True r')
ax.set_xlabel('r')
ax.set_ylabel('FoM')
ax.set_title(f'Slice at K = {K_vals[j_best]/1e6:.2f}M')
ax.legend()

# (c) Slice through K at best r
ax = axes[2]
i_best = best_idx[0]
ax.plot(K_vals/1e6, fom_grid[i_best, :], color='teal', linewidth=2)
ax.axvline(K_true/1e6, color='crimson', linestyle='--', label='True K')
ax.set_xlabel('K (millions)')
ax.set_ylabel('FoM')
ax.set_title(f'Slice at r = {r_vals[i_best]:.4f}')
ax.legend()

plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight')
plt.show()

print(f"Grid search best: r = {best_r:.5f}, K = {best_K:.0f}")
print(f"True values:      r = {r_true:.5f}, K = {K_true:.0f}")
print(f"Grid evaluations: {Nr*NK:,}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

5. Grid Search vs. MCMC

Criterion	Grid Search	MCMC
Scaling	$O(N^d)$ — curse of dimensionality	$O(n)$ per sample, independent of $d$
Output	Full FoM landscape	Posterior samples
Best for	1–3 parameters, cheap models	Any dimensionality, expensive models
Resolution	Fixed by grid spacing	Adaptive (concentrates on high-probability regions)
Parallelism	Embarrassingly parallel	Requires ensemble methods (e.g., parallel tempering)

Practical Recommendation

Use grid search for initial exploration (to understand the landscape shape and identify rough parameter bounds), then refine with MCMC for proper uncertainty quantification. The Fortran grid search is particularly valuable when the model itself is a Fortran simulation: the inner loop stays entirely in compiled code with no Python overhead.

← KPZ Equation Part 8: Fractals & Quantum →

Criterion	Grid Search	MCMC
Scaling	\(O(N^d)\) — curse of dimensionality	\(O(n)\) per sample, independent of \(d\)
Output	Full FoM landscape	Posterior samples
Best for	1–3 parameters, cheap models	Any dimensionality, expensive models
Resolution	Fixed by grid spacing	Adaptive (concentrates on high-probability regions)
Parallelism	Embarrassingly parallel	Requires ensemble methods (e.g., parallel tempering)