<commit_before>/*
 * Copyright (C) 2015-2018 S[&]T, The Netherlands.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * 3. Neither the name of the copyright holder nor the names of its
 *    contributors may be used to endorse or promote products derived from
 *    this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */

#include "harp-internal.h"

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

#define MAX_NAME_LENGTH 128
#define LATLON_BLOCK_SIZE 1024

typedef enum binning_type_enum
{
    binning_skip,
    binning_remove,
    binning_average,
    binning_sum,
    binning_angle,
    binning_time_min,
    binning_time_max,
    binning_datetime
} binning_type;


static binning_type get_binning_type(harp_variable *variable)
{
    int i;

    /* variables with enumeration values get removed */
    if (variable->num_enum_values > 0)
    {
        return binning_remove;
    }

    /* any variable with a time dimension that is not the first dimension gets removed */
    for (i = 1; i < variable->num_dimensions; i++)
    {
        if (variable->dimension[i] == harp_dimension_time)
        {
            return binning_remove;
        }
    }

    /* we only bin variables with a time dimension */
    if (variable->num_dimensions == 0 || variable->dimension_type[0] != harp_dimension_time)
    {
        return binning_skip;
    }

    /* we can't bin string values */
    if (variable->data_type == harp_type_string)
    {
        return binning_remove;
    }

    /* we can't bin values that have no unit */
    if (variable->unit == NULL)
    {
        long variable_name_length = strlen(variable->name);

        /* '...count' variables are just summed up, but only if they are unitless and use an int32 data type */
        if (variable->data_type == harp_type_int32 && variable_name_length >= 5 &&
            strcmp(&variable->name[variable_name_length - 5], "count") == 0)
        {
            return binning_sum;
        }

        return binning_remove;
    }

    /* uncertainty propagation needs to be handled differently (remove for now) */
    if (strstr(variable->name, "_uncertainty") != NULL)
    {
        return binning_remove;
    }

    /* we can't bin averaging kernels */
    if (strstr(variable->name, "_avk") != NULL)
    {
        return binning_remove;
    }

    /* we can't bin latitude/longitude bounds if they define an area */
    if (strcmp(variable->name, "latitude_bounds") == 0 || strcmp(variable->name, "longitude_bounds") == 0)
    {
        if (variable->num_dimensions > 0 &&
            variable->dimension_type[variable->num_dimensions - 1] == harp_dimension_independent &&
            variable->dimension[variable->num_dimensions - 1] > 2)
        {
            return binning_remove;
        }
    }

    /* note that latitude does not require special angular treatment since there is no wrap-around */
    if (strstr(variable->name, "longitude") != NULL || strstr(variable->name, "angle") != NULL ||
        strstr(variable->name, "direction") != NULL)
    {
        return binning_angle;
    }

    /* use minimum/maximum for datetime start/stop */
    if (variable->num_dimensions == 1)
    {
        if (strcmp(variable->name, "datetime_start") == 0)
        {
            return binning_time_min;
        }
        if (strcmp(variable->name, "datetime_stop") == 0)
        {
            return binning_time_max;
        }
    }

    /* use average by default */
    return binning_average;
}

static binning_type get_spatial_binning_type(harp_variable *variable, int area_binning)
{
    binning_type type = get_binning_type(variable);

    if (type != binning_remove)
    {
        /* remove all latitude/longitude variables */
        if (strstr(variable->name, "latitude") != NULL || strstr(variable->name, "longitude") != NULL)
        {
            return binning_remove;
        }
        if (area_binning && strcmp(variable->name, "count") == 0)
        {
            return binning_remove;
        }
        if (strcmp(variable->name, "datetime_length") == 0)
        {
            return binning_remove;
        }
        if (strcmp(variable->name, "datetime") == 0)
        {
            return binning_datetime;
        }
    }

    return type;
}

/* find a <variable->name>_count variable.
 * If the variable exists but is invalid its entry in the bintype array will be set to binning_remove.
 */
static int get_count_for_variable(harp_product *product, harp_variable *variable, binning_type *bintype, int32_t *count)
{
    char variable_name[MAX_NAME_LENGTH];
    harp_variable *count_variable = NULL;
    int index = -1;
    long i, j;

    snprintf(variable_name, MAX_NAME_LENGTH, "%s_count", variable->name);
    if (harp_product_has_variable(product, variable_name))
    {
        if (harp_product_get_variable_index_by_name(product, variable_name, &index) != 0)
        {
            return -1;
        }
        if (bintype[index] == binning_sum)
        {
            if (product->variable[index]->data_type != harp_type_int32 ||
                product->variable[index]->num_dimensions != variable->num_dimensions)
            {
                bintype[index] = binning_remove;
            }
            else
            {
                count_variable = product->variable[index];
            }
        }
    }

    if (count_variable == NULL && harp_product_has_variable(product, "count"))
    {
        if (harp_product_get_variable_index_by_name(product, "count", &index) != 0)
        {
            return -1;
        }
        if (bintype[index] == binning_sum)
        {
            if (product->variable[index]->data_type != harp_type_int32 ||
                product->variable[index]->num_dimensions > variable->num_dimensions)
            {
                bintype[index] = binning_remove;
            }
            else
            {
                count_variable = product->variable[index];
            }
        }
    }

    if (count_variable == NULL)
    {
        return 0;
    }

    /* make sure that the dimensions of the count variable match the first dimensions of the given variable */
    for (i = 0; i < count_variable->num_dimensions; i++)
    {
        if (count_variable->dimension_type[i] != variable->dimension_type[i] ||
            count_variable->dimension[i] != variable->dimension[i])
        {
            bintype[index] = binning_remove;
            return 0;
        }
    }

    /* store data into count parameter */
    if (variable->num_elements == count_variable->num_elements)
    {
        memcpy(count, count_variable->data.int32_data, count_variable->num_elements * sizeof(int32_t));
    }
    else
    {
        long num_sub_elements = variable->num_elements / count_variable->num_elements;

        for (i = 0; i < count_variable->num_elements; i++)
        {
            for (j = 0; j < num_sub_elements; j++)
            {
                count[i * num_sub_elements + j] = count_variable->data.int32_data[i];
            }
        }
    }

    return 1;
}

