Ruby Fiddle - Function behaves differently between C and Ruby

Question

I am using Ruby Fiddle to access a C function to perform some heavy calculations. The C function works perfectly well when called directly, but when used via Fiddle it returns various rows of nans and infs in an unpredictable way. The function operates on matrices that are passed as pointers to arrays.

I have debugged the C code and everything works fine. I have also saved the various parameters passed to the C function to file to be sure that Fiddle was not passing some weird values, but there is no obvious (at least for me) problem.

Additionally it seems that for 'smaller' matrices this does not happen. Apologies in advance for the code being very long, but this is the only way to exactly reproduce what it is happening. Data for testing is in this file. (gist). You can just copy and paste the Ruby and C test data to the code below.

If it helps, I am working on macos Catalina and Ruby 2.2.4 (requirement).

The C code is the following:

#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

double *mrt(
    double *person_total_shortwave,
    int rows_pts, int cols_pts,
    double *dry_bulb_temperatures,
    double *ground_temperatures,
    double *atmospheric_longwave,
    double *sky_view_factors,
    double *shading_view_factors_matrix,
    int rows_svf, int cols_svf,
    double *shading_temperatures_matrix,
    int rows_st, int cols_st,
    int size_dbt,
    double person_emissivity_shortwave, double person_emissivity_longwave,
    double surroundings_emissivity, double ground_emissivity, double ground_person_view_factor);

int save_to_file(char *filename, int m, int n, double *mat)
{
    FILE *fp;
    int i, j;

    if ((fp = freopen(filename, "w", stdout)) == NULL)
    {
        printf("Cannot open file.\n");
        exit(1);
    }
    for (i = 0; i < m; i++)
        for (j = 0; j < n; j++)
            if (j == n - 1)
                printf("%.17g \n", mat[i * n + j]);
            else
                printf("%.17g ", mat[i * n + j]);
    fclose(fp);
    return 0;
}

int main(void)
{

    // REPLACE WITH C DATA FROM FILE

    double *mrt_result;
    mrt_result = mrt(person_total_shortwave[0],
                     rows_pts, cols_pts,
                     dry_bulb_temperatures,
                     ground_temperatures,
                     atmospheric_longwave,
                     sky_view_factors,
                     shading_view_factors_matrix[0],
                     rows_svf, cols_svf,
                     shading_temperatures_matrix[0],
                     rows_st, cols_st,
                     size_dbt,
                     person_emissivity_shortwave, person_emissivity_longwave,
                     surroundings_emissivity, ground_emissivity, ground_person_view_factor);

    // save_to_file(rows_pts, cols_pts, mrt_result);

    return 0;
}
// https://www.dropbox.com/s/ix06edrrctad421/Calculation%20of%20Mean%20Radiant%20Temperature3.odt?dl=0

double *mrt(
    double *person_total_shortwave,
    int rows_pts, int cols_pts,
    double *dry_bulb_temperatures,
    double *ground_temperatures,
    double *atmospheric_longwave,
    double *sky_view_factors,
    double *shading_view_factors_matrix,
    int rows_svf, int cols_svf,
    double *shading_temperatures_matrix,
    int rows_st, int cols_st,
    int size_dbt,
    double person_emissivity_shortwave, double person_emissivity_longwave,
    double surroundings_emissivity, double ground_emissivity, double ground_person_view_factor)
{
    save_to_file("PTS.txt", rows_pts, cols_pts, person_total_shortwave);
    save_to_file("SVF.txt", 1, cols_svf, sky_view_factors);
    save_to_file("SVSM.txt", rows_svf, cols_svf, shading_view_factors_matrix);
    save_to_file("STM.txt", rows_st, cols_st, shading_temperatures_matrix);

    double *mrt_mat = (double *)calloc(rows_pts * cols_pts, sizeof(double));
    save_to_file("MRT_MAT0.txt", rows_pts, cols_pts, mrt_mat);

    int t, c, k;
    double sigma = 5.67E-8;
    double body_area = 1.94;

    double tmrt4_shortwave, tmrt4_longwave_ground, tmrt4_atm_longwave, tmrt4_surroundings;
    double tmrt4_shading[rows_svf];
    memset(tmrt4_shading, 0.0, rows_svf * sizeof(double));
    double tmrt4;
    double surroundings_view_factor;

    // t runs through the timesteps
    // c runs through the points in the mesh
    for (t = 0; t < rows_pts; t++)
    {
        for (c = 0; c < cols_pts; c++)
        {
            tmrt4_shortwave = 1.0 / (sigma * body_area) * (person_emissivity_shortwave / person_emissivity_longwave) * person_total_shortwave[t * cols_pts + c];

            // We are assuming that the ground is at ambient temperature
            tmrt4_longwave_ground = ground_person_view_factor * ground_emissivity * pow((273.15 + ground_temperatures[t]), 4);

            // Here we are using the actual long wave radiation from the sky
            tmrt4_atm_longwave = (1.0 - ground_person_view_factor) / sigma * sky_view_factors[c] * atmospheric_longwave[t];

            // We need to remove the contribution of all the shading devices
            // k runs through the shading devices
            surroundings_view_factor = 1.0 - sky_view_factors[c];
            for (k = 0; k < rows_svf; k++)
            {
                surroundings_view_factor -= shading_view_factors_matrix[k * cols_svf + c];
            }
            tmrt4_surroundings = (1.0 - ground_person_view_factor) * surroundings_view_factor * surroundings_emissivity * pow((273.15 + dry_bulb_temperatures[t] - 0.0), 4);

            // We now need to account for all contributions of the shading devices
            for (k = 0; k < rows_svf; k++)
            {
                tmrt4_shading[k] = (1.0 - ground_person_view_factor) * (shading_view_factors_matrix[k * cols_svf + c]) * surroundings_emissivity * pow((273.15 + shading_temperatures_matrix[k * cols_svf + t]), 4);
            }

            // Finally we add them all (see paper) for the total contribution
            tmrt4 = tmrt4_shortwave + tmrt4_longwave_ground + tmrt4_atm_longwave + tmrt4_surroundings;
            for (k = 0; k < rows_svf; k++)
            {
                tmrt4 += tmrt4_shading[k];
            }

            // Just convert to celsius
            mrt_mat[t * cols_pts + c] = pow(tmrt4, 0.25) - 273.15;
        }
    }
    save_to_file("MRT_MAT.txt", rows_pts, cols_pts, mrt_mat);

    // double x = 1.5;
    // while (1)
    // {
    //     x *= sin(x) / atan(x) * tanh(x) * sqrt(x);
    // }
    return mrt_mat;
}

