Deepseek and Grok LLMs for Fortran coding

[Beliavsky]

Asking Grok 3 to do something related,

Write a function that has an input argument of a 1D array of reals and a format string such as “(f0.6)” and returns a deferred-length character array of sufficient LEN to hold the real numbers converted to strings.

first gave a function where it used only the first element of the array to determine the LEN, with Grok mentioning this limitation. Asked to use a LEN sufficient to hold all the numbers converted to character variables, Grok gave

MODULE string_conversion
  IMPLICIT NONE
  PRIVATE
  PUBLIC :: real_to_string_array  ! Expose only the function to users

CONTAINS

  FUNCTION real_to_string_array(real_array, fmt) RESULT(str_array)
    REAL, INTENT(IN) :: real_array(:)          ! Input 1D array of reals
    CHARACTER(LEN=*), INTENT(IN) :: fmt        ! Format string, e.g., "(f0.6)"
    CHARACTER(LEN=:), ALLOCATABLE :: str_array(:)  ! Deferred-length character array output
    INTEGER :: i, n, len_str, max_len
    CHARACTER(LEN=256) :: temp_str             ! Temporary string for formatting

    ! Get the size of the input array
    n = SIZE(real_array)

    ! Determine the maximum length needed by checking all elements
    max_len = 0
    DO i = 1, n
      WRITE(temp_str, fmt) real_array(i)       ! Format each real number
      len_str = LEN_TRIM(ADJUSTL(temp_str))    ! Get trimmed length of formatted string
      IF (len_str > max_len) max_len = len_str ! Update max length if current is larger
    END DO

    ! Allocate the deferred-length character array with the maximum length
    ALLOCATE(CHARACTER(LEN=max_len) :: str_array(n))

    ! Convert each real number to a string using the provided format
    DO i = 1, n
      WRITE(temp_str, fmt) real_array(i)       ! Write real to temporary string
      str_array(i) = TRIM(ADJUSTL(temp_str))   ! Trim and left-adjust the string
    END DO

  END FUNCTION real_to_string_array

END MODULE string_conversion

! Example program to test the module
PROGRAM test_real_to_string
  USE string_conversion    ! Use the module
  IMPLICIT NONE
  REAL :: numbers(5) = [1.234567, 89.123, -0.005, 12345.6789, -123456.789]
  CHARACTER(LEN=:), ALLOCATABLE :: strings(:)
  INTEGER :: i

  ! Call the function with format "(f0.6)"
  strings = real_to_string_array(numbers, "(f0.6)")

  ! Print the results
  DO i = 1, SIZE(strings)
    PRINT *, "Number ", i, ": '", strings(i), "' (Length: ", LEN(strings(i)), ")"
  END DO

  ! Clean up
  DEALLOCATE(strings)
END PROGRAM test_real_to_string

with output

 Number            1 : '1.234567      ' (Length:           14 )
 Number            2 : '89.123001     ' (Length:           14 )
 Number            3 : '-.005000      ' (Length:           14 )
 Number            4 : '12345.678711  ' (Length:           14 )
 Number            5 : '-123456.789062' (Length:           14 )

I think having coders who are at least mid-level proficient in every programming language, available 24/7, for free, is a boon.

[certik]

@themos

Pick your favourite Fortran-aware LLM and ask it to write a minimal test for calling a function that takes a deferred-length character array with an actual argument that is a deferred-length character array component of a derived type.

Grok can do it: Fortran Function with Deferred-Length Array | Shared Grok Conversation. It made a tiny mistake (forgot to put allocatable on one of the declaration), I just posted a compiler error, and it fixed it. The last printout compiles and runs with both GFortran and LFortran as is:

$ gfortran b.f90 && ./a.out
 Function returned:            3
$ lfortran b.f90
Function returned:     3

[themos]

Grok’s output is in a different league from what I tried, more impressive. But when I see the “each” in “The length of each character string is deferred and set at allocation time” I wonder if it would verify whether “two elements of the same character array can/cannot have different lengths”.

[certik]

Sure, I just posted your comment to Grok: Fortran Function with Deferred-Length Array | Shared Grok Conversation.

Thus, in response to your query: two elements of the same deferred-length character array cannot have different lengths.

I think that’s correct.

[Beliavsky]