static int add_count_variable(harp_product *product, binning_type *bintype, const char *variable_name,
                              int num_dimensions, harp_dimension_type *dimension_type, long *dimension, int32_t *count)
{
    char count_variable_name[MAX_NAME_LENGTH];
    harp_variable *variable;

    if (variable_name != NULL)
    {
        snprintf(count_variable_name, MAX_NAME_LENGTH, "%s_count", variable_name);
    }
    else
    {
        strcpy(count_variable_name, "count");
    }

    if (!harp_product_has_variable(product, count_variable_name))
    {
        if (harp_variable_new(count_variable_name, harp_type_int32, num_dimensions, dimension_type, dimension,
                              &variable) != 0)
        {
            return -1;
        }
        memcpy(variable->data.int32_data, count, variable->num_elements * sizeof(int32_t));
        if (harp_product_add_variable(product, variable) != 0)
        {
            harp_variable_delete(variable);
            return -1;
        }
        bintype[product->num_variables - 1] = binning_sum;
    }
    else
    {
        int index;

        if (harp_product_get_variable_index_by_name(product, count_variable_name, &index) != 0)
        {
            return -1;
        }
        if (bintype[index] == binning_remove)
        {
            /* if the existing count variable was scheduled for removal than replace it with a new one */
            if (harp_variable_new(count_variable_name, harp_type_int32, num_dimensions, dimension_type, dimension,
                                  &variable) != 0)
            {
                return -1;
            }
            memcpy(variable->data.int32_data, count, variable->num_elements * sizeof(int32_t));
            if (harp_product_replace_variable(product, variable) != 0)
            {
                harp_variable_delete(variable);
                return -1;
            }
            bintype[index] = binning_sum;
        }
    }

    return 0;
}

static int filter_spatial_binable_variables(harp_product *product, int area_binning)
{
    int i;

    for (i = product->num_variables - 1; i >= 0; i--)
    {
        if (get_spatial_binning_type(product->variable[i], area_binning) == binning_remove)
        {
            if (harp_product_remove_variable(product, product->variable[i]) != 0)
            {
                return -1;
            }
        }
    }

    return 0;
}

static int find_matching_cells_and_weights_for_bounds(harp_variable *latitude_bounds, harp_variable *longitude_bounds,
                                                      long num_latitude_edges, double *latitude_edges,
                                                      long num_longitude_edges, double *longitude_edges,
                                                      long *num_latlon_index, long **latlon_cell_index,
                                                      double **latlon_weight)
{
    long num_elements;
    long i;

    num_elements = latitude_bounds->dimension[0];
    for (i = 0; i < num_elements; i++)
    {
        num_latlon_index[i] = 0;
    }

    return 0;
}

static int find_matching_cells_for_points(harp_variable *latitude, harp_variable *longitude, long num_latitude_edges,
                                          double *latitude_edges, long num_longitude_edges, double *longitude_edges,
                                          long *num_latlon_index, long **latlon_cell_index)
{
    long latitude_index = -1;
    long longitude_index = -1;
    long cumsum_index = 0;
    long num_elements;
    long i;

    num_elements = latitude->dimension[0];
    for (i = 0; i < num_elements; i++)
    {
        harp_interpolate_find_index(num_latitude_edges, latitude_edges, latitude->data.double_data[i], &latitude_index);
        if (latitude_index < 0 || latitude_index >= num_latitude_edges)
        {
            num_latlon_index[i] = 0;
            continue;
        }
        harp_interpolate_find_index(num_longitude_edges, longitude_edges, longitude->data.double_data[i],
                                    &longitude_index);
        if (longitude_index < 0 || longitude_index >= num_longitude_edges)
        {
            num_latlon_index[i] = 0;
            continue;
        }
        num_latlon_index[i] = 1;
        if (cumsum_index % LATLON_BLOCK_SIZE == 0)
        {
            long *new_latlon_cell_index;

            new_latlon_cell_index = realloc(*latlon_cell_index, (cumsum_index + LATLON_BLOCK_SIZE) * sizeof(long));
            if (new_latlon_cell_index == NULL)
            {
                harp_set_error(HARP_ERROR_OUT_OF_MEMORY, "out of memory (could not allocate %lu bytes) (%s:%u)",
                               (cumsum_index + LATLON_BLOCK_SIZE) * sizeof(long), __FILE__, __LINE__);
                return -1;
            }
            *latlon_cell_index = new_latlon_cell_index;
        }
        (*latlon_cell_index)[cumsum_index] = latitude_index * (num_longitude_edges - 1) + longitude_index;
        cumsum_index++;
    }

    return 0;
}

/** \addtogroup harp_product
 * @{
 */

/** Bin the product's variables.
 * This will bin all variables in the time dimension. Each time sample will be put in the bin defined by bin_index.
 * All variables with a time dimension will then be resampled using these bins.
 * The resulting value for each variable will be the average of all values for the bin.
 * Variables with multiple dimensions will have all elements in the sub dimensions averaged on an element by element
 * basis.
 *
 * Variables that have a time dimension but no unit (or using a string data type) will be removed.
 *
 * All variables that are binned (except existing 'count' variables) are converted to a double data type.
 * Bins that have no samples will end up with a NaN value.
 *
 * If the product did not already have a 'count' variable then a 'count' variable will be added to the product that
 * will contain the number of samples per bin.
 *
 * \param product Product to regrid.
 * \param num_bins Number of target bins.
 * \param num_elements Length of bin_index array (should equal the length of the time dimension)
 * \param bin_index Array of target bin index numbers (0 .. num_bins-1) for each sample in the time dimension.
 *
 * \return
 *   \arg \c 0, Success.
 *   \arg \c -1, Error occurred (check #harp_errno).
 */