and I compile is with clang -g --extra-warnings utils.c -o utils.

The Ruby code is the following

require "fiddle"
require "fiddle/import"

module RG
  extend Fiddle::Importer
  @handler.handlers.each { |h| h.close unless h.close_enabled? } unless @handler.nil?
  GC.start

  dlload File.join("utils")
  extern "double* mrt(double*, int, int, double*, double*, double*, double*, double*, int, int, double*, int, int, int, double, double, double, double, double)"

  def self.mat_to_ptr(matrix)
    return Fiddle::Pointer[matrix.flatten.pack("E*")]
  end

  def self.ptr_to_mat(ptr, rows, cols)
    length = rows * cols * Fiddle::SIZEOF_DOUBLE
    mat = ptr[0, length]
    return mat.unpack("E*").each_slice(cols).to_a
  end

  def self.mean_radiant_temperature(
    person_total_shortwave,
    dry_bulb_temperatures,
    ground_temperatures,
    atmospheric_longwave,
    sky_view_factors,
    shading_view_factors_matrix,
    shading_temperatures_matrix,
    person_emissivity_shortwave,
    person_emissivity_longwave,
    surroundings_emissivity,
    ground_emissivity,
    ground_person_view_factor
  )
    rows_pts = person_total_shortwave.size
    cols_pts = person_total_shortwave[0].size
    person_total_shortwave_pointer = RG.mat_to_ptr(person_total_shortwave)

    dry_bulb_temperatures_pointer = RG.mat_to_ptr(dry_bulb_temperatures)
    ground_temperatures_pointer = RG.mat_to_ptr(ground_temperatures)
    size_dbt = dry_bulb_temperatures.size

    atmospheric_longwave_pointer = RG.mat_to_ptr(atmospheric_longwave)
    sky_view_factors_pointer = RG.mat_to_ptr(sky_view_factors)

    rows_svf = shading_view_factors_matrix.size
    if rows_svf > 0
      cols_svf = shading_view_factors_matrix[0].size
    else
      cols_svf = 0
    end

    shading_view_factors_matrix_pointer = RG.mat_to_ptr(shading_view_factors_matrix)

    rows_st = shading_temperatures_matrix.size
    if rows_st > 0
      cols_st = shading_temperatures_matrix[0].size
    else
      cols_st = 0
    end

    shading_temperatures_matrix_pointer = RG.mat_to_ptr(shading_temperatures_matrix)

    mrt_pointer = mrt(
      person_total_shortwave_pointer,
      rows_pts, cols_pts,
      dry_bulb_temperatures_pointer,
      ground_temperatures_pointer,
      atmospheric_longwave_pointer,
      sky_view_factors_pointer,
      shading_view_factors_matrix_pointer,
      rows_svf, cols_svf,
      shading_temperatures_matrix_pointer,
      rows_st, cols_st,
      size_dbt,
      person_emissivity_shortwave,
      person_emissivity_longwave,
      surroundings_emissivity,
      ground_emissivity,
      ground_person_view_factor
    )
    return RG.ptr_to_mat(mrt_pointer, rows_pts, cols_pts)
  end
end

// REPLACE WITH RUBY DATA FROM FILE

mean_radiant_temperatures = RG.mean_radiant_temperature(
  person_total_shortwave,
  dry_bulb_temperatures,
  ground_temperatures,
  atmospheric_longwave,
  sky_view_factors,
  shading_view_factors_matrix,
  shading_temperatures_matrix,
  person_emissivity_shortwave,
  person_emissivity_longwave,
  surroundings_emissivity,
  ground_emissivity,
  ground_person_view_factor
)
File.open("MRT_MAT_RUBY.txt", "w") do |f|
  mean_radiant_temperatures.each do |row|
    f.puts row.join(' ')
  end
end

If you want to test it, first launch ./utils. It will save a few files on the local folder. Have a look at MRT_MAT.txt.

Now launch the ruby code several times. It will generate the same file, but you will notice that "often" the file will contain random rows with nan and inf. The same data is returned to Ruby and saved in the file MRT_MAT_RUBY.txt in the local directory.

I am quite familiar with Ruby, but C is not really my strength. It would be great being able to debug the C code when called from Ruby, but I don't really know how to do it.

Answer 1

I guess the string created by pack method in mat_to_ptr method is collected by GC.

You should keep the string until pte_to_mat is called.