Mercury Coder is a new LLM that uses diffusion (used in image generation) rather than autoregression to generate code, as explained in this video. It is extremely fast since it generates the whole code at once, but it did get stuck on the prompt about simulating and fitting a mixture of normals. For the prompt

Define a Fortran derived type for a dataframe with character variable column names and integer indices. Simulate such a dataframe with random uniform data and print statistics on the columns. Declare real variables as real(kind=dp), where dp is a module parmeter set to kind(1.0d0)

it generated a working code once I told it to use fixed-length rather than deferred length character variables for the column names.

module DataFrameModule
    implicit none
    integer, parameter :: dp = kind(1.0d0)
    integer, parameter :: nlen = 10  ! Define the length of column names

    type :: DataFrame
        character(len=nlen), allocatable :: column_names(:)
        real(kind=dp), allocatable :: data(:, :)
        integer :: n_rows, n_cols
    contains
        procedure :: initialize => initialize_dataframe
        procedure :: simulate_data => simulate_dataframe_data
        procedure :: print_stats => print_dataframe_stats
    end type DataFrame

contains

    subroutine initialize_dataframe(this, n_rows, n_cols, column_names)
        class(DataFrame), intent(inout) :: this
        integer, intent(in) :: n_rows, n_cols
        character(len=*), intent(in) :: column_names(n_cols)
        this%n_rows = n_rows
        this%n_cols = n_cols
        allocate(this%column_names(n_cols))
        this%column_names = column_names
        allocate(this%data(n_rows, n_cols))
    end subroutine initialize_dataframe

    subroutine simulate_dataframe_data(this)
        class(DataFrame), intent(inout) :: this
        call random_number(this%data)
    end subroutine simulate_dataframe_data

    subroutine print_dataframe_stats(this)
        class(DataFrame), intent(in) :: this
        real(kind=dp) :: mean, std_dev
        integer :: j
        do j = 1, this%n_cols
            call calculate_mean_std(this%data(:, j), mean, std_dev)
            print *, "Column: ", trim(this%column_names(j))
            print *, "Mean: ", mean
            print *, "Standard Deviation: ", std_dev
            print *
        end do
    end subroutine print_dataframe_stats

    subroutine calculate_mean_std(data, mean, std_dev)
        real(kind=dp), intent(in) :: data(:)
        real(kind=dp), intent(out) :: mean, std_dev
        real(kind=dp) :: sum, sum_sq
        integer :: n, i
        n = size(data)
        sum = 0.0_dp
        sum_sq = 0.0_dp
        do i = 1, n
            sum = sum + data(i)
            sum_sq = sum_sq + data(i)**2
        end do
        mean = sum / n
        std_dev = sqrt((sum_sq / n) - mean**2)
    end subroutine calculate_mean_std

end module DataFrameModule

program Main
    use DataFrameModule
    implicit none
    type(DataFrame) :: df
    character(len=nlen), dimension(3) :: column_names = ["A", "B", "C"]
    integer :: n_rows, n_cols

    n_rows = 100
    n_cols = 3

    call df%initialize(n_rows, n_cols, column_names)
    call df%simulate_data()
    call df%print_stats()

end program Main

giving sample output

 Column: A
 Mean:   0.47517106404222736     
 Standard Deviation:   0.31178191704286989     

 Column: B
 Mean:   0.49751105740387586     
 Standard Deviation:   0.30548776910706921     

 Column: C
 Mean:   0.50007198678270182     
 Standard Deviation:   0.29607418178862621

It is fun to experiment with given its speed, but I think other LLMs are a bit better for more involved coding problems.

[kimala]

@themos

My first experience with one, just now, is that it is incapable of saying “I don’t know”

A colleague (from the informatics department) once told me that he does not use LLMs because “if I need a wrong answer given with extreme confidence, I will just ask an engineer”.

[Beliavsky]

Today’s featured LLM is qwen, described here and the subject of this video. It is an open-source model that is relatively small and which can run on a PC. It presents a lot of “thinking” before it shows code. For the normal mixture simulation and estimation prompt, feeding it gfortran error messages twice gave the following code in 5 minutes.

module mixture_module
  implicit none
  integer, parameter :: dp = kind(1.0d0)

  type :: MixtureModel
    integer :: K
    real(dp), allocatable :: weights(:)
    real(dp), allocatable :: means(:)
    real(dp), allocatable :: sds(:)
  end type MixtureModel