LIBHARP_API int harp_product_bin(harp_product *product, long num_bins, long num_elements, long *bin_index)
{
    harp_dimension_type dimension_type[HARP_MAX_NUM_DIMS];
    binning_type *bintype = NULL;
    long filtered_count_size = 0;
    int32_t *filtered_count = NULL;
    int32_t *count = NULL;
    long *index = NULL;
    long i, j, k;

    if (num_elements != product->dimension[harp_dimension_time])
    {
        harp_set_error(HARP_ERROR_INVALID_ARGUMENT, "num_elements (%ld) does not match time dimension length (%ld) "
                       "(%s:%u)", num_elements, product->dimension[harp_dimension_time], __FILE__, __LINE__);
        return -1;
    }

    for (i = 0; i < num_elements; i++)
    {
        if (bin_index[i] < 0 || bin_index[i] >= num_bins)
        {
            harp_set_error(HARP_ERROR_INVALID_ARGUMENT, "bin_index[%ld] (%ld) should be in the range [0..%ld) (%s:%u)",
                           i, bin_index[i], num_bins, __FILE__, __LINE__);
            return -1;
        }
    }

    /* make 'bintype' big enough to also store any count variables that we may want to add (i.e. 1 + factor 2) */
    bintype = malloc((2 * product->num_variables + 1) * sizeof(binning_type));
    if (bintype == NULL)
    {
        harp_set_error(HARP_ERROR_OUT_OF_MEMORY, "out of memory (could not allocate %lu bytes) (%s:%u)",
                       (2 * product->num_variables + 1) * sizeof(binning_type), __FILE__, __LINE__);
        goto error;
    }
    for (k = 0; k < product->num_variables; k++)
    {
        bintype[k] = get_binning_type(product->variable[k]);

        /* determine the maximum number of elements (as size for the 'filtered_count' array) */
        if (bintype[k] != binning_remove && bintype[k] != binning_skip)
        {
            long total_num_elements = product->variable[k]->num_elements;

            if (num_bins > num_elements)
            {
                /* use longest time dimension (before vs. after binning) */
                total_num_elements = num_bins * (total_num_elements / num_elements);
            }
            if (total_num_elements > filtered_count_size)
            {
                filtered_count_size = total_num_elements;
            }
        }
    }

    index = malloc(num_bins * sizeof(long));
    if (index == NULL)
    {
        harp_set_error(HARP_ERROR_OUT_OF_MEMORY, "out of memory (could not allocate %lu bytes) (%s:%u)",
                       num_bins * sizeof(long), __FILE__, __LINE__);
        goto error;
    }
    count = malloc(num_bins * sizeof(int32_t));
    if (count == NULL)
    {
        harp_set_error(HARP_ERROR_OUT_OF_MEMORY, "out of memory (could not allocate %lu bytes) (%s:%u)",
                       num_bins * sizeof(int32_t), __FILE__, __LINE__);
        goto error;
    }
    filtered_count = malloc(filtered_count_size * sizeof(int32_t));
    if (filtered_count == NULL)
    {
        harp_set_error(HARP_ERROR_OUT_OF_MEMORY, "out of memory (could not allocate %lu bytes) (%s:%u)",
                       filtered_count_size * sizeof(int32_t), __FILE__, __LINE__);
        goto error;
    }

    /* for each bin, store the index of the first sample that contributes to the bin */
    for (i = 0; i < num_bins; i++)
    {
        index[i] = 0;
        count[i] = 0;
    }
    for (i = 0; i < num_elements; i++)
    {
        if (count[bin_index[i]] == 0)
        {
            index[bin_index[i]] = i;
        }
        count[bin_index[i]]++;
    }

    /* pre-process all variables */
    for (k = 0; k < product->num_variables; k++)
    {
        harp_variable *variable;
        long num_sub_elements;

        if (bintype[k] == binning_skip || bintype[k] == binning_remove)
        {
            continue;
        }

        variable = product->variable[k];
        num_sub_elements = variable->num_elements / num_elements;

        /* convert variables to double */
        if (bintype[k] != binning_sum)
        {
            if (harp_variable_convert_data_type(variable, harp_type_double) != 0)
            {
                goto error;
            }
        }

        /* convert all angles to degrees */
        if (bintype[k] == binning_angle)
        {
            if (harp_convert_unit(variable->unit, "degrees", variable->num_elements, variable->data.double_data) != 0)
            {
                goto error;
            }
            /* wrap all angles to within [x-180,x+180] where x is the first sample of the bin */
            for (i = 0; i < num_elements; i++)
            {
                long target_index = index[bin_index[i]];

                if (target_index != i)
                {
                    for (j = 0; j < num_sub_elements; j++)
                    {
                        double min = variable->data.double_data[target_index * num_sub_elements + j] - 180;
                        double max = min + 360;

                        variable->data.double_data[i * num_sub_elements + j] =
                            harp_wrap(variable->data.double_data[i * num_sub_elements + j], min, max);
                    }
                }
            }
        }

        /* pre-multiply variables by existing counts */
        if (bintype[k] == binning_average || bintype[k] == binning_angle)
        {
            int result;

            result = get_count_for_variable(product, variable, bintype, filtered_count);
            if (result < 0)
            {
                goto error;
            }
            if (result == 1)
            {
                /* multiply by the count */
                for (i = 0; i < variable->num_elements; i++)
                {
                    variable->data.double_data[i] *= filtered_count[i];
                }
            }
        }
    }

    /* sum up all samples into bins (in place) and create count variables where needed */
    for (k = 0; k < product->num_variables; k++)
    {
        harp_variable *variable;
        long num_sub_elements;

        if (bintype[k] == binning_skip || bintype[k] == binning_remove)
        {
            continue;
        }

        variable = product->variable[k];
        assert(variable->dimension[0] == num_elements);
        num_sub_elements = variable->num_elements / num_elements;

        /* TODO: use special approach for lat/lon (bounds) averaging (use spatial averaging + bounds) */

        if (bintype[k] == binning_time_min)
        {
            /* take minimum of all values per bin */
            assert(num_sub_elements == 1);
            for (i = 0; i < num_elements; i++)
            {
                long target_index = index[bin_index[i]];

                if (variable->data.double_data[i] < variable->data.double_data[target_index])
                {
                    variable->data.double_data[target_index] = variable->data.double_data[i];
                }
            }
        }
        else if (bintype[k] == binning_time_max)
        {
            /* take maximum of all values per bin */
            assert(num_sub_elements == 1);
            for (i = 0; i < num_elements; i++)
            {
                long target_index = index[bin_index[i]];

                if (variable->data.double_data[i] > variable->data.double_data[target_index])
                {
                    variable->data.double_data[target_index] = variable->data.double_data[i];
                }
            }
        }
        else
        {
            int store_count_variable = 0;

            /* sum up all values of a bin into the location of the first sample */

            if (variable->data_type != harp_type_int32)
            {
                memset(filtered_count, 0, variable->num_elements * sizeof(int32_t));
            }

            for (i = 0; i < num_elements; i++)
            {
                long target_index = index[bin_index[i]];

                if (target_index != i)
                {
                    if (variable->data_type == harp_type_int32)
                    {
                        for (j = 0; j < num_sub_elements; j++)
                        {
                            variable->data.int32_data[target_index * num_sub_elements + j] +=
                                variable->data.int32_data[i * num_sub_elements + j];
                        }
                    }
                    else
                    {
                        for (j = 0; j < num_sub_elements; j++)
                        {
                            if (harp_isnan(variable->data.double_data[i * num_sub_elements + j]))
                            {
                                filtered_count[i * num_sub_elements + j] = 0;
                                store_count_variable = 1;
                            }
                            else
                            {
                                filtered_count[i * num_sub_elements + j] = 1;
                                variable->data.double_data[target_index * num_sub_elements + j] +=
                                    variable->data.double_data[i * num_sub_elements + j];
                            }
                        }
                    }
                }
                else if (variable->data_type != harp_type_int32)
                {
                    for (j = 0; j < num_sub_elements; j++)
                    {
                        if (harp_isnan(variable->data.double_data[target_index * num_sub_elements + j]))
                        {
                            filtered_count[target_index * num_sub_elements + j] = 0;
                            variable->data.double_data[target_index * num_sub_elements + j] = 0;
                            store_count_variable = 1;
                        }
                        else
                        {
                            filtered_count[target_index * num_sub_elements + j] = 1;
                        }
                    }
                }
            }

            if (store_count_variable)
            {
                if (add_count_variable(product, bintype, variable->name, variable->num_dimensions,
                                       variable->dimension_type, variable->dimension, filtered_count) != 0)
                {
                    goto error;
                }
            }
        }
    }

    /* resample variables */
    for (k = 0; k < product->num_variables; k++)
    {
        if (bintype[k] == binning_skip || bintype[k] == binning_remove)
        {
            continue;
        }

        /* resample the time dimension to the target bins */
        if (harp_variable_rearrange_dimension(product->variable[k], 0, num_bins, index) != 0)
        {
            goto error;
        }
    }

    /* update product dimensions */
    product->dimension[harp_dimension_time] = num_bins;

    /* add global count variable if it didn't exist yet */
    dimension_type[0] = harp_dimension_time;
    if (add_count_variable(product, bintype, NULL, 1, dimension_type, &num_bins, count) != 0)
    {
        goto error;
    }

    /* post-process all variables */
    for (k = 0; k < product->num_variables; k++)
    {
        harp_variable *variable;
        long num_sub_elements;
        int count_applied = 0;

        if (bintype[k] == binning_skip || bintype[k] == binning_remove)
        {
            continue;
        }

        variable = product->variable[k];
        num_sub_elements = variable->num_elements / num_elements;

        /* divide variables by the sample count */
        if (bintype[k] == binning_average || bintype[k] == binning_angle)
        {
            double nan_value = harp_nan();
            int result;

            result = get_count_for_variable(product, variable, bintype, filtered_count);
            if (result < 0)
            {
                goto error;
            }
            if (result == 1)
            {
                /* divide by the count (or set value to NaN if count==0) */
                for (i = 0; i < variable->num_elements; i++)
                {
                    assert(i < filtered_count_size);
                    if (filtered_count[i] == 0)
                    {
                        variable->data.double_data[i] = nan_value;
                    }
                    else if (filtered_count[i] > 1)
                    {
                        variable->data.double_data[i] /= filtered_count[i];
                    }
                }
                count_applied = 1;
            }
        }

        /* set all empty bins to NaN (for double) or 0 (for int32) */
        if (!count_applied)
        {
            for (i = 0; i < num_bins; i++)
            {
                if (count[i] == 0)
                {
                    double nan_value = harp_nan();

                    if (variable->data_type == harp_type_int32)
                    {
                        for (j = 0; j < num_sub_elements; j++)
                        {
                            variable->data.int32_data[i * num_sub_elements + j] = 0;
                        }
                    }
                    else
                    {
                        for (j = 0; j < num_sub_elements; j++)
                        {
                            variable->data.double_data[i * num_sub_elements + j] = nan_value;
                        }
                    }
                }
            }
        }

        /* convert angle variables back to original unit */
        if (bintype[k] == binning_angle)
        {
            /* convert all angles back to the original unit */
            if (harp_convert_unit("degrees", variable->unit, variable->num_elements, variable->data.double_data) != 0)
            {
                goto error;
            }
        }
    }

    /* remove all variables that need to be removed (in reverse order!) */
    for (k = product->num_variables - 1; k >= 0; k--)
    {
        if (bintype[k] == binning_remove)
        {
            if (harp_product_remove_variable(product, product->variable[k]) != 0)
            {<commit_msg>Don't use existing counts for averaging angles<commit_after>/*
 * Copyright (C) 2015-2018 S[&]T, The Netherlands.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * 3. Neither the name of the copyright holder nor the names of its
 *    contributors may be used to endorse or promote products derived from
 *    this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */

#include "harp-internal.h"

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

#define MAX_NAME_LENGTH 128
#define LATLON_BLOCK_SIZE 1024

typedef enum binning_type_enum
{
    binning_skip,
    binning_remove,
    binning_average,
    binning_sum,
    binning_angle,
    binning_time_min,
    binning_time_max,
    binning_datetime
} binning_type;


static binning_type get_binning_type(harp_variable *variable)
{
    int i;

    /* variables with enumeration values get removed */
    if (variable->num_enum_values > 0)
    {
        return binning_remove;
    }

    /* any variable with a time dimension that is not the first dimension gets removed */
    for (i = 1; i < variable->num_dimensions; i++)
    {
        if (variable->dimension[i] == harp_dimension_time)
        {
            return binning_remove;
        }
    }

    /* we only bin variables with a time dimension */
    if (variable->num_dimensions == 0 || variable->dimension_type[0] != harp_dimension_time)
    {
        return binning_skip;
    }

    /* we can't bin string values */
    if (variable->data_type == harp_type_string)
    {
        return binning_remove;
    }

    /* we can't bin values that have no unit */
    if (variable->unit == NULL)
    {
        long variable_name_length = strlen(variable->name);

        /* '...count' variables are just summed up, but only if they are unitless and use an int32 data type */
        if (variable->data_type == harp_type_int32 && variable_name_length >= 5 &&
            strcmp(&variable->name[variable_name_length - 5], "count") == 0)
        {
            if (strstr(variable->name, "latitude") != NULL || strstr(variable->name, "longitude") != NULL ||
                strstr(variable->name, "angle") != NULL || strstr(variable->name, "direction") != NULL)
            {
                /* we can't propagate average of angles (since we would also need the 'magnitudes' for this) */
                /* so just remove any counts for angles */
                return binning_remove;
            }
            return binning_sum;
        }

        return binning_remove;
    }

    /* uncertainty propagation needs to be handled differently (remove for now) */
    if (strstr(variable->name, "_uncertainty") != NULL)
    {
        return binning_remove;
    }

    /* we can't bin averaging kernels */
    if (strstr(variable->name, "_avk") != NULL)
    {
        return binning_remove;
    }

    /* we can't bin latitude/longitude bounds if they define an area */
    if (strcmp(variable->name, "latitude_bounds") == 0 || strcmp(variable->name, "longitude_bounds") == 0)
    {
        if (variable->num_dimensions > 0 &&
            variable->dimension_type[variable->num_dimensions - 1] == harp_dimension_independent &&
            variable->dimension[variable->num_dimensions - 1] > 2)
        {
            return binning_remove;
        }
    }

    if (strstr(variable->name, "latitude") != NULL || strstr(variable->name, "longitude") != NULL ||
        strstr(variable->name, "angle") != NULL || strstr(variable->name, "direction") != NULL)
    {
        return binning_angle;
    }

    /* use minimum/maximum for datetime start/stop */
    if (variable->num_dimensions == 1)
    {
        if (strcmp(variable->name, "datetime_start") == 0)
        {
            return binning_time_min;
        }
        if (strcmp(variable->name, "datetime_stop") == 0)
        {
            return binning_time_max;
        }
    }

    /* use average by default */
    return binning_average;
}

static binning_type get_spatial_binning_type(harp_variable *variable, int area_binning)
{
    binning_type type = get_binning_type(variable);

    if (type != binning_remove)
    {
        /* remove all latitude/longitude variables */
        if (strstr(variable->name, "latitude") != NULL || strstr(variable->name, "longitude") != NULL)
        {
            return binning_remove;
        }
        if (area_binning && strcmp(variable->name, "count") == 0)
        {
            return binning_remove;
        }
        if (strcmp(variable->name, "datetime_length") == 0)
        {
            return binning_remove;
        }
        if (strcmp(variable->name, "datetime") == 0)
        {
            return binning_datetime;
        }
    }

    return type;
}

static int get_count_variable_for_variable(harp_product *product, harp_variable *variable, binning_type *bintype,
                                           harp_variable **count_variable)
{
    char variable_name[MAX_NAME_LENGTH];
    int index;
    int i;

    *count_variable = NULL;

    snprintf(variable_name, MAX_NAME_LENGTH, "%s_count", variable->name);
    if (!harp_product_has_variable(product, variable_name))
    {
        return 0;
    }

    if (harp_product_get_variable_index_by_name(product, variable_name, &index) != 0)
    {
        return -1;
    }
    if (bintype[index] != binning_sum)
    {
        return 0;
    }
    if (product->variable[index]->data_type != harp_type_int32 ||
        product->variable[index]->num_dimensions != variable->num_dimensions)
    {
        bintype[index] = binning_remove;
        return 0;
    }

    /* make sure that the dimensions of the count variable match the dimensions of the given variable */
    for (i = 0; i < product->variable[index]->num_dimensions; i++)
    {
        if (product->variable[index]->dimension_type[i] != variable->dimension_type[i] ||
            product->variable[index]->dimension[i] != variable->dimension[i])
        {
            bintype[index] = binning_remove;
            return 0;
        }
    }

    *count_variable = product->variable[index];

    return 0;
}

/* find a <variable->name>_count variable.
 * If the variable exists but is invalid its entry in the bintype array will be set to binning_remove.
 */
static int get_count_for_variable(harp_product *product, harp_variable *variable, binning_type *bintype, int32_t *count)
{
    harp_variable *count_variable = NULL;
    long i, j;

    if (get_count_variable_for_variable(product, variable, bintype, &count_variable) != 0)
    {
        return -1;
    }

    if (count_variable == NULL && harp_product_has_variable(product, "count"))
    {
        int index;

        if (harp_product_get_variable_index_by_name(product, "count", &index) != 0)
        {
            return -1;
        }
        if (bintype[index] == binning_sum)
        {
            if (product->variable[index]->data_type != harp_type_int32 ||
                product->variable[index]->num_dimensions > variable->num_dimensions)
            {
                bintype[index] = binning_remove;
            }
            else
            {
                count_variable = product->variable[index];

                /* make sure that the dimensions of the count variable match the first dimensions of the variable */
                for (i = 0; i < count_variable->num_dimensions; i++)
                {
                    if (count_variable->dimension_type[i] != variable->dimension_type[i] ||
                        count_variable->dimension[i] != variable->dimension[i])
                    {
                        bintype[index] = binning_remove;
                        return 0;
                    }
                }
            }
        }

    }

    if (count_variable == NULL)
    {
        return 0;
    }

    /* store data into count parameter */
    if (variable->num_elements == count_variable->num_elements)
    {
        memcpy(count, count_variable->data.int32_data, count_variable->num_elements * sizeof(int32_t));
    }
    else
    {
        long num_sub_elements = variable->num_elements / count_variable->num_elements;

        for (i = 0; i < count_variable->num_elements; i++)
        {
            for (j = 0; j < num_sub_elements; j++)
            {
                count[i * num_sub_elements + j] = count_variable->data.int32_data[i];
            }
        }
    }

    return 1;
}

static int add_count_variable(harp_product *product, binning_type *bintype, const char *variable_name,
                              int num_dimensions, harp_dimension_type *dimension_type, long *dimension, int32_t *count)
{
    char count_variable_name[MAX_NAME_LENGTH];
    harp_variable *variable;
    int index;

    if (variable_name != NULL)
    {
        snprintf(count_variable_name, MAX_NAME_LENGTH, "%s_count", variable_name);
    }
    else
    {
        strcpy(count_variable_name, "count");
    }

    if (!harp_product_has_variable(product, count_variable_name))
    {
        if (harp_variable_new(count_variable_name, harp_type_int32, num_dimensions, dimension_type, dimension,
                              &variable) != 0)
        {
            return -1;
        }
        memcpy(variable->data.int32_data, count, variable->num_elements * sizeof(int32_t));
        if (harp_product_add_variable(product, variable) != 0)
        {
            harp_variable_delete(variable);
            return -1;
        }
        bintype[product->num_variables - 1] = binning_sum;
    }
    else
    {
        if (harp_product_get_variable_index_by_name(product, count_variable_name, &index) != 0)
        {
            return -1;
        }
        if (bintype[index] == binning_remove)
        {
            /* if the existing count variable was scheduled for removal than replace it with a new one */
            if (harp_variable_new(count_variable_name, harp_type_int32, num_dimensions, dimension_type, dimension,
                                  &variable) != 0)
            {
                return -1;
            }
            memcpy(variable->data.int32_data, count, variable->num_elements * sizeof(int32_t));
            if (harp_product_replace_variable(product, variable) != 0)
            {
                harp_variable_delete(variable);
                return -1;
            }
            bintype[index] = binning_sum;
        }
        assert(bintype[index] = binning_sum);
        /* if the count variable already exists and does not get removed then we assume it is correct/consistent
         * (i.e. existing count=0 <-> variable=NaN) */
    }

    return 0;
}

static int filter_spatial_binable_variables(harp_product *product, int area_binning)
{
    int i;

    for (i = product->num_variables - 1; i >= 0; i--)
    {
        if (get_spatial_binning_type(product->variable[i], area_binning) == binning_remove)
        {
            if (harp_product_remove_variable(product, product->variable[i]) != 0)
            {
                return -1;
            }
        }
    }

    return 0;
}

static int find_matching_cells_and_weights_for_bounds(harp_variable *latitude_bounds, harp_variable *longitude_bounds,
                                                      long num_latitude_edges, double *latitude_edges,
                                                      long num_longitude_edges, double *longitude_edges,
                                                      long *num_latlon_index, long **latlon_cell_index,
                                                      double **latlon_weight)
{
    long num_elements;
    long i;

    num_elements = latitude_bounds->dimension[0];
    for (i = 0; i < num_elements; i++)
    {
        num_latlon_index[i] = 0;
    }

    return 0;
}

static int find_matching_cells_for_points(harp_variable *latitude, harp_variable *longitude, long num_latitude_edges,
                                          double *latitude_edges, long num_longitude_edges, double *longitude_edges,
                                          long *num_latlon_index, long **latlon_cell_index)
{
    long latitude_index = -1;
    long longitude_index = -1;
    long cumsum_index = 0;
    long num_elements;
    long i;

    num_elements = latitude->dimension[0];
    for (i = 0; i < num_elements; i++)
    {
        harp_interpolate_find_index(num_latitude_edges, latitude_edges, latitude->data.double_data[i], &latitude_index);
        if (latitude_index < 0 || latitude_index >= num_latitude_edges)
        {
            num_latlon_index[i] = 0;
            continue;
        }
        harp_interpolate_find_index(num_longitude_edges, longitude_edges, longitude->data.double_data[i],
                                    &longitude_index);
        if (longitude_index < 0 || longitude_index >= num_longitude_edges)
        {
            num_latlon_index[i] = 0;
            continue;
        }
        num_latlon_index[i] = 1;
        if (cumsum_index % LATLON_BLOCK_SIZE == 0)
        {
            long *new_latlon_cell_index;

            new_latlon_cell_index = realloc(*latlon_cell_index, (cumsum_index + LATLON_BLOCK_SIZE) * sizeof(long));
            if (new_latlon_cell_index == NULL)
            {
                harp_set_error(HARP_ERROR_OUT_OF_MEMORY, "out of memory (could not allocate %lu bytes) (%s:%u)",
                               (cumsum_index + LATLON_BLOCK_SIZE) * sizeof(long), __FILE__, __LINE__);
                return -1;
            }
            *latlon_cell_index = new_latlon_cell_index;
        }
        (*latlon_cell_index)[cumsum_index] = latitude_index * (num_longitude_edges - 1) + longitude_index;
        cumsum_index++;
    }

    return 0;
}

/** \addtogroup harp_product
 * @{
 */

/** Bin the product's variables.
 * This will bin all variables in the time dimension. Each time sample will be put in the bin defined by bin_index.
 * All variables with a time dimension will then be resampled using these bins.
 * The resulting value for each variable will be the average of all values for the bin.
 * Variables with multiple dimensions will have all elements in the sub dimensions averaged on an element by element
 * basis.
 *
 * Variables that have a time dimension but no unit (or using a string data type) will be removed.
 *
 * All variables that are binned (except existing 'count' variables) are converted to a double data type.
 * Bins that have no samples will end up with a NaN value.
 *
 * If the product did not already have a 'count' variable then a 'count' variable will be added to the product that
 * will contain the number of samples per bin.
 *
 * \param product Product to regrid.
 * \param num_bins Number of target bins.
 * \param num_elements Length of bin_index array (should equal the length of the time dimension)
 * \param bin_index Array of target bin index numbers (0 .. num_bins-1) for each sample in the time dimension.
 *
 * \return
 *   \arg \c 0, Success.
 *   \arg \c -1, Error occurred (check #harp_errno).
 */
LIBHARP_API int harp_product_bin(harp_product *product, long num_bins, long num_elements, long *bin_index)
{
    harp_dimension_type dimension_type[HARP_MAX_NUM_DIMS];
    binning_type *bintype = NULL;
    long filtered_count_size = 0;
    int32_t *filtered_count = NULL;
    int32_t *count = NULL;
    long *index = NULL;
    long i, j, k;

    if (num_elements != product->dimension[harp_dimension_time])
    {
        harp_set_error(HARP_ERROR_INVALID_ARGUMENT, "num_elements (%ld) does not match time dimension length (%ld) "
                       "(%s:%u)", num_elements, product->dimension[harp_dimension_time], __FILE__, __LINE__);
        return -1;
    }

    for (i = 0; i < num_elements; i++)
    {
        if (bin_index[i] < 0 || bin_index[i] >= num_bins)
        {
            harp_set_error(HARP_ERROR_INVALID_ARGUMENT, "bin_index[%ld] (%ld) should be in the range [0..%ld) (%s:%u)",
                           i, bin_index[i], num_bins, __FILE__, __LINE__);
            return -1;
        }
    }

    /* make 'bintype' big enough to also store any count variables that we may want to add (i.e. 1 + factor 2) */
    bintype = malloc((2 * product->num_variables + 1) * sizeof(binning_type));
    if (bintype == NULL)
    {
        harp_set_error(HARP_ERROR_OUT_OF_MEMORY, "out of memory (could not allocate %lu bytes) (%s:%u)",
                       (2 * product->num_variables + 1) * sizeof(binning_type), __FILE__, __LINE__);
        goto error;
    }
    for (k = 0; k < product->num_variables; k++)
    {
        bintype[k] = get_binning_type(product->variable[k]);

        /* determine the maximum number of elements (as size for the 'filtered_count' array) */
        if (bintype[k] != binning_remove && bintype[k] != binning_skip)
        {
            long total_num_elements = product->variable[k]->num_elements;

            if (num_bins > num_elements)
            {
                /* use longest time dimension (before vs. after binning) */
                total_num_elements = num_bins * (total_num_elements / num_elements);
            }
            if (total_num_elements > filtered_count_size)
            {
                filtered_count_size = total_num_elements;
            }
        }
    }

    index = malloc(num_bins * sizeof(long));
    if (index == NULL)
    {
        harp_set_error(HARP_ERROR_OUT_OF_MEMORY, "out of memory (could not allocate %lu bytes) (%s:%u)",
                       num_bins * sizeof(long), __FILE__, __LINE__);
        goto error;
    }
    count = malloc(num_bins * sizeof(int32_t));
    if (count == NULL)
    {
        harp_set_error(HARP_ERROR_OUT_OF_MEMORY, "out of memory (could not allocate %lu bytes) (%s:%u)",
                       num_bins * sizeof(int32_t), __FILE__, __LINE__);
        goto error;
    }
    filtered_count = malloc(filtered_count_size * sizeof(int32_t));
    if (filtered_count == NULL)
    {
        harp_set_error(HARP_ERROR_OUT_OF_MEMORY, "out of memory (could not allocate %lu bytes) (%s:%u)",
                       filtered_count_size * sizeof(int32_t), __FILE__, __LINE__);
        goto error;
    }

    /* for each bin, store the index of the first sample that contributes to the bin */
    for (i = 0; i < num_bins; i++)
    {
        index[i] = 0;
        count[i] = 0;
    }
    for (i = 0; i < num_elements; i++)
    {
        if (count[bin_index[i]] == 0)
        {
            index[bin_index[i]] = i;
        }
        count[bin_index[i]]++;
    }

    /* pre-process all variables */
    for (k = 0; k < product->num_variables; k++)
    {
        harp_variable *variable;
        long num_sub_elements;

        if (bintype[k] == binning_skip || bintype[k] == binning_remove)
        {
            continue;
        }

        variable = product->variable[k];
        num_sub_elements = variable->num_elements / num_elements;

        /* convert variables to double */
        if (bintype[k] != binning_sum)
        {
            if (harp_variable_convert_data_type(variable, harp_type_double) != 0)
            {
                goto error;
            }
        }

        if (bintype[k] == binning_angle)
        {
            /* convert all angles to complex values [cos(x),sin(x)] */
            if (harp_convert_unit(variable->unit, "rad", variable->num_elements, variable->data.double_data) != 0)
            {
                goto error;
            }
            if (harp_variable_add_dimension(variable, variable->num_dimensions, harp_dimension_independent, 2) != 0)
            {
                goto error;
            }
            for (i = 0; i < variable->num_elements; i += 2)
            {
                variable->data.double_data[i] = cos(variable->data.double_data[i]);
                variable->data.double_data[i + 1] = sin(variable->data.double_data[i]);
            }
        }

        /* pre-multiply variables by existing counts */
        if (bintype[k] == binning_average)
        {
            int result;

            result = get_count_for_variable(product, variable, bintype, filtered_count);
            if (result < 0)
            {
                goto error;
            }
            if (result == 1)
            {
                /* multiply by the count */
                for (i = 0; i < variable->num_elements; i++)
                {
                    variable->data.double_data[i] *= filtered_count[i];
                }
            }
        }
    }

    /* sum up all samples into bins (in place) and create count variables where needed */
    for (k = 0; k < product->num_variables; k++)
    {
        harp_variable *variable;
        long num_sub_elements;

        if (bintype[k] == binning_skip || bintype[k] == binning_remove)
        {
            continue;
        }

        variable = product->variable[k];
        assert(variable->dimension[0] == num_elements);
        num_sub_elements = variable->num_elements / num_elements;

        if (bintype[k] == binning_time_min)
        {
            /* take minimum of all values per bin */
            assert(num_sub_elements == 1);
            for (i = 0; i < num_elements; i++)
            {
                long target_index = index[bin_index[i]];

                if (variable->data.double_data[i] < variable->data.double_data[target_index])
                {
                    variable->data.double_data[target_index] = variable->data.double_data[i];
                }
            }
        }
        else if (bintype[k] == binning_time_max)
        {
            /* take maximum of all values per bin */
            assert(num_sub_elements == 1);
            for (i = 0; i < num_elements; i++)
            {
                long target_index = index[bin_index[i]];

                if (variable->data.double_data[i] > variable->data.double_data[target_index])
                {
                    variable->data.double_data[target_index] = variable->data.double_data[i];
                }
            }
        }
        else
        {
            int store_count_variable = 0;

            /* sum up all values of a bin into the location of the first sample */

            if (variable->data_type != harp_type_int32)
            {
                memset(filtered_count, 0, variable->num_elements * sizeof(int32_t));
            }

            for (i = 0; i < num_elements; i++)
            {
                long target_index = index[bin_index[i]];

                if (target_index != i)
                {
                    if (bintype[k] == binning_angle)
                    {
                        /* for angle variables we use 1 filtered_count element per complex pair */
                        for (j = 0; j < num_sub_elements; j += 2)
                        {
                            if (harp_isnan(variable->data.double_data[i * num_sub_elements + j]))
                            {
                                filtered_count[(i * num_sub_elements + j) / 2] = 0;
                                store_count_variable = 1;
                            }
                            else
                            {
                                filtered_count[(i * num_sub_elements + j) / 2] = 1;
                                variable->data.double_data[target_index * num_sub_elements + j] +=
                                    variable->data.double_data[i * num_sub_elements + j];
                                variable->data.double_data[target_index * num_sub_elements + j + 1] +=
                                    variable->data.double_data[i * num_sub_elements + j + 1];
                            }
                        }
                    }
                    else if (variable->data_type == harp_type_int32)
                    {
                        for (j = 0; j < num_sub_elements; j++)
                        {
                            variable->data.int32_data[target_index * num_sub_elements + j] +=
                                variable->data.int32_data[i * num_sub_elements + j];
                        }
                    }
                    else
                    {
                        for (j = 0; j < num_sub_elements; j++)
                        {
                            if (harp_isnan(variable->data.double_data[i * num_sub_elements + j]))
                            {
                                filtered_count[i * num_sub_elements + j] = 0;
                                store_count_variable = 1;
                            }
                            else
                            {
                                filtered_count[i * num_sub_elements + j] = 1;
                                variable->data.double_data[target_index * num_sub_elements + j] +=
                                    variable->data.double_data[i * num_sub_elements + j];
                            }
                        }
                    }
                }
                else if (variable->data_type != harp_type_int32)
                {
                    if (bintype[k] == binning_angle)
                    {
                        for (j = 0; j < num_sub_elements; j += 2)
                        {
                            if (harp_isnan(variable->data.double_data[target_index * num_sub_elements + j]))
                            {
                                filtered_count[(target_index * num_sub_elements + j) / 2] = 0;
                                variable->data.double_data[target_index * num_sub_elements + j] = 0;
                                variable->data.double_data[target_index * num_sub_elements + j + 1] = 0;
                                store_count_variable = 1;
                            }
                            else
                            {
                                filtered_count[(target_index * num_sub_elements + j) / 2] = 1;
                            }
                        }
                    }
                    else
                    {
                        for (j = 0; j < num_sub_elements; j++)
                        {
                            if (harp_isnan(variable->data.double_data[target_index * num_sub_elements + j]))
                            {
                                filtered_count[target_index * num_sub_elements + j] = 0;
                                variable->data.double_data[target_index * num_sub_elements + j] = 0;
                                store_count_variable = 1;
                            }
                            else
                            {
                                filtered_count[target_index * num_sub_elements + j] = 1;
                            }
                        }
                    }
                }
            }
            if (store_count_variable)
            {
                if (bintype[k] == binning_angle)
                {
                    /* we store the count for angles (temporarily) to be able to set the count=0 values to NaN */
                    /* don't include the 'complex' dimension for the count variable for angles */
                    if (add_count_variable(product, bintype, variable->name, variable->num_dimensions - 1,
                                           variable->dimension_type, variable->dimension, filtered_count) != 0)
                    {
                        goto error;
                    }
                }
                else
                {
                    if (add_count_variable(product, bintype, variable->name, variable->num_dimensions,
                                           variable->dimension_type, variable->dimension, filtered_count) != 0)
                    {
                        goto error;
                    }
                }
            }
        }
    }

    /* resample variables */
    for (k = 0; k < product->num_variables; k++)
    {
        if (bintype[k] == binning_skip || bintype[k] == binning_remove)
        {
            continue;
        }

        /* resample the time dimension to the target bins */
        if (harp_variable_rearrange_dimension(product->variable[k], 0, num_bins, index) != 0)
        {
            goto error;
        }
    }

    /* update product dimensions */
    product->dimension[harp_dimension_time] = num_bins;

    /* add global count variable if it didn't exist yet */
    dimension_type[0] = harp_dimension_time;
    if (add_count_variable(product, bintype, NULL, 1, dimension_type, &num_bins, count) != 0)
    {
        goto error;
    }

    /* post-process all variables */
    for (k = 0; k < product->num_variables; k++)
    {
        harp_variable *variable;
        long num_sub_elements;
        int count_applied = 0;

        if (bintype[k] == binning_skip || bintype[k] == binning_remove)
        {
            continue;
        }

        variable = product->variable[k];
        num_sub_elements = variable->num_elements / num_elements;

        if (bintype[k] == binning_angle)
        {
            harp_variable *count_variable;

            /* convert angle variables back from complex values to angles */
            for (i = 0; i < variable->num_elements; i += 2)
            {
                variable->data.double_data[i] = atan2(variable->data.double_data[i + 1], variable->data.double_data[i]);
            }
            if (harp_variable_remove_dimension(variable, variable->num_dimensions - 1, 0) != 0)
            {
                goto error;
            }
            /* convert all angles back to the original unit */
            if (harp_convert_unit("rad", variable->unit, variable->num_elements, variable->data.double_data) != 0)
            {
                goto error;
            }

            /* set values to NaN if count==0 */
            if (get_count_variable_for_variable(product, variable, bintype, &count_variable) != 0)
            {
                goto error;
            }
            if (count_variable != NULL)
            {
                double nan_value = harp_nan();
                int count_index;

                for (i = 0; i < variable->num_elements; i++)
                {
                    if (count_variable->data.int32_data[i] == 0)
                    {
                        variable->data.double_data[i] = nan_value;
                    }
                }
                count_applied = 1;

                /* remove the count variable for angles, since it is meaningless for further propagation of averages */
                if (harp_product_get_variable_index_by_name(product, count_variable->name, &count_index) != 0)
                {
                    goto error;
                }
                bintype[count_index] = binning_remove;
            }
        }

        /* divide variables by the sample count and/or set values to NaN if count==0 */
        if (bintype[k] == binning_average)
        {
            double nan_value = harp_nan();
            int result;

            result = get_count_for_variable(product, variable, bintype, filtered_count);
            if (result < 0)
            {
                goto error;
            }
            if (result == 1)
            {
                /* divide by the count (or set value to NaN if count==0) */
                for (i = 0; i < variable->num_elements; i++)
                {
                    if (filtered_count[i] == 0)
                    {
                        variable->data.double_data[i] = nan_value;
                    }
                    else if (filtered_count[i] > 1)
                    {
                        variable->data.double_data[i] /= filtered_count[i];
                    }
                }
                count_applied = 1;
            }
        }

        /* set all empty bins to NaN (for double) or 0 (for int32) */
        if (!count_applied)
        {
            for (i = 0; i < num_bins; i++)
            {
                if (count[i] == 0)
                {
                    double nan_value = harp_nan();

                    if (variable->data_type == harp_type_int32)
                    {
                        for (j = 0; j < num_sub_elements; j++)
                        {
                            variable->data.int32_data[i * num_sub_elements + j] = 0;
                        }
                    }
                    else
                    {
                        for (j = 0; j < num_sub_elements; j++)
                        {
                            variable->data.double_data[i * num_sub_elements + j] = nan_value;
                        }
                    }
                }
            }
        }
    }

    /* remove all variables that need to be removed (in reverse order!) */
    for (k = product->num_variables - 1; k >= 0; k--)
    {
        if (bintype[k] == binning_remove)
        {
            if (harp_product_remove_variable(product, product->variable[k]) != 0)
            {
{'fraction_non_alphanumeric': 0.09785469335672702, 'fraction_numerical': 0.007943217162453145, 'mean_word_length': 1.6856286142642964, 'pattern_counts': {'":': 0, '<': 87, '<?xml version=': 0, '>': 345, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 41, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4516373', 'n_tokens_mistral': 18546, 'n_tokens_neox': 17259, 'n_words': 4411}