contains

  pure function normal_pdf(x, mu, sigma) result(pdf)
    real(dp), intent(in) :: x, mu, sigma
    real(dp) :: pdf
    real(dp) :: z
    z = (x - mu) / sigma
    pdf = exp(-0.5_dp * z**2) / (sigma * sqrt(2.0_dp * acos(-1.0_dp)))
  end function normal_pdf

  subroutine fit_mixture(data, K, max_iter, tol, model)
    real(dp), intent(in) :: data(:)
    integer, intent(in) :: K, max_iter
    real(dp), intent(in) :: tol
    type(MixtureModel), intent(out) :: model

    real(dp), allocatable :: responsibilities(:,:), new_weights(:), new_means(:), new_sds(:)
    real(dp) :: prev_loglik, current_loglik, sample_mean, sample_sd, r, component_sum
    integer :: iter, n, comp, i, idx
    logical :: converged

    n = size(data)
    model%K = K
    allocate(model%weights(K), model%means(K), model%sds(K))
    
    ! Initialize parameters
    model%weights = 1.0_dp / K
    call random_seed()
    do comp = 1, K
      call random_number(r)
      idx = floor(r * n) + 1
      model%means(comp) = data(idx)
    end do
    sample_mean = sum(data) / n
    sample_sd = sqrt(sum((data - sample_mean)**2) / (n - 1))
    model%sds = sample_sd

    allocate(responsibilities(n, K), new_weights(K), new_means(K), new_sds(K))
    
    prev_loglik = -huge(1.0_dp)
    converged = .false.

    do iter = 1, max_iter
      ! E-step
      do i = 1, n
        component_sum = 0.0_dp
        do comp = 1, K
          responsibilities(i, comp) = model%weights(comp) * normal_pdf(data(i), model%means(comp), model%sds(comp))
          component_sum = component_sum + responsibilities(i, comp)
        end do
        responsibilities(i, :) = responsibilities(i, :) / component_sum
      end do

      ! M-step
      do comp = 1, K
        new_weights(comp) = sum(responsibilities(:, comp)) / n
        new_means(comp) = sum(responsibilities(:, comp) * data) / sum(responsibilities(:, comp))
        new_sds(comp) = sqrt(sum(responsibilities(:, comp) * (data - new_means(comp))**2) / sum(responsibilities(:, comp)))
      end do

      model%weights = new_weights
      model%means = new_means
      model%sds = new_sds

      ! Compute log-likelihood
      current_loglik = 0.0_dp
      do i = 1, n
        component_sum = 0.0_dp
        do comp = 1, K
          component_sum = component_sum + model%weights(comp) * normal_pdf(data(i), model%means(comp), model%sds(comp))
        end do
        current_loglik = current_loglik + log(component_sum)
      end do

      if (abs(current_loglik - prev_loglik) < tol) then
        converged = .true.
        exit
      end if
      prev_loglik = current_loglik
    end do

    deallocate(responsibilities, new_weights, new_means, new_sds)
  end subroutine fit_mixture

  subroutine print_mixture(model)
    type(MixtureModel), intent(in) :: model
    integer :: comp
    print *, "Component | Weight     | Mean       | Std Dev"
    print *, "--------------------------------------------"
    do comp = 1, model%K
      print '(I9, 3F12.6)', comp, model%weights(comp), model%means(comp), model%sds(comp)
    end do
  end subroutine print_mixture

end module mixture_module

program test_mixture
  use mixture_module
  implicit none

  integer, parameter :: K_true = 2
  integer, parameter :: N = 1000
  real(dp) :: data(N)
  type(MixtureModel) :: true_model, estimated_model
  integer :: i
  real(dp) :: u

  ! Set true parameters
  true_model%K = K_true
  allocate(true_model%weights(K_true), true_model%means(K_true), true_model%sds(K_true))
  true_model%weights = [0.3_dp, 0.7_dp]
  true_model%means = [0.0_dp, 5.0_dp]
  true_model%sds = [1.0_dp, 1.0_dp]

  ! Generate data
  call random_seed()
  do i = 1, N
    call random_number(u)
    if (u < true_model%weights(1)) then
      data(i) = true_model%means(1) + true_model%sds(1) * normal_rand()
    else
      data(i) = true_model%means(2) + true_model%sds(2) * normal_rand()
    end if
  end do

  ! Fit the mixture model
  call fit_mixture(data, K_true, 100, 1.0e-6_dp, estimated_model)

  ! Print results
  print *, "True parameters:"
  call print_mixture(true_model)
  print *, "Estimated parameters:"
  call print_mixture(estimated_model)

contains

  function normal_rand() result(z)
    real(dp) :: z, u1, u2, r
    call random_number(u1)
    call random_number(u2)
    r = sqrt(-2.0_dp * log(u1))
    z = r * cos(2.0_dp * acos(-1.0_dp) * u2)
  end function normal_rand

end program test_mixture

which compiles and gives sample output

True parameters:
Component | Weight | Mean | Std Dev

    1    0.300000    0.000000    1.000000
    2    0.700000    5.000000   10.000000

Estimated parameters:
Component | Weight | Mean | Std Dev

    1    0.711504    5.011159    9.874043
    2    0.288496    0.073350    1.129037

[Beliavsky]

Google Gemini and other LLMs can compare two versions of a Fortran source file, describe the differences, and give a commit message. The Python script using Gemini

import subprocess
import sys
import google.generativeai as genai
import os
import textwrap

def print_wrapped_squeeze_blanks(text, width=80):
    # Split into lines and process each
    lines = text.split('\n')
    result = []
    prev_blank = False
    
    for line in lines:
        # Wrap long lines
        wrapped = [line] if len(line) <= width else textwrap.wrap(line, width)
        
        if line.strip():  # Non-blank line
            result.extend(wrapped)
            prev_blank = False
        elif not prev_blank:  # Blank line, but not after another blank
            result.append('')
            prev_blank = True
            
    print('\n'.join(result))

def get_git_diff(filename):
    """Get the diff between current file and last committed version."""
    try:
        # Get the diff using git command
        diff = subprocess.check_output(
            ['git', 'diff', 'HEAD', filename],
            text=True,
            stderr=subprocess.STDOUT
        )
        return diff if diff else "No changes detected in the file."
    except subprocess.CalledProcessError as e:
        return f"Error getting diff: {e.output}"
    except Exception as e:
        return f"Error: {str(e)}"

def generate_commit_message(diff):
    """Generate a commit message using Gemini API."""
    with open(r"c:\python\code\gemini_key.txt", "r") as key_file:
        api_key = key_file.read().strip()
    
    genai.configure(api_key=api_key)
    model_name = "gemini-1.5-flash"
    model = genai.GenerativeModel(model_name)  # Using a fast model
    
    # Create prompt for Gemini
    prompt = (
        "Describe the substantive differences between the local and committed versions of the file"
        f"{diff}. Then provide a one-line commit message, with commit message appearing on its own line."
    )
    
    try:
        response = model.generate_content(prompt)
        return response.text.strip()
    except Exception as e:
        return f"Error generating commit message: {str(e)}"

def main():
    # Check if filename is provided as command line argument
    if len(sys.argv) != 2:
        print("Usage: python script.py <filename>")
        sys.exit(1)
    
    filename = sys.argv[1]
    
    # Verify file exists
    if not os.path.exists(filename):
        print(f"Error: File '{filename}' not found")
        sys.exit(1)
    
    # Get the diff
    diff = get_git_diff(filename)
    
    # Generate and print the commit message
    commit_message = generate_commit_message(diff)
    print_wrapped_squeeze_blanks(commit_message)

if __name__ == "__main__":
    main()

run on a source file basic_stats.f90 gave

The committed version adds functions and subroutines for computing and printing covariance matrices (cov, cov_mat, print_cov_mat). It also modifies the basic_stats module to include covariance-related variables.

Add covariance functions and subroutines

I previously mentioned a simple Groq cloud Fortran coding agent. There are now analogs for Gemini and OpenAI. The latter requires a credit card, and the repo gives an example prompt and the cost ($0.07) of iterating until compilable code was generated.

[aledinola]

Nice example! Is the line

 deallocate(responsibilities, new_weights

at the end of the subroutine fit_mixture really needed in modern fortran? I was wondering if it is there on purpose or if the LLM does not know this feature of F2003, namely that allocated arrays are automatically deallocated when they go out of scope

[Beliavsky]

It is not required, but I have noticed that many LLMs add such code, as do many human Fortran programmers. The LLMs have been trained on more C and C++ code than Fortran code, since more code in the former languages are available, so I think they sometimes do in Fortran what they would do in C and C++. Some LLMs like to intersperse declarations with executable code in Fortran (without using block), as you can do in C and C++.

[Beliavsky]

Now OpenAI has released ChatGPT-5 and removed the other models. For the mixture-of-normals prompt above, after a few iterations it gave the working code

module constants_mod
implicit none
private
public :: dp, twopi, logtwopi

integer, parameter :: dp = kind(1.0d0)
real(kind=dp), parameter :: twopi    = 6.283185307179586476925286766559_dp
real(kind=dp), parameter :: logtwopi = 1.837877066409345339081937709_dp

end module constants_mod

module mixture_mod
use constants_mod, only: dp, twopi, logtwopi
implicit none
private
public :: mixture_params, init_mixture, em_fit_mixture, print_mixture

type mixture_params
  integer :: k = 0
  real(kind=dp), allocatable :: w(:)
  real(kind=dp), allocatable :: mu(:)
  real(kind=dp), allocatable :: sigma(:)
end type mixture_params

contains

pure elemental function normal_logpdf(x, m, s) result(lp)
real(kind=dp), intent(in) :: x, m, s
real(kind=dp) :: lp
lp = -0.5_dp*(logtwopi + 2.0_dp*log(s) + ((x-m)/s)**2)
end function normal_logpdf

pure function logsumexp_vec(v) result(ls)
real(kind=dp), intent(in) :: v(:)
real(kind=dp) :: ls
real(kind=dp) :: vmax
real(kind=dp), allocatable :: tmp(:)
integer :: i, n
n = size(v)
vmax = v(1)
do i = 2, n
  if (v(i) > vmax) then
    vmax = v(i)
  end if
end do
allocate(tmp(n))
do i = 1, n
  tmp(i) = exp(v(i) - vmax)
end do
ls = vmax + log(sum(tmp))
deallocate(tmp)
end function logsumexp_vec

pure elemental function clamp(x, a, b) result(y)
real(kind=dp), intent(in) :: x, a, b
real(kind=dp) :: y
if (x < a) then
  y = a
else if (x > b) then
  y = b
else
  y = x
end if
end function clamp

subroutine init_mixture(y, k, mix)
real(kind=dp), intent(in) :: y(:)
integer, intent(in) :: k
type(mixture_params), intent(out) :: mix
real(kind=dp) :: ymean, ystd
integer :: i, n
real(kind=dp), allocatable :: ys(:), q(:)
n = size(y)
mix%k = k
allocate(mix%w(k), mix%mu(k), mix%sigma(k))
ymean = sum(y)/real(n,dp)
ystd = sqrt(max( (sum((y - ymean)**2)/real(max(n-1,1),dp)), 1.0e-6_dp ))
mix%w = 1.0_dp/real(k,dp)
allocate(ys(n), q(k))
ys = y
call sort_inplace(ys)
do i = 1, k
  q(i) = percentile_of_sorted(ys, (real(i,dp)-0.5_dp)/real(k,dp))
end do
mix%mu = q
mix%sigma = ystd
deallocate(ys, q)
end subroutine init_mixture

subroutine em_fit_mixture(y, mix, max_iter, tol, min_sigma, seed, converged, n_iter)
real(kind=dp), intent(in) :: y(:)
type(mixture_params), intent(inout) :: mix
integer, intent(in), optional :: max_iter
real(kind=dp), intent(in), optional :: tol
real(kind=dp), intent(in), optional :: min_sigma
integer, intent(in), optional :: seed
logical, intent(out), optional :: converged
integer, intent(out), optional :: n_iter
integer :: it, itmax, k, n, i, j
real(kind=dp) :: eps, minsig
real(kind=dp) :: ll, ll_prev, diff
real(kind=dp), allocatable :: r(:,:), logw(:), logpdf(:,:), logpost(:,:), denom(:)
real(kind=dp), allocatable :: nk(:), mu_new(:), sigma_new(:), w_new(:)
real(kind=dp) :: num, den, s2
itmax = merge(max_iter, 500, present(max_iter))
eps   = merge(tol, 1.0e-8_dp, present(tol))
minsig = merge(min_sigma, 1.0e-6_dp, present(min_sigma))
if (present(seed)) then
  continue
end if
k = mix%k
n = size(y)
allocate(r(n,k), logw(k), logpdf(n,k), logpost(n,k), denom(n))
allocate(nk(k), mu_new(k), sigma_new(k), w_new(k))
ll_prev = -huge(1.0_dp)

do it = 1, itmax
  do j = 1, k
    logw(j) = log(clamp(mix%w(j), 1.0e-300_dp, 1.0e+300_dp))
  end do
  do j = 1, k
    do i = 1, n
      logpdf(i,j) = normal_logpdf(y(i), mix%mu(j), mix%sigma(j))
    end do
  end do
  do i = 1, n
    do j = 1, k
      logpost(i,j) = logw(j) + logpdf(i,j)
    end do
    denom(i) = logsumexp_vec(logpost(i,:))
    do j = 1, k
      r(i,j) = exp(logpost(i,j) - denom(i))
    end do
  end do
  ll = sum(denom)

  do j = 1, k
    nk(j) = sum(r(:,j))
  end do
  do j = 1, k
    if (nk(j) <= 0.0_dp) nk(j) = 1.0e-12_dp
  end do

  do j = 1, k
    num = 0.0_dp
    do i = 1, n
      num = num + r(i,j)*y(i)
    end do
    den = nk(j)
    mu_new(j) = num / den
  end do

  do j = 1, k
    s2 = 0.0_dp
    do i = 1, n
      s2 = s2 + r(i,j)*(y(i) - mu_new(j))**2
    end do
    s2 = s2 / nk(j)
    sigma_new(j) = sqrt(max(s2, minsig))
  end do

  w_new = nk / real(n,dp)
  do j = 1, k
    sigma_new(j) = max(sigma_new(j), minsig)
    w_new(j) = max(w_new(j), 1.0e-12_dp)
  end do
  w_new = w_new / sum(w_new)

  diff = abs(ll - ll_prev)
  mix%w = w_new
  mix%mu = mu_new
  mix%sigma = sigma_new
  if (diff < eps) exit
  ll_prev = ll
end do

if (present(converged)) then
  converged = (it < itmax)
end if
if (present(n_iter)) n_iter = it

deallocate(r, logw, logpdf, logpost, denom, nk, mu_new, sigma_new, w_new)
end subroutine em_fit_mixture

subroutine print_mixture(mix, label)
type(mixture_params), intent(in) :: mix
character(*), intent(in) :: label
integer :: j
write(*,'(a)') repeat("-", 58)
write(*,'(a)') trim(label)
write(*,'(a)') repeat("-", 58)
write(*,'(a)') "   j           weight              mean             sigma"
write(*,'(a)') repeat("-", 58)
do j = 1, mix%k
  write(*,'(i4,3(2x,es15.7))') j, mix%w(j), mix%mu(j), mix%sigma(j)
end do
write(*,'(a)') repeat("-", 58)
end subroutine print_mixture

subroutine sort_inplace(x)
real(kind=dp), intent(inout) :: x(:)
integer :: i, j, n
real(kind=dp) :: key
n = size(x)
do i = 2, n
  key = x(i)
  j = i - 1
  do
    if (j < 1) exit
    if (x(j) > key) then
      x(j+1) = x(j)
      j = j - 1
    else
      exit
    end if
  end do
  x(j+1) = key
end do
end subroutine sort_inplace

pure function percentile_of_sorted(xs, p) result(q)
real(kind=dp), intent(in) :: xs(:)
real(kind=dp), intent(in) :: p
real(kind=dp) :: q
integer :: n
real(kind=dp) :: pp, idx, frac
integer :: i0, i1
n = size(xs)
pp = clamp(p, 0.0_dp, 1.0_dp)
idx = 1.0_dp + pp*real(max(n-1,0),dp)
i0 = int(floor(idx))
i1 = int(ceiling(idx))
if (i0 < 1) i0 = 1
if (i1 < 1) i1 = 1
if (i0 > n) i0 = n
if (i1 > n) i1 = n
frac = idx - real(i0,dp)
q = (1.0_dp - frac)*xs(i0) + frac*xs(i1)
end function percentile_of_sorted

end module mixture_mod

program test_mixture_em
use constants_mod, only: dp
use mixture_mod
implicit none

integer :: n, k_true, k_fit, i, j, iters
type(mixture_params) :: mix_true, mix_est
real(kind=dp), allocatable :: y(:)
real(kind=dp), allocatable :: u(:), z(:)
real(kind=dp), allocatable :: cumw(:)
integer, allocatable :: comp(:)
real(kind=dp) :: tol, minsig
logical :: ok
integer :: nseed, clock
integer, allocatable :: seed_put(:)

n = 5000
k_true = 3
k_fit = 3
allocate(y(n), u(n), z(n), comp(n))
allocate(mix_true%w(k_true), mix_true%mu(k_true), mix_true%sigma(k_true))
mix_true%k = k_true
mix_true%w = [0.55_dp, 0.35_dp, 0.10_dp]
mix_true%mu = [-1.0_dp, 1.5_dp, 4.0_dp]
mix_true%sigma = [0.7_dp, 0.6_dp, 1.2_dp]

! portable seeding: query seed length and fill from system_clock
call random_seed(size=nseed)
allocate(seed_put(nseed))
call system_clock(count=clock)
do i = 1, nseed
  seed_put(i) = clock + 37*i
end do
call random_seed(put=seed_put)
deallocate(seed_put)

allocate(cumw(k_true))
cumw(1) = mix_true%w(1)
do j = 2, k_true
  cumw(j) = cumw(j-1) + mix_true%w(j)
end do
call random_number(u)
call randn_array(z)

do i = 1, n
  if (u(i) < cumw(1)) then
    comp(i) = 1
  else if (u(i) < cumw(2)) then
    comp(i) = 2
  else
    comp(i) = 3
  end if
  j = comp(i)
  y(i) = mix_true%mu(j) + mix_true%sigma(j)*z(i)
end do

call print_mixture(mix_true, "true mixture parameters")

call init_mixture(y, k_fit, mix_est)
tol = 1.0e-8_dp
minsig = 1.0e-6_dp
call em_fit_mixture(y, mix_est, max_iter=1000, tol=tol, min_sigma=minsig, converged=ok, n_iter=iters)

call print_mixture(mix_est, "estimated mixture parameters (em)")

contains

subroutine randn_array(z)
use constants_mod, only: dp, twopi
real(kind=dp), intent(out) :: z(:)
integer :: m, i
real(kind=dp), allocatable :: u1(:), u2(:)
real(kind=dp) :: r, theta
m = size(z)
allocate(u1((m+1)/2), u2((m+1)/2))
call random_number(u1)
call random_number(u2)
do i = 1, size(u1)
  u1(i) = max(u1(i), 1.0e-12_dp)
  r = sqrt(-2.0_dp*log(u1(i)))
  theta = twopi*u2(i)
  z(2*i-1) = r*cos(theta)
  if (2*i <= m) then
    z(2*i) = r*sin(theta)
  end if
end do
deallocate(u1, u2)
end subroutine randn_array

end program test_mixture_em

with sample output

----------------------------------------------------------
true mixture parameters
----------------------------------------------------------
   j           weight              mean             sigma
----------------------------------------------------------
   1    5.5000000E-01   -1.0000000E+00    7.0000000E-01
   2    3.5000000E-01    1.5000000E+00    6.0000000E-01
   3    1.0000000E-01    4.0000000E+00    1.2000000E+00
----------------------------------------------------------
----------------------------------------------------------
estimated mixture parameters (em)
----------------------------------------------------------
   j           weight              mean             sigma
----------------------------------------------------------
   1    5.5237690E-01   -1.0101226E+00    6.8874118E-01
   2    3.5657542E-01    1.5156374E+00    6.2589453E-01
   3    9.1047677E-02    4.0244089E+00    1.1813303E+00
----------------------------------------------------------

I can’t say yet whether GPT-5 is better than o3.

[msz59]

@Beliavsky

Now OpenAI has released ChatGPT-5 and removed the other models.

It defaults to 5 but in Legacy models menu I can still choose 4o. NB. I have the Plus subscription.