Function index (Python)

`choices`

`choice_file_combined`

Aggregate trial-level comparison judgments into a consolidated choice file.

This script is the downstream aggregation step following rs_detailed_choice.py, which produces trial-level triadic (or tetradic) distance comparisons from ranking experiment click-order data. Unlike the preceding utilities, which emit one row per comparison per trial, this module collapses across trials to produce a summary of unique comparisons.

The script does the following:

Takes as input a detailed choices file in which the same comparison (defined by two pairs of stimulus IDs) may appear multiple times.
Groups comparisons that are identical under canonicalization (standardized stimulus ordering).
Tallies how many times each unique comparison appears, stored in a new column N(comparisons).
Sums judgments across occurrences of each comparison, yielding the total number of times the first stimulus pair was judged more dissimilar than the second.
Drops all trial-level information; the trial field is not retained.

Notes: The resulting file is intended for downstream modeling and inference methods that operate on aggregated choice counts rather than individual trials, such as distance-based embedding and psychophysical scaling approaches (e.g., Waraich and Victor, 2022; 2024).

This assumes comparisons have already been standardized and, if
applicable, stimulus labels replaced with numeric IDs.

`build_combine_choice_mat(input_mat_path, output_dir, exp_name, subject)`

Combine repeated triadic or tetradic comparison judgments.

This function inspects the response column names in the input .mat file to determine whether the file contains triadic comparisons (ref, s1, s2) or tetradic comparisons (s1, s2, s3, s4). It then dispatches to the appropriate aggregation function to combine repeated judgments and save the resulting .mat file.

Parameters:

Name	Type	Description	Default
`input_mat_path`	`str`	Path to the input `.mat` file containing the fields `responses`, `responses_colnames`, and `metadata`.	required
`output_dir`	`str`	Directory where the combined `.mat` file will be saved.	required
`exp_name`	`str`	Experiment name used to build the output filename.	required
`subject`	`str`	Subject identifier used to build the output filename.	required

Returns:

Name	Type	Description
`None`		Calls the appropriate aggregation function and writes the combined `.mat` file to `output_dir`.

`build_combine_choice_mat_triadic_format(input_mat_path, output_dir, exp_name, subject)`

Combine judgments across repeated triadic or tetradic comparisons.

This function loads a .mat file containing trial-level comparison judgments, aggregates repeated comparisons by summing the judgment counts, and records how many times each unique comparison was repeated. The result is saved as a new combined .mat file in the specified output directory.

Parameters:

Name	Type	Description	Default
`input_mat_path`	`str`	Path to the input .mat file containing the fields responses, responses_colnames, and metadata.	required
`output_dir`	`str`	Directory where the combined .mat file will be saved.	required
`exp_name`	`str`	Experiment name used to build the output filename.	required
`subject`	`str`	Subject identifier used to build the output filename.	required

Returns:

Name	Type	Description
`None`		Saves a .mat file named {exp_name}combined_choices.mat in output_dir. Prints the full output path after saving.

Source code in src/rs_py/choices/choice_file_combined.py

def build_combine_choice_mat_triadic_format(input_mat_path, output_dir, exp_name, subject):
    """
        Combine judgments across repeated triadic or tetradic comparisons.

        This function loads a .mat file containing trial-level comparison judgments,
        aggregates repeated comparisons by summing the judgment counts, and records how
        many times each unique comparison was repeated. The result is saved as a new
        combined .mat file in the specified output directory.

        Args:
            input_mat_path (str): Path to the input .mat file containing the fields
                responses, responses_colnames, and metadata.
            output_dir (str): Directory where the combined .mat file will be saved.
            exp_name (str): Experiment name used to build the output filename.
            subject (str): Subject identifier used to build the output filename.

        Returns:
            None:
                Saves a .mat file named {exp_name}_combined_choices_{subject}.mat in
                output_dir. Prints the full output path after saving.
    """
    data = loadmat(input_mat_path, squeeze_me=True)

    responses = data['responses']
    colnames = [name.strip() for name in data['responses_colnames']]
    metadata = data['metadata']
    stim_list = data['stim_list']

    COL_REF = colnames.index("ref")
    COL_S1 = colnames.index("s1")
    COL_S2 = colnames.index("s2")
    COL_JUDGMENT = colnames.index("N(D(ref, s1) > D(ref, s2))")

    aggregated_responses = {}
    comparison_repeats = {}

    for row in responses:
        key = (
            int(row[COL_REF]),
            int(row[COL_S1]),
            int(row[COL_S2])
        )
        if key not in aggregated_responses:
            aggregated_responses[key] = row[COL_JUDGMENT]
            comparison_repeats[key] = 1
        else:
            aggregated_responses[key] += row[COL_JUDGMENT]
            comparison_repeats[key] += 1

    comb_response_colnames = ['ref', 's1', 's2', 'N(D(ref, s1) > D(ref, s2))', 'N_Repeats(D(ref, s1) > D(ref, s2))']
    col_ref = 0
    col_s1 = 1
    col_s2 = 2
    col_judgment = 3
    col_count = 4

    comb_responses = np.zeros((len(comparison_repeats), len(comb_response_colnames)), dtype=int)
    i = 0
    for key in aggregated_responses:
        comb_responses[i, col_ref] = key[0]
        comb_responses[i, col_s1] = key[1]
        comb_responses[i, col_s2] = key[2]
        comb_responses[i, col_judgment] = aggregated_responses[key]
        comb_responses[i, col_count] = comparison_repeats[key]
        i += 1

    combined_choices = {
        'metadata': metadata,
        'stim_list': stim_list,
        'responses_colnames': comb_response_colnames,
        'responses': comb_responses
    }
    output_path = os.path.join(output_dir, f"{exp_name}_choices_{subject}.mat")
    savemat(output_path, combined_choices)
    print(f"Saved results to {output_path}")

`build_combined_mat_tetradic_format(input_mat_path, output_dir, exp_name, subject)`

Combine judgments across repeated tetradic comparisons.

This function loads a .mat file containing trial-level comparison judgments, aggregates repeated comparisons by summing the judgment counts, and records how many times each unique comparison was repeated. The result is saved as a new combined .mat file in the specified output directory.

Parameters:

Name	Type	Description	Default
`input_mat_path`	`str`	Path to the input `.mat` file containing the fields `responses`, `responses_colnames`, and `metadata`.	required
`output_dir`	`str`	Directory where the combined `.mat` file will be saved.	required
`exp_name`	`str`	Experiment name used to build the output filename.	required
`subject`	`str`	Subject identifier used to build the output filename.	required

Returns:

Name	Type	Description
`None`		Saves a `.mat` file named `{exp_name}_combined_choices_{subject}.mat` in `output_dir` and prints the full output path after saving.

Source code in src/rs_py/choices/choice_file_combined.py

def build_combined_mat_tetradic_format(input_mat_path, output_dir, exp_name, subject):
    """
        Combine judgments across repeated tetradic comparisons.

        This function loads a `.mat` file containing trial-level comparison judgments,
        aggregates repeated comparisons by summing the judgment counts, and records how
        many times each unique comparison was repeated. The result is saved as a new
        combined `.mat` file in the specified output directory.

        Args:
            input_mat_path (str): Path to the input `.mat` file containing the fields
                `responses`, `responses_colnames`, and `metadata`.
            output_dir (str): Directory where the combined `.mat` file will be saved.
            exp_name (str): Experiment name used to build the output filename.
            subject (str): Subject identifier used to build the output filename.

        Returns:
            None:
                Saves a `.mat` file named `{exp_name}_combined_choices_{subject}.mat` in
                `output_dir` and prints the full output path after saving.
    """

    data = loadmat(input_mat_path, squeeze_me=True)

    responses = data['responses']
    colnames = [name.strip() for name in data['responses_colnames']]
    metadata = data['metadata']
    stim_list = data['stim_list']

    COL_S1 = colnames.index("s1")
    COL_S2 = colnames.index("s2")
    COL_S3 = colnames.index("s3")
    COL_S4 = colnames.index("s4")
    COL_JUDGMENT = colnames.index("N(D(s1, s2) > D(s3, s4))")

    aggregated_responses = {}
    comparison_repeats = {}

    for row in responses:
        key = (
            int(row[COL_S1]),
            int(row[COL_S2]),
            int(row[COL_S3]),
            int(row[COL_S4])
        )
        if key not in aggregated_responses:
            aggregated_responses[key] = row[COL_JUDGMENT]
            comparison_repeats[key] = 1
        else:
            aggregated_responses[key] += row[COL_JUDGMENT]
            comparison_repeats[key] += 1

    comb_response_colnames = ['s1', 's2', 's3', 's4', 'N(D(s1, s2) > D(s3, s4))',
                              'N_Repeats(D(s1, s2) > D(s3, s4))']
    COL_S1 = 0
    COL_S2 = 1
    COL_S3 = 2
    COL_S4 = 3
    COL_JUDGMENT = 4
    COL_COUNT = 5

    comb_responses = np.zeros((len(comparison_repeats), len(comb_response_colnames)), dtype=int)
    i = 0
    for key in aggregated_responses:
        comb_responses[i, COL_S1] = key[0]
        comb_responses[i, COL_S2] = key[1]
        comb_responses[i, COL_S3] = key[2]
        comb_responses[i, COL_S4] = key[3]
        comb_responses[i, COL_JUDGMENT] = aggregated_responses[key]
        comb_responses[i, COL_COUNT] = comparison_repeats[key]
        i += 1

    combined_choices = {
        'metadata': metadata,
        'stim_list': stim_list,
        'responses_colnames': comb_response_colnames,
        'responses': comb_responses
    }
    output_path = os.path.join(output_dir, f"{exp_name}_choices_{subject}.mat")
    savemat(output_path, combined_choices)
    print(f"Saved results to {output_path}")
    return

`choice_file_detailed`

Utilities for converting ranking experiment click-order data into triadic or tetradic distance comparisons.

This module reads response CSV files from ranking experiments, parses per-trial click orderings relative to a reference stimulus, and generates all corresponding pairwise distance comparisons. It provides functionality to represent pairwise comparisons in a standardized way, replace stimulus labels with numeric IDs, and prepare data for downstream analysis (e.g., MATLAB .mat files), following the framework of Waraich and Victor (2022).

Each row in the output (.mat) file contains information about a judgment between two pairs of stimuli, and a trial number. It contains a matlab struct, with one field for the trial-by-trial judgments, and one field for metadata. The metadata in turn contains information regarding the stimuli, paradigm, number of trials, or any optional information a user may find useful to add. In our case, (Waraich and Victor, 2024), metadata contains a list of stimuli, number of sessions, the subject identifier, and the task.

`generate_comparisons(reference, clicks, trial_num)`

Generate all pairwise triadic distance comparisons for a single trial.

Given a reference stimulus and an ordered list of non-reference stimuli, generate all pairwise triadic comparisons implied by the click order.

Parameters:

Name	Type	Description	Default
`reference`	`str`	Label of the reference stimulus.	required
`clicks`	`list[str]`	Ordered list of non-reference stimuli in the order they were clicked.	required
`trial_num`	`int`	Trial index assigned to all generated comparisons.	required

Returns:

Type	Description
	comparisons (list[dict]) - List of comparison dictionaries, one for each unordered stimulus pair.
	Each dictionary contains: - trial - s1 - s2 - s3 - s4 - operator - judgment

Notes

The number of generated comparisons is $n# choose 2, where $n$ is the number of comparison stimuli. Each unordered stimulus pair appears exactly once per trial. Canonicalization of comparison keys and stimulus-ID remapping are performed by downstream processing steps.

`parse_click_sequence(row)`

Parse the click order for a trial into a sequence of stimulus labels.

Given a trial row, read the stored click-order indices from row['clicks'] and return the corresponding stimuli in the order they were clicked.

Parameters:

Name	Type	Description	Default
`row`	`dict`	A single trial row containing a `clicks` field and a list of stimuli. The `clicks` field must encode the clicked stimulus indices in the order they were selected.	required

Returns:

Type	Description
	sequence (list[str]) - Ordered list of stimulus labels corresponding to the clicked stimuli for that trial.

Notes

Assumes data are being read in from csv files created by Waraich and Victor (2022) paradigm.

Source code in src/rs_py/choices/choice_file_detailed.py

def parse_click_sequence(row):
    """
    Parse the click order for a trial into a sequence of stimulus labels.

    Given a trial row, read the stored click-order indices from `row['clicks']`
    and return the corresponding stimuli in the order they were clicked.

    Args:
        row (dict):
            A single trial row containing a `clicks` field and a list of stimuli.
            The `clicks` field must encode the clicked stimulus indices in the
            order they were selected.

    Returns:
        - sequence (list[str]) - Ordered list of stimulus labels corresponding
            to the clicked stimuli for that trial.

    Notes:
        - Assumes data are being read in from csv files created by Waraich and Victor (2022) paradigm.
    """
    sequence = ast.literal_eval(row['clicks'])
    return [row[stim_num] for stim_num in sequence]

`process_subject_data(input_directory)`

Read response CSV files for a single subject and convert each trial’s click order into all implied pairwise triadic distance comparisons relative to the reference stimulus.

For every row in every matching response CSV file, the function parses the observed click order, generates all (n choose 2) comparisons among the comparison stimuli, and assigns trial numbers sequentially across all files and rows.

Parameters:

Name	Type	Description	Default
`input_directory`	`str`	Path to a directory containing one or more response CSV files for a single subject. Subdirectories are searched recursively for files matching the response-file pattern defined by `get_response_files`.	required

Returns:

Type	Description
	all_comparisons (list[dict]) - Flat list of comparison dictionaries, one for each generated triadic comparison. Each dictionary contains: trial s1 s2 operator s3 s4 judgment
	stimuli (set) - Labels of all stimuli encountered across all processed trials, including both reference and non-reference stimuli.

Notes

Trial numbering starts at 1 and increases sequentially across all files and rows; no session boundaries are inferred.
To ensure compatibility with experiments that utilize multiple exemplars per stimulus, stimuli are generated from the stimulus names in the CSV files after ignoring the text after an underscore. If a stimulus name contains no underscore, then this is not a problem. This affects experiments where an experimenter uses multiple images/ exemplars per stimulus which are coded as 'stim23_002', 'stim23_003' etc. In this case both instances will be tied to one stimulus called 'stim23'.
CSV files are read using UTF-8 with BOM (utf-8-sig) encoding to match experimental data exports.
This function does not canonicalize comparison keys, remap stimulus IDs, or aggregate across trials; those steps are handled downstream.
The interpretation of click order and judgment semantics is handled by generate_comparisons.

`replace_stimuli_with_ids(comparisons, stimuli_set)`

Replace stimulus labels in comparison dictionaries with integer stimulus IDs.

Constructs a deterministic mapping from stimulus labels to integer IDs based on the sorted order of stimuli_set, then replaces the values of the stimulus fields (s1, s2, s3, s4) in each comparison dictionary with their corresponding integer IDs.

Parameters:

Name	Type	Description	Default
`comparisons`	`list[dict]`	List of comparison dictionaries. Each dictionary must contain the keys `s1`, `s2`, `s3`, and `s4`, whose values are stimulus labels present in `stimuli_set`. The list is modified in place.	required
`stimuli_set`	`set`	Set of all stimulus labels to be mapped to integer IDs. IDs are assigned starting from 1, in alphabetical order of the labels.	required

Returns:

Type	Description
	comparisons (list[dict]) - The input list of comparison dictionaries with stimulus labels replaced by integer stimulus IDs. The mapping of IDs to stimuli is also returned or stored separately, depending on downstream use.

Notes

The function mutates the input comparisons list in place.
Judgment values and non-stimulus fields are not modified.
If a comparison contains a stimulus label not present in stimuli_set, a KeyError will be raised.
The function does not perform canonicalization, reordering of stimulus keys, or judgment flipping.

`standardize_comparison_keys(comparisons, comparison_type='triadic')`

Standardize comparison field order and canonicalize comparison keys.

Takes in a list of comparisons, edits their fields (s1, s2, s3, s4, judgment), and returns the same list with comparison keys ordered consistently. Equivalent comparisons are therefore represented identically, which supports downstream tallying and aggregation.

Parameters:

Name	Type	Description	Default
`comparisons`	`list[dict]`	A list of comparisons and judgments with the keys `s1`, `s2`, `s3`, `s4`, and `judgment`.	required
`comparison_type`	`str`	Comparison type to standardize. `'triadic'` by default, and the only type implemented as of 02/17/2026.	`'triadic'`

Returns:

Type	Description
	comparisons (list[dict]) - The same list of comparisons with field names and comparison order standardized. Each dictionary contains: s1 s2 s3 s4 judgment

Notes

If (s2, s4) appears in reverse order, the elements are swapped and the judgment is flipped.
For triadic comparisons, the reference appears in both pairs of stimuli and is always assigned as s1 = s3 = ref.
The ordering of the two stimulus pairs is determined by alphabetical order of the non-reference elements.
Example: for pairs (k, l) and (c, k), the standardized form is s1 = k, s2 = c, s3 = k, s4 = l.
Tetradic comparisons are not yet implemented.
Within each pair, stimulus order is also determined alphabetically.
Example: for pairs (k, l) and (h, w), the standardized form is s1 = h, s2 = w, s3 = k, s4 = l.

`choice_likelihoods`

Utilities for computing log-likelihoods for geometry-based comparison data.

This file provides methods for evaluating the log-likelihood of aggregated comparison counts under a distance-geometry response model. The observed data are of the form $N((i, j) > (k, l))$, meaning the number of times one pair of stimuli was judged more different than another pair.

The model computes choice probabilities from interstimulus distances using an error-function approximation, then evaluates the corresponding log-likelihood. It also includes reference calculations for random-choice and empirical best-case benchmarks.

Note: The likelihood calculations use log2, so the returned values are in bits. When sigma_point = 0, the Gaussian-noise assumption matches the erf-based probability exactly. When sigma_point > 0, the erf probability is only an approximation because the full noise model is not strictly Gaussian. (By sigma_point, we mean the standard deviation of additive noise on stimulus coordinates before distances are computed.)

`best_model_ll(judgments, judgment_repeats)`

Compute the log-likelihood under an empirical best-case benchmark.

This benchmark assigns each comparison the probability estimated from its observed frequency in the data.

Parameters:

Name	Type	Description	Default
`judgments`	`dict`	Mapping from comparison keys to observed judgment counts.	required
`judgment_repeats`	`dict`	Mapping from comparison keys to total repeats.	required

Returns:

Type	Description
	tuple[float, bool]: The log-likelihood and infeasibility flag returned by `calculate_ll`.

Source code in src/rs_py/choices/choice_likelihoods.py

def best_model_ll(judgments, judgment_repeats):
    """
        Compute the log-likelihood under an empirical best-case benchmark.

        This benchmark assigns each comparison the probability estimated from its
        observed frequency in the data.

        Args:
            judgments (dict): Mapping from comparison keys to observed judgment counts.
            judgment_repeats (dict): Mapping from comparison keys to total repeats.

        Returns:
            tuple[float, bool]: The log-likelihood and infeasibility flag returned by `calculate_ll`.
    """
    counts = []
    probs = []
    repeats = []
    for (k, v) in judgments.items():
        counts.append(v)
        probs.append(v / judgment_repeats[k])
        repeats.append(judgment_repeats[k])
    return calculate_ll(np.array(counts), np.array(probs), np.array(repeats))

`calculate_ll(counts, probs, num_repeats)`

Calculate the log-likelihood of observed counts given model probabilities.

This function augments each comparison with its reverse outcome, checks whether the model assigns zero probability to any observed event, and returns the log-likelihood using a base-2 logarithm.

Parameters:

Name	Type	Description	Default
`counts`	`ndarray`	Observed counts for each comparison direction.	required
`probs`	`ndarray`	Model probabilities for the corresponding comparisons.	required
`num_repeats`	`ndarray`	Total number of repeats for each comparison.	required

Returns:

Type	Description
	tuple[float, bool]: The log-likelihood and a boolean flag indicating whether the model is infeasible.
	The flag is `True` when an observed event has zero model probability.

Note

The function uses a variable called EPSILON set to 1e-30. This very small positive constant is added only when a model assigns zero probability to an observed comparison. It prevents log2(0) from causing an error during likelihood computation, while keeping the adjustment negligible.

Source code in src/rs_py/choices/choice_likelihoods.py

def calculate_ll(counts, probs, num_repeats):
    """
        Calculate the log-likelihood of observed counts given model probabilities.

        This function augments each comparison with its reverse outcome, checks whether
        the model assigns zero probability to any observed event, and returns the
        log-likelihood using a base-2 logarithm.

        Args:
            counts (numpy.ndarray): Observed counts for each comparison direction.
            probs (numpy.ndarray): Model probabilities for the corresponding
                comparisons.
            num_repeats (numpy.ndarray): Total number of repeats for each comparison.

        Returns:
            tuple[float, bool]: The log-likelihood and a boolean flag indicating whether the model is infeasible.
            The flag is `True` when an observed event has zero model probability.

        Note:
            The function uses a variable called `EPSILON` set to `1e-30.` This very small positive constant is added only when a model assigns
            zero probability to an observed comparison. It prevents `log2(0)` from causing an error
            during likelihood computation, while keeping the adjustment negligible.
    """

    reverse_counts = num_repeats - counts
    reverse_probs = 1 - probs

    probs = concatenate((probs, reverse_probs))
    counts = concatenate((counts, reverse_counts))

    model_bad = False
    # check if geometry is bad, i.e. prob = 0 but count > 0
    prob_zero = probs == 0
    if (counts[prob_zero] > 0).any():
        model_bad = True

    # make sure there are no zero probabilities (avoid log(0) error)
    probs[prob_zero] += EPSILON
    log_likelihood = counts.dot(log2(probs))
    return log_likelihood, model_bad

`dist_model_ll_vectorized(pair_a, pair_b, judgment_counts, judgment_repeats, params, stimuli)`

Compute the log-likelihood under a specific perceptual space model with given distances between points.

This function converts stimulus coordinates into a full pairwise distance matrix, computes comparison probabilities with find_probabilities, and then evaluates the log-likelihood against the observed judgment counts.

Parameters:

Name	Type	Description	Default
`pair_a`	`ndarray`	Array of index pairs defining the first distance in each comparison.	required
`pair_b`	`ndarray`	Array of index pairs defining the second distance in each comparison.	required
`judgment_counts`	`ndarray`	Observed counts for the first pair being judged more dissimilar than the second.	required
`judgment_repeats`	`ndarray`	Total number of repeats for each comparison.	required
`params`	`dict`	Model parameters. Must contain `noise_st_dev`.	required
`stimuli`	`ndarray`	Stimulus coordinates used to compute pairwise distances.	required

Returns:

Type	Description
	tuple[float, bool]: The log-likelihood and the infeasibility flag returned by `calculate_ll`.

`find_probabilities(distances, pair_a, pair_b, noise_st_dev)`

Compute choice probabilities from pairwise distance differences.

For each comparison, this function subtracts the second pair distance from the first pair distance and converts that difference into a probability that the first pair is judged more different.

Parameters:

Name	Type	Description	Default
`distances`	`ndarray`	Square matrix of pairwise stimulus distances.	required
`pair_a`	`ndarray`	Array of index pairs defining the first distance in each comparison.	required
`pair_b`	`ndarray`	Array of index pairs defining the second distance in each comparison.	required
`noise_st_dev`	`float`	Standard deviation of the combined Gaussian noise term.	required

Returns:

Type	Description
	numpy.ndarray: Probability that the first pair is judged more different than the second.

Note

When noise_st_dev = 0, the output is deterministic, with probabilities 0, 0.5, or 1 depending on the sign of the distance difference. When noise_st_dev > 0, the probability is computed using an approximation based on the erf function.

Source code in src/rs_py/choices/choice_likelihoods.py

def find_probabilities(distances, pair_a, pair_b, noise_st_dev):
    """
        Compute choice probabilities from pairwise distance differences.

        For each comparison, this function subtracts the second pair distance from the
        first pair distance and converts that difference into a probability that the
        first pair is judged more different.

        Args:
            distances (numpy.ndarray): Square matrix of pairwise stimulus distances.
            pair_a (numpy.ndarray): Array of index pairs defining the first distance in
                each comparison.
            pair_b (numpy.ndarray): Array of index pairs defining the second distance in each comparison.
            noise_st_dev (float): Standard deviation of the combined Gaussian noise term.

        Returns:
            numpy.ndarray: Probability that the first pair is judged more different than the second.

        Note:
            When `noise_st_dev = 0`, the output is deterministic, with probabilities `0`, `0.5`, or `1`
            depending on the sign of the distance difference. When `noise_st_dev > 0`, the probability
             is computed using an approximation based on the erf function.
    """
    difference = distances[pair_a[:, 0], pair_a[:, 1]] - distances[pair_b[:, 0], pair_b[:, 1]]
    if noise_st_dev == 0:
        # if (sigmas['compare'] + sigmas['dist'] == 0) or no_noise is True:
        probabilities = (difference < 0) * 0 + (difference > 0) * 1 + (difference == 0) * 0.5
    else:
        # total_st_dev = sqrt((sigmas['dist'] ** 2) + sigmas['compare'] ** 2)
        probabilities = 0.5 * (1 + erf(difference / float(2 * noise_st_dev)))
    return probabilities

`random_choice_ll(judgments, judgment_repeats)`

Compute the log-likelihood under a random-choice benchmark.

This benchmark assumes every comparison direction is equally likely, so each observed judgment has probability 0.5.

Parameters:

Name	Type	Description	Default
`judgments`	`dict`	Mapping from comparison keys to observed judgment counts.	required
`judgment_repeats`	`dict`	Mapping from comparison keys to total repeats.	required

Returns:

Type	Description
	tuple[float, bool]: The log-likelihood and infeasibility flag returned by `calculate_ll`.

Source code in src/rs_py/choices/choice_likelihoods.py

def random_choice_ll(judgments, judgment_repeats):
    """
    Compute the log-likelihood under a random-choice benchmark.

    This benchmark assumes every comparison direction is equally likely, so each
    observed judgment has probability 0.5.

    Args:
        judgments (dict): Mapping from comparison keys to observed judgment counts.
        judgment_repeats (dict): Mapping from comparison keys to total repeats.

    Returns:
        tuple[float, bool]: The log-likelihood and infeasibility flag returned by `calculate_ll`.
 """
    counts = []
    probs = []
    repeats = []
    for (k, v) in judgments.items():
        counts.append(v)
        probs.append(0.5)
        repeats.append(judgment_repeats[k])
    return calculate_ll(np.array(counts), np.array(probs), np.array(repeats))

`demos`

`demo_combined_choices`

To demo the script that aggregates detailed trial-level choice judgments into a combined (unique-comparison) choice file.

This demos

combine_choice_mat(input_mat_path, output_dir, exp_name, subject)

Enter 0 to use default values.

`demo_detailed_choices`

Demo: create a detailed choices .mat file from raw ranking judgments.

This script

Parses raw subject CSV files from a ranking experiment.
Generates triadic or tetradic comparisons.
Replaces stimulus labels with numeric IDs.
Canonicalizes comparison keys.
Saves a detailed trial-level .mat file.

Enter 0 to use default values.

`demo_fit_euclidean`

Demo to fit models of Euclidean spaces of different dimensions, from tallied similarity judgments. Enter '0' to use default values.

`demo_inputs()`

Populate demo defaults. Adjust the default filepath/outdir to wherever your sample materials live.

Source code in src/rs_py/demos/demo_fit_euclidean.py

def demo_inputs():
    """
    Populate demo defaults.
    Adjust the default filepath/outdir to wherever your sample materials live.
    """
    base_dir = Path(__file__).resolve().parent.parent
    model_defaults = CONFIG["inputs"]["model_fit"]

    defaults = {
        "filepath": (base_dir / "samples/choice_files/image_choices_S4.mat").resolve(),
        "exp_name": "image",
        "subject": "S4",
        "output_dir": (base_dir / "samples/outputs").resolve(),
        "sigma": model_defaults['sigma'],
        "model_dimensions": model_defaults['model_dimensions'],
        "learning_rate": model_defaults['learning_rate'],
        "tolerance": model_defaults['tolerance'],
        "max_trials": model_defaults['max_trials'],
        "max_iterations": model_defaults['max_iterations'],
        "minimization": model_defaults['minimization']
    }
    # Creates the folder if missing; does nothing if it already exists
    defaults["output_dir"].mkdir(parents=True, exist_ok=True)

    return defaults

`geometry`

`euclidean`

Sample points from Euclidean spaces of configurable dimension.

This module provides utilities for drawing stimulus locations in an n-dimensional Euclidean space. Points can be sampled from the surface of a sphere, from a Gaussian distribution, or uniformly from inside a sphere.

`EuclideanSpace`

Represent an n-dimensional Euclidean space for sampling points.

Source code in src/rs_py/geometry/euclidean.py

class EuclideanSpace:
    """Represent an n-dimensional Euclidean space for sampling points."""

    def __init__(self, num_dim):
        """
            Set the dimensionality of the space.
            Args:
                num_dim (int): Number of dimensions in the Euclidean space. """
        # NOTE: "Sampling strategy fails for space with less than 2 dimensions."
        self.dimensions = num_dim

    def sample_space(self, magnitude, method="spherical_shell"):
        """
            Sample one point from the space.
            Args:
                magnitude (float): Target magnitude or scale of the sample.
                method (str): Sampling method. Options are `spherical_shell`, `gaussian`, and `uniform`.

            Returns:
                numpy.ndarray or float: A sampled point. For `spherical_shell`, returns a scalar when `dimensions == 1`.
        """
        if method == "spherical_shell":
            if self.dimensions == 1:
                return np.random.normal(0, 1)
            else:
                sample = np.random.normal(0, 1, self.dimensions)
                # get length of vector
                length = np.sqrt(sample.dot(sample))
                # normalize and scale vector so it has magnitude passed in as arg
                scaled_sample = np.array([(float(x) / length) * magnitude for x in sample])
                return scaled_sample
        # sample from inside a Gaussian, not limited to the surface
        elif method == "gaussian":
            sample = np.random.normal(0, magnitude, self.dimensions)
            return sample
        # sample each dimension from a uniform distribution on [0, magnitude]
        # but only take points that lie in a sphere of radius magnitude
        elif method == "uniform":
            sample = np.random.uniform(-magnitude, magnitude, self.dimensions)
            length = np.sqrt(sample.dot(sample))
            while length > magnitude:
                sample = np.random.uniform(-magnitude, magnitude, self.dimensions)
                length = np.sqrt(sample.dot(sample))
            return sample

    def get_samples(self, num_stimuli, magnitude=1, method="spherical_shell"):
        """
            Sample multiple points from the space.
            Args:
                num_stimuli (int): Number of points to sample.
                magnitude (float): Target magnitude or scale of each sample.
                method (str): Sampling method. Options are `spherical_shell`, `gaussian`, and `uniform`.

            Returns:
                numpy.ndarray: Array of sampled points, one per stimulus.
        """
        return np.array([self.sample_space(magnitude, method) for _ in range(num_stimuli)])

`init(num_dim)`

Set the dimensionality of the space. Args: num_dim (int): Number of dimensions in the Euclidean space.

Source code in src/rs_py/geometry/euclidean.py

def __init__(self, num_dim):
    """
        Set the dimensionality of the space.
        Args:
            num_dim (int): Number of dimensions in the Euclidean space. """
    # NOTE: "Sampling strategy fails for space with less than 2 dimensions."
    self.dimensions = num_dim

`get_samples(num_stimuli, magnitude=1, method='spherical_shell')`

Sample multiple points from the space. Args: num_stimuli (int): Number of points to sample. magnitude (float): Target magnitude or scale of each sample. method (str): Sampling method. Options are spherical_shell, gaussian, and uniform.

Returns:

Type	Description
	numpy.ndarray: Array of sampled points, one per stimulus.

Source code in src/rs_py/geometry/euclidean.py

def get_samples(self, num_stimuli, magnitude=1, method="spherical_shell"):
    """
        Sample multiple points from the space.
        Args:
            num_stimuli (int): Number of points to sample.
            magnitude (float): Target magnitude or scale of each sample.
            method (str): Sampling method. Options are `spherical_shell`, `gaussian`, and `uniform`.

        Returns:
            numpy.ndarray: Array of sampled points, one per stimulus.
    """
    return np.array([self.sample_space(magnitude, method) for _ in range(num_stimuli)])

`sample_space(magnitude, method='spherical_shell')`

Sample one point from the space. Args: magnitude (float): Target magnitude or scale of the sample. method (str): Sampling method. Options are spherical_shell, gaussian, and uniform.

Returns:

Type	Description
	numpy.ndarray or float: A sampled point. For `spherical_shell`, returns a scalar when `dimensions == 1`.

Source code in src/rs_py/geometry/euclidean.py

def sample_space(self, magnitude, method="spherical_shell"):
    """
        Sample one point from the space.
        Args:
            magnitude (float): Target magnitude or scale of the sample.
            method (str): Sampling method. Options are `spherical_shell`, `gaussian`, and `uniform`.

        Returns:
            numpy.ndarray or float: A sampled point. For `spherical_shell`, returns a scalar when `dimensions == 1`.
    """
    if method == "spherical_shell":
        if self.dimensions == 1:
            return np.random.normal(0, 1)
        else:
            sample = np.random.normal(0, 1, self.dimensions)
            # get length of vector
            length = np.sqrt(sample.dot(sample))
            # normalize and scale vector so it has magnitude passed in as arg
            scaled_sample = np.array([(float(x) / length) * magnitude for x in sample])
            return scaled_sample
    # sample from inside a Gaussian, not limited to the surface
    elif method == "gaussian":
        sample = np.random.normal(0, magnitude, self.dimensions)
        return sample
    # sample each dimension from a uniform distribution on [0, magnitude]
    # but only take points that lie in a sphere of radius magnitude
    elif method == "uniform":
        sample = np.random.uniform(-magnitude, magnitude, self.dimensions)
        length = np.sqrt(sample.dot(sample))
        while length > magnitude:
            sample = np.random.uniform(-magnitude, magnitude, self.dimensions)
            length = np.sqrt(sample.dot(sample))
        return sample

`hyperbolic`

Add curvature to a Euclidean fit and test whether the fit improves.

This module takes a Euclidean embedding that already explains psychophysical judgments, maps the points onto the top sheet of a two-sheet hyperboloid using Loid geometry, and then computes pairwise hyperbolic distances. The resulting distances can be used to evaluate whether adding negative curvature improves log-likelihood relative to the Euclidean model.

The mapping is controlled by a curvature parameter, which determines how close the projected points lie to the hyperboloid center. When the curvature parameter approaches 0, the geometry approaches Euclidean space, so improved fit at positive curvature provides evidence for curvature in the space, after accounting for the extra parameter.

`hyperbolic_distances(X, curvature)`

Compute pairwise hyperbolic distances between points on the hyperboloid.

Parameters:

Name	Type	Description	Default
`X`	`ndarray`	Array of shape (d + 1, n) containing points already projected onto the hyperboloid.	required
`curvature`	`float`	Positive curvature parameter used to scale the resulting distances.	required

Returns:

Type	Description
	numpy.ndarray: Pairwise hyperbolic distance matrix for the points in X.

Notes

Distances are computed as $arccosh(-[X, X])$.

Source code in src/rs_py/geometry/hyperbolic.py

def hyperbolic_distances(X, curvature):
    """
    Compute pairwise hyperbolic distances between points on the hyperboloid.

    Args:
        X (numpy.ndarray): Array of shape (d + 1, n) containing points already
            projected onto the hyperboloid.
        curvature (float): Positive curvature parameter used to scale the
            resulting distances.

    Returns:
        numpy.ndarray: Pairwise hyperbolic distance matrix for the points in X.

    Notes:
        Distances are computed as $arccosh(-[X, X])$.
    """
    # test entries along the diagonal should equal 1
    # test all entries should be less than or equal to -1 or what notes say
    H = np.eye(X.shape[0])
    H[0, 0] = -1
    inner_product = X.T @ H @ X
    # return interstimulus_distances
    return np.arccosh(-round(inner_product, 6)) / curvature

`loid_map(X, degree_curvature)`

Map Euclidean points onto the top sheet of a two-sheet hyperboloid.

Parameters:

Name	Type	Description	Default
`X`	`ndarray`	Array of shape (d, n), with n points in d Euclidean dimensions.	required
`degree_curvature`	`float`	Curvature parameter controlling how strongly the points are projected toward hyperbolic geometry.	required

Returns:

Type	Description
	numpy.ndarray: Array of shape (d + 1, n) containing the projected points on the hyperboloid.

Source code in src/rs_py/geometry/hyperbolic.py

def loid_map(X, degree_curvature):
    """
        Map Euclidean points onto the top sheet of a two-sheet hyperboloid.

        Args:
            X (numpy.ndarray): Array of shape (d, n), with n points in d Euclidean
                dimensions.
            degree_curvature (float): Curvature parameter controlling how strongly
                the points are projected toward hyperbolic geometry.

        Returns:
            numpy.ndarray: Array of shape (d + 1, n) containing the projected points on the hyperboloid.
    """
    # retain coordinates of X but add a 0-th coordinate which is a function of the d-dimensional coordinate values
    d, n = X.shape
    Y = np.zeros((d + 1, n))
    Y[1:, :] = degree_curvature * X
    dot_prods = np.einsum('ij,ij->j', X, X)
    Y[0, :] = np.sqrt(1 + (degree_curvature ** 2) * dot_prods)
    return Y

`spherical`

Sample points on a sphere and compute pairwise spherical distances.

This module provides helpers for mapping Euclidean points onto a spherical surface and measuring pairwise spherical distances. It includes legacy (deprecated) versions of the mapping and distance functions, plus the current implementations used in the codebase.

`sphere_map(X, radius)`

Map Euclidean points to the surface of a sphere.

This function rescales the input by the sphere radius and applies the current spherical mapping used for visualization and analysis.

Parameters:

Name	Type	Description	Default
`X`	`ndarray`	Array of shape (d, n), with n points in d Euclidean dimensions.	required
`radius`	`float`	Radius of the sphere.	required

Returns:

Type	Description
	numpy.ndarray: Array of shape (d + 1, n) containing points projected onto the sphere.

Notes

mu (curvature) is 1/ radius.

Source code in src/rs_py/geometry/spherical.py

def sphere_map(X, radius):
    """
    Map Euclidean points to the surface of a sphere.

    This function rescales the input by the sphere radius and applies the
    current spherical mapping used for visualization and analysis.

    Args:
        X (numpy.ndarray): Array of shape (d, n), with n points in d Euclidean
            dimensions.
        radius (float): Radius of the sphere.

    Returns:
        numpy.ndarray: Array of shape (d + 1, n) containing points projected onto the sphere.

    Notes:
        mu (curvature) is 1/ radius.
    """
    X = X/radius
    # retain coordinates of X but add a 0-th coordinate which is a function of the d-dimensional coordinate values
    d, n = X.shape
    Y = np.zeros((d + 1, n))
    # squared norms of all vectors = dot products
    dot_prods = np.einsum('ij,ij->j', X, X)  # https://stackoverflow.com/questions/6229519/numpy-column-wise-dot-product
    for k in range(1, d+1):
        for p in range(n):
            Y[k, p] = 2 * X[k-1, p]/(1 + dot_prods[p])
    for p in range(n):
        Y[0, p] = (1 - dot_prods[p])/(1 + dot_prods[p])
    return Y

`spherical_distances(X, radius)`

Compute pairwise spherical distances between points on the sphere.

Parameters:

Name	Type	Description	Default
`X`	`ndarray`	Array of shape (d, n) containing points already projected onto the sphere.	required
`radius`	`float`	Radius of the sphere.	required

Returns:

Type	Description
	numpy.ndarray: Pairwise spherical distance matrix.

Notes

Assumes the points passed in are not off the surface - already projected! Distance is $arccos(-[X, X])$ for all points X

Source code in src/rs_py/geometry/spherical.py

def spherical_distances(X, radius):
    """
    Compute pairwise spherical distances between points on the sphere.

    Args:
        X (numpy.ndarray): Array of shape (d, n) containing points already
            projected onto the sphere.
        radius (float): Radius of the sphere.

    Returns:
        numpy.ndarray: Pairwise spherical distance matrix.

    Notes:
        Assumes the points passed in are not off the surface - already projected!
        Distance is $arccos(-[X, X])$ for all points X
    """
    # standard inner product
    inner_product = X.T @ X  # do not make it negative
    interstimulus_distances = (radius/2) * np.arccos(round(inner_product, 6))
    return interstimulus_distances

`model`

`compare_curved_space_models`

Here, given a data set, generate surrogate datasets that are drawn from the original choice probablities. Then, for each surrogate dataset, essentially calculate LLs for a range of curvatures and (chosen/ adjusted) sigma values.

This should yield a scatterplot of model LLs by scripts values. We start by doing so for the 2D case and we can progress to higher dimensions as needed.

## DEPRECATED - DO NOT USE

The inputs are - a range of scripts values - dimensions for which to apply the analysis. Minimum dimension is 2. - path to json files for individual datasets. - optionally, can and probably should pass in a legend for when multiple curves are drawn - number of iterations (surrogates) - max number of iterations to run the model_fitting pipeline with.

`fit_geometric_models`

`hyperbolic_points_of_best_fit(judgments, number_repeats, args, start_points=None)`

Fit hyperbolic coordinates to the observed judgments.

Parameters:

Name	Type	Description	Default
`judgments`	`dict`	Pairwise comparison judgments.	required
`number_repeats`	`dict`	Number of repeats for each comparison key in `judgments`.	required
`args`	`dict`	Model settings. Must include: - `n_dim`: embedding dimension - `num_stimuli`: number of stimuli - `sigma`: noise scale	required
`start_points`	`array - like`	Initial coordinates for optimization. If not provided, MDS-derived starting points are used.	`None`

Returns:

Name	Type	Description
`tuple`		(coordinates, solution_ll) - coordinates: fitted hyperbolic coordinates - solution_ll: negative log-likelihood of the fitted solution

Source code in src/rs_py/model/fit_geometric_models.py

def hyperbolic_points_of_best_fit(judgments, number_repeats, args, start_points=None):
    """
        Fit hyperbolic coordinates to the observed judgments.

        Args:
            judgments (dict):
                Pairwise comparison judgments.
            number_repeats (dict):
                Number of repeats for each comparison key in `judgments`.
            args (dict):
                Model settings. Must include:
                - `n_dim`: embedding dimension
                - `num_stimuli`: number of stimuli
                - `sigma`: noise scale
            start_points (array-like, optional):
                Initial coordinates for optimization. If not provided, MDS-derived
                starting points are used.

        Returns:
            tuple:
                (coordinates, solution_ll)
                - coordinates: fitted hyperbolic coordinates
                - solution_ll: negative log-likelihood of the fitted solution
    """
    def hyperbolic_cost(stimulus_params, pair_a, pair_b, counts, repeats, parameters):
        curvature = args['scripts']
        vectors = ac.params_to_points(stimulus_params, parameters['num_stimuli'], parameters['n_dim'])
        # mean center points
        vectors = vectors.T
        center = mean(vectors, 1)
        center = center.reshape((parameters['n_dim'], 1)) * ones((1, parameters['num_stimuli']))
        vectors = vectors - center
        # map Euclidean points to hyperboloid
        hyperboloid_points = loid_map(vectors, curvature)
        # compute distances
        distances = hyperbolic_distances(hyperboloid_points, curvature)
        probs = find_probabilities(distances, pair_a, pair_b, parameters['noise_st_dev'])
        # calculate log-likelihood, is_bad flag
        ll, is_bad = calculate_ll(counts, probs, repeats)
        LOG.debug('geometry is good: {}'.format(not is_bad))
        # fmin_costs.append(-1 * ll)  # debugging fmin
        return -1 * ll

    # calculate noise before continuing
    args['noise_st_dev'] = args['sigma']
    if start_points is None:
        # if not specified start minimization at coordiates returned by MDS after calculation of win-loss distances
        start_0 = mds.get_coordinates(args['n_dim'], judgments, number_repeats)[0]
    else:
        start_0 = start_points
    start = ac.anchor_points(start_0)

    LOG.debug("########  Procrustes distance between start and anchored start: {}".format(
        procrustes(start, start_0)[2]))
    # turn points to params
    start_params = ac.points_to_params(start)
    LOG.info('######## Run minimization on MDS start points (scipy minimize)')

    pairs_a, pairs_b, response_counts, comp_repeats = util.judgments_to_arrays(judgments, number_repeats)
    solution = gradient_descent(hyperbolic_cost, start_params, pairs_a, pairs_b, response_counts, comp_repeats, args)
    solution_ll = hyperbolic_cost(solution, pairs_a, pairs_b, response_counts, comp_repeats, args)

    coordinates = ac.params_to_points(solution, args['num_stimuli'], args['n_dim'])
    LOG.debug('########  Procrustes distance between anchored start and final solution: {}'.format(
        procrustes(start, coordinates)[2])
    )
    return coordinates, solution_ll  # , sum_residual_squares

`points_of_best_fit(judgments, number_repeats, args, start_points=None)`

Fit Euclidean coordinates to the observed judgments.

Parameters:

Name	Type	Description	Default
`judgments`	`dict`	Pairwise comparison judgments	required
`number_repeats`	`dict`	Number of repeats for each comparison key in `judgments`.	required
`args`	`dict`	Model settings. Must include: - `n_dim`: embedding dimension - `num_stimuli`: number of stimuli - `sigma`: noise scale - `tolerance`: stopping tolerance for optimization - `minimization`: optimization method, either `gradient-descent` or `nelder-mead`	required
`start_points`	`array - like`	Initial coordinates for optimization. If not provided, MDS-derived starting points are used.	`None`

Returns:

Name	Type	Description
`tuple`		(coordinates, solution_ll) - coordinates: fitted Euclidean coordinates - solution_ll: negative log-likelihood of the fitted solution

Source code in src/rs_py/model/fit_geometric_models.py

def points_of_best_fit(judgments, number_repeats, args, start_points=None):
    """
    Fit Euclidean coordinates to the observed judgments.

    Args:
        judgments (dict):
            Pairwise comparison judgments
        number_repeats (dict):
            Number of repeats for each comparison key in `judgments`.
        args (dict):
            Model settings. Must include:
            - `n_dim`: embedding dimension
            - `num_stimuli`: number of stimuli
            - `sigma`: noise scale
            - `tolerance`: stopping tolerance for optimization
            - `minimization`: optimization method, either `gradient-descent`
              or `nelder-mead`
        start_points (array-like, optional):
            Initial coordinates for optimization. If not provided, MDS-derived
            starting points are used.

    Returns:
        tuple:
            (coordinates, solution_ll)
            - coordinates: fitted Euclidean coordinates
            - solution_ll: negative log-likelihood of the fitted solution
    """

    def cost(stimulus_params, pair_a, pair_b, counts, repeats, parameters):
        vectors = ac.params_to_points(stimulus_params, parameters['num_stimuli'], parameters['n_dim'])
        ll, is_bad = dist_model_ll_vectorized(pair_a, pair_b, counts, repeats, parameters, vectors)
        LOG.debug('geometry is good: {}'.format(not is_bad))
        return -1 * ll

    # calculate noise before continuing
    args['noise_st_dev'] = args['sigma']
    if start_points is None:
        # if not specified start minimization at coordiates returned by MDS after calculation of win-loss distances
        start_0 = mds.get_coordinates(args['n_dim'], judgments, number_repeats)[0]
    else:
        start_0 = start_points
    start = ac.anchor_points(start_0)
    LOG.info("########  Procrustes distance between start and anchored start: {}".format(
        procrustes(start, start_0)[2]))
    # turn points to params
    start_params = ac.points_to_params(start)
    LOG.info('######## Run minimization on MDS start points (scipy minimize)')
    # make maxiter 60000 for 5D geometry
    options_min = {
        'disp': True,
        'fatol': args['tolerance']}
    if args['n_dim'] < 4:
        options_min['maxiter'] = 85000
    elif args['n_dim'] >= 4:
        options_min['maxiter'] = 110000
    pairs_a, pairs_b, response_counts, comp_repeats = util.judgments_to_arrays(judgments, number_repeats)
    if args['minimization'] == 'nelder-mead':
        optimal = optimize.minimize(cost, start_params,
                                    args=(pairs_a, pairs_b, response_counts, comp_repeats, args),
                                    method='Nelder-Mead',
                                    options=options_min
                                    )
        LOG.info(
            '######## {} Iterations completed: {}. Model Dim: {}.'.format(optimal.message, optimal.nit,
                                                                          args['n_dim'])
        )
        solution = optimal.model_coords
        solution_ll = optimal.fun
    else:
        solution = gradient_descent(cost, start_params, pairs_a, pairs_b, response_counts, comp_repeats, args)
        stim = ac.params_to_points(solution, args['num_stimuli'], args['n_dim'])
        ll_final, is_model_bad = dist_model_ll_vectorized(pairs_a, pairs_b, response_counts, comp_repeats, args, stim)
        solution_ll = -1 * ll_final
        LOG.debug("Final Model is good/ feasible: {}".format(not is_model_bad))

    coordinates = ac.params_to_points(solution, args['num_stimuli'], args['n_dim'])

    try:
        procr_dist = procrustes(start, coordinates)[2]
    except ValueError:
        procr_dist = 'WARNING - problem with the coordinates. Nans or infs possible.'
    LOG.debug('########  Procrustes distance between anchored start and final solution: {}'.format(
        procr_dist)
    )
    return coordinates, solution_ll  # , sum_residual_squares

`spherical_points_of_best_fit(judgments, number_repeats, args, start_points=None)`

Fit spherical coordinates to the observed judgments.

Parameters:

Name	Type	Description	Default
`judgments`	`dict`	Pairwise comparison judgments	required
`number_repeats`	`dict`	Number of repeats for each comparison key in `judgments`.	required
`args`	`dict`	Model settings. Must include at least: - `n_dim`: embedding dimension - `num_stimuli`: number of stimuli - `sigma`: noise scale	required
`start_points`	`array - like`	Initial coordinates for optimization. If not provided, MDS-derived starting points are used.	`None`

Returns:

Name	Type	Description
`tuple`		(coordinates, solution_ll) - coordinates: fitted spherical coordinates - solution_ll: negative log-likelihood of the fitted solution

Source code in src/rs_py/model/fit_geometric_models.py

def spherical_points_of_best_fit(judgments, number_repeats, args, start_points=None):
    """
    Fit spherical coordinates to the observed judgments.

    Args:
        judgments (dict):
            Pairwise comparison judgments
        number_repeats (dict):
            Number of repeats for each comparison key in `judgments`.
        args (dict):
            Model settings. Must include at least:
            - `n_dim`: embedding dimension
            - `num_stimuli`: number of stimuli
            - `sigma`: noise scale
        start_points (array-like, optional):
            Initial coordinates for optimization. If not provided, MDS-derived
            starting points are used.

    Returns:
        tuple:
            (coordinates, solution_ll)
            - coordinates: fitted spherical coordinates
            - solution_ll: negative log-likelihood of the fitted solution

    """
    # for debugging
    fmin_costs = []

    def spherical_cost(stimulus_params, pair_a, pair_b, counts, repeats, parameters):
        curvature = args['scripts']
        vectors = ac.params_to_points(stimulus_params, parameters['num_stimuli'], parameters['n_dim'])
        # mean center points
        vectors = vectors.T
        center = mean(vectors, 1)
        center = center.reshape((parameters['n_dim'], 1)) * ones((1, parameters['num_stimuli']))
        vectors = vectors - center
        # map Euclidean points to sphere
        sph_points = sphere_map(vectors, 1/curvature)
        # compute distances
        distances = spherical_distances(sph_points, 1/curvature)
        probs = find_probabilities(distances, pair_a, pair_b, parameters['noise_st_dev'])
        # calculate log-likelihood, is_bad flag
        ll, is_bad = calculate_ll(counts, probs, repeats)
        LOG.debug('geometry is good: {}'.format(not is_bad))
        # fmin_costs.append(-1 * ll)  # debugging fmin
        return -1 * ll

    # calculate noise before continuing
    args['noise_st_dev'] = args['sigma']
    if start_points is None:
        # if not specified start minimization at coordiates returned by MDS after calculation of win-loss distances
        start_0 = mds.get_coordinates(args['n_dim'], judgments, number_repeats)[0]
    else:
        start_0 = start_points
    start = ac.anchor_points(start_0)

    LOG.debug("########  Procrustes distance between start and anchored start: {}".format(
        procrustes(start, start_0)[2]))
    # turn points to params
    start_params = ac.points_to_params(start)
    LOG.info('######## Run minimization on MDS start points (scipy minimize)')

    pairs_a, pairs_b, response_counts, comp_repeats = util.judgments_to_arrays(judgments, number_repeats)
    solution = gradient_descent(spherical_cost, start_params, pairs_a, pairs_b, response_counts, comp_repeats, args)
    solution_ll = spherical_cost(solution, pairs_a, pairs_b, response_counts, comp_repeats, args)

    coordinates = ac.params_to_points(solution, args['num_stimuli'], args['n_dim'])
    try:
        procr_dist = procrustes(start, coordinates)[2]
    except ValueError:
        procr_dist = 'WARNING - problem with the coordinates. Nans or infs possible.'
    LOG.debug('########  Procrustes distance between anchored start and final solution: {}'.format(
        procr_dist)
    )
    return coordinates, solution_ll  # , sum_residual_squares

`scripts`

`model_fitting`

Script wrapper for the geometric model fit demo.

`build_model_fit_args(user_params)`

Parameters:

Name	Type	Description	Default
`user_params`	`dict`	Partial parameter dictionary supplied by the caller. Any keys present here override the defaults in CONFIG["inputs"][]. If provided, this may include required keys such as "filepath" and "output_dir", along with optional settings like "sigma" or "model_dimensions".	required

Returns: args (dict): A complete parameter dictionary for the model-fit runner. This is built from CONFIG["inputs"][] plus any overrides from user_params.

Notes

The returned dict is intended to be passed directly to the run function. Required keys must still be present after merging.

Source code in src/rs_py/scripts/model_fitting.py

def build_model_fit_args(user_params):
    """
    Args:
        user_params (dict):
            Partial parameter dictionary supplied by the caller.
            Any keys present here override the defaults in CONFIG["inputs"]["model_fit"].
            If provided, this may include required keys such as "filepath" and "output_dir",
            along with optional settings like "sigma" or "model_dimensions".
    Returns:
        args (dict):
        A complete parameter dictionary for the model-fit runner.
        This is built from CONFIG["inputs"]["model_fit"] plus any overrides from user_params.

    Notes:
        The returned dict is intended to be passed directly to the run function.
         Required keys must still be present after merging.
    """
    params = deepcopy(CONFIG["inputs"]["model_fit"])

    # Merge a partial user dict over the defaults.
    if user_params:
        params.update(user_params)

    missing = [k for k in REQUIRED_KEYS if not params.get(k)]
    if missing:
        raise ValueError(f"Missing required parameter(s): {', '.join(missing)}")

    args = {
        "filepath": params["filepath"],
        "outdir": params["output_dir"],
        "exp_name": params["exp_name"],
        "subject": params["subject"],
        "sigma": params["sigma"],
        "model_dimensions": params["model_dimensions"],
        "learning_rate": params["learning_rate"],
        "tolerance": params["tolerance"],
        "max_iterations": params["max_iterations"],
        "minimization": params["minimization"],
        "filter_trials": params["max_trials"],
    }

    # Convenience alias used by the fitting code.
    args["noise_st_dev"] = args["sigma"]

    return args

`fit(params)`

Fit Euclidean geometric models to aggregated choice data and save model outputs.

This function loads a combined choice file, optionally subsamples a fixed number of unique triadic judgments, fits Euclidean models across one or more requested dimensions, compares them against best-choice and random-choice baselines, and saves both a CSV summary and a .mat file of fitted coordinates.

Parameters:

Name	Type	Description	Default
`params`	`dict`	Dictionary of run parameters. Common keys include: Required: filepath (str): Path to the combined choices `.mat` file. outdir (str): Directory where output files will be written. Recommended / commonly used: exp_name (str): Experiment name used in output filenames. subject (str): Subject ID used in output filenames. model_dimensions (list[int]): Dimensions to fit, for example `[1, 2, 3, 4, 5]`. The default is typically a small set of low-dimensional models, but this can be changed depending on the number of stimuli and the modeling goal. Tested values may extend beyond the default range. sigma (float): Noise standard deviation used in the choice model. This is a required scale parameter and is usually kept near 1. learning_rate (float): Learning rate for optimization. tolerance (float): Convergence tolerance used by the optimizer. This is most relevant for gradient-descent-based fitting. max_iterations (int): Maximum number of optimization iterations. minimization (str): Optimization method. Supported values include `gradient-descent` and `nelder-mead`. Gradient descent is typically faster and the recommended default. filter_trials (int): Maximum number of unique triadic judgments to use. Set to 0 or omit to use all available data. If this is smaller than the total number of unique comparisons, the function samples a random subset.	required

Returns:

Name	Type	Description
`None`		Writes output files to `outdir` and prints a summary to the console.

Raises:

Type	Description
`FileNotFoundError`	If the input file does not exist.
`ValueError`	If no triadic judgments are found after loading or filtering.

Notes

The input file is expected to be a combined choice file produced by the earlier aggregation step.
sigma is a required noise scale and controls the uncertainty in distance comparisons.
model_dimensions should be chosen with the number of stimuli and the intended model complexity in mind.
The function saves fitted coordinates for each dimension, along with a CSV comparing model log-likelihoods.

Source code in src/rs_py/scripts/model_fitting.py

def fit(params):
    """
    Fit Euclidean geometric models to aggregated choice data and save model outputs.

    This function loads a combined choice file, optionally subsamples a fixed
    number of unique triadic judgments, fits Euclidean models across one or more
    requested dimensions, compares them against best-choice and random-choice
    baselines, and saves both a CSV summary and a `.mat` file of fitted
    coordinates.

    Args:
        params (dict):
            Dictionary of run parameters. Common keys include:

            Required:
                filepath (str):
                    Path to the combined choices `.mat` file.
                outdir (str):
                    Directory where output files will be written.

            Recommended / commonly used:
                exp_name (str):
                    Experiment name used in output filenames.
                subject (str):
                    Subject ID used in output filenames.
                model_dimensions (list[int]):
                    Dimensions to fit, for example `[1, 2, 3, 4, 5]`.
                    The default is typically a small set of low-dimensional models,
                    but this can be changed depending on the number of stimuli and
                    the modeling goal. Tested values may extend beyond the default
                    range.
                sigma (float):
                    Noise standard deviation used in the choice model. This is a
                    required scale parameter and is usually kept near 1.
                learning_rate (float):
                    Learning rate for optimization.
                tolerance (float):
                    Convergence tolerance used by the optimizer. This is most
                    relevant for gradient-descent-based fitting.
                max_iterations (int):
                    Maximum number of optimization iterations.
                minimization (str):
                    Optimization method. Supported values include
                    `gradient-descent` and `nelder-mead`.
                    Gradient descent is typically faster and the recommended
                    default.
                filter_trials (int):
                    Maximum number of unique triadic judgments to use. Set to 0
                    or omit to use all available data. If this is smaller than the
                    total number of unique comparisons, the function samples a
                    random subset.

    Returns:
        None:
            Writes output files to `outdir` and prints a summary to the console.

    Raises:
        FileNotFoundError:
            If the input file does not exist.
        ValueError:
            If no triadic judgments are found after loading or filtering.

    Notes:
        - The input file is expected to be a combined choice file produced by the
          earlier aggregation step.
        - `sigma` is a required noise scale and controls the uncertainty in
          distance comparisons.
        - `model_dimensions` should be chosen with the number of stimuli and the
          intended model complexity in mind.
        - The function saves fitted coordinates for each dimension, along with a
          CSV comparing model log-likelihoods.
    """
    print("\n" + "=" * 70)
    print("GEOMETRIC MODEL FIT")
    print("=" * 70)

    print("DATA")
    print("-" * 70)
    print(f"Filepath:            {params['filepath']}")
    print(f"Experiment:          {params['exp_name']}")
    print(f"Subject:             {params['subject']}")
    print(f"Output directory:    {params['outdir']}")
    print(f"Max trials used:     {params['filter_trials']}")

    print("\nOPTIMIZATION SETTINGS")
    print("-" * 70)
    print(f"Max iterations:      {params['max_iterations']}")
    print(f"Learning rate:       {params['learning_rate']}")
    print(f"Tolerance:           {params['tolerance']}")

    print("\nNOISE PARAMETERS")
    print("-" * 70)
    print(f"Sigma (compare):     {params['sigma']:.6f}")

    print("=" * 70)

    # break up ranking responses into pairwise judgments
    pairwise_responses, pairwise_num_repeats, metadata, stim_list = read_combined_choices(params['filepath'])
    params["num_stimuli"] = len(stim_list)

    print("\nLoaded pairwise judgments")
    print("-" * 60)
    print(f"Number of unique comparisons: {len(pairwise_responses)}")
    print(f"Total triads (including repeats): {sum(pairwise_num_repeats.values())}")
    print("-" * 60)

    # only consider a subset of trials
    subset = {}
    if params['filter_trials']:
        if params['filter_trials'] < len(pairwise_responses):
            all_keys = list(pairwise_responses.keys())
            chosen = random.sample(all_keys, params['filter_trials'])
            subset = {k: pairwise_responses[k] for k in chosen}

            print("\nUsing subset of trials")
            print("-" * 60)
            print(f"Trials used: {len(subset)}")
            print(f"Total triads used: {sum(pairwise_num_repeats[k] for k in subset)}")
            print("=" * 60)
        else:
            subset = pairwise_responses
    else:
        subset = pairwise_responses

    # initialize results dataframe
    total_num_triads = sum([pairwise_num_repeats[k] for k in subset.keys()])
    if total_num_triads == 0:
        raise ValueError("No triads found (total_num_triads == 0). Check input data or arguments.")

    # initialize results
    result = {
        "Model": [],
        "Log Likelihood": [],
        "number of points": [],
        "Experiment": [],
        "Subject": []
    }

    # Euclidean models across dimensions
    coords_by_dim = {}
    lls_by_dim = {}

    dimensions = params["model_dimensions"]
    n_workers = min(os.cpu_count() or 1, len(dimensions))

    with mp.get_context("spawn").Pool(n_workers) as pool:
        fitted_models = pool.starmap(
            fit_one_dimension,
            [
                (dim, subset, pairwise_num_repeats, params)
                for dim in dimensions
            ],
        )

    for dim, model_coords, ll_nd in fitted_models:
        print("\n" + "=" * 60)
        print(f"FITTED {dim}D EUCLIDEAN MODEL")
        print("=" * 60)

        model_name = f"{dim}D"

        # (Distances computed previously; keep if useful for debugging)
        _ = pdist(model_coords)
        coords_by_dim[dim] = model_coords
        lls_by_dim[dim] = ll_nd

        print("\nOptimized embedding:")
        print("-" * 60)
        print(f"Shape: {model_coords.shape}")
        print(f"Mean coordinate value: {np.mean(model_coords):.4f}")
        print(f"Std of coordinates:    {np.std(model_coords):.4f}")
        print(f"Min/Max coordinate:    {np.min(model_coords):.4f} / {np.max(model_coords):.4f}")
        print("-" * 60)

        outfilename = '{}/{}_{}_anchored_points_sigma_{}_dim_{}'.format(
            params['outdir'], params['subject'], params['exp_name'], params['sigma'], dim
        )
        np.save(outfilename, model_coords)

        ll_nd = -ll_nd / float(total_num_triads)
        print(f"Negative Log Likelihood per triad of the model: {ll_nd:.4f}")

        result['Model'].append(model_name)
        result['Log Likelihood'].append(ll_nd)
        result['number of points'].append(params['num_stimuli'])

    # ---- Best and random baselines ----
    ll_best = an.best_model_ll(
        subset, pairwise_num_repeats)[0] / float(total_num_triads)
    result['Model'].append('best')
    result['Log Likelihood'].append(ll_best)
    result['number of points'].append(params['num_stimuli'])

    ll_random = an.random_choice_ll(
        subset, pairwise_num_repeats)[0] / float(total_num_triads)
    result['Model'].append('random')
    result['Log Likelihood'].append(ll_random)
    result['number of points'].append(params['num_stimuli'])

    print("\n" + "=" * 60)
    print("BASELINE COMPARISON")
    print("=" * 60)
    print(f"Best possible model LL:   {ll_best:.4f}")
    print(f"Random choice model LL:   {ll_random:.4f}")
    print("=" * 60)

    # Fill Experiment/Subject columns to match row count
    n_rows = len(result["Model"])
    result["Experiment"] = [params['exp_name']] * n_rows
    result["Subject"] = [params['subject']] * n_rows

    # ---- Output results ----
    data_frame = pd.DataFrame(result)

    print("\nFINAL MODEL COMPARISON")
    print("=" * 60)
    print(data_frame.to_string(index=False))
    print("=" * 60)

    data_frame.to_csv('{}/{}-{}-geometry-likelihoods_with_{}_trials_sigma_{}_{}_pts_anchored.csv'
                      .format(params['outdir'],
                              params['subject'],
                              params['exp_name'],
                              'all' if not params['filter_trials'] else params['filter_trials'],
                              params['sigma'],
                              params['num_stimuli']
                              ), index=False)

    # write combined file with coords and lls
    stimuli = stim_list
    lls_by_dim['best'] = ll_best
    lls_by_dim['random'] = ll_random

    mat_path = create_coords_file(
        outdir=params['outdir'],
        exp=params['exp_name'],
        subject=params['subject'],
        model_dimensions=params['model_dimensions'],
        points=coords_by_dim,
        lls=lls_by_dim,
        stim_labels=stimuli
    )
    print("Saved:", mat_path)

`write_choice_file`

Script wrapper for the combined-choice demo.

`write_choice_file_detailed`

To demo the script choice_file_detailed.py

`run(user_params)`

Run the detailed choice-file builder on a set of user parameters.

Parameters:

Name	Type	Description	Default
`user_params`	`dict`	Dictionary of parameters used to build the detailed choice file. Required keys: - input_path: Path to the input directory containing raw data. - output_dir: Path to the directory where outputs should be saved. Optional key: - metadata: Dictionary with optional fields such as: - exp_name: Experiment name. - subject: Subject identifier. - stim_list: List of stimuli. - num_sessions: Number of sessions. - num_trials: Number of trials. - total_judgments: Total number of judgments. - judgment_type: Type of judgment task, default is "triadic".	required

Returns:

Type	Description
	None

Raises: - ValueError: If a required parameter is missing or empty. - TypeError: If metadata is provided but is not a dictionary. - FileNotFoundError: If input_path does not exist.

Source code in src/rs_py/scripts/write_choice_file_detailed.py

def run(user_params):
    """
    Run the detailed choice-file builder on a set of user parameters.

    Args:
        user_params (dict): Dictionary of parameters used to build the detailed choice
            file.
            Required keys:
                - input_path: Path to the input directory containing raw data.
                - output_dir: Path to the directory where outputs should be saved.
            Optional key:
                - metadata: Dictionary with optional fields such as:
                - exp_name: Experiment name.
                - subject: Subject identifier.
                - stim_list: List of stimuli.
                - num_sessions: Number of sessions.
                - num_trials: Number of trials.
                - total_judgments: Total number of judgments.
                - judgment_type: Type of judgment task, default is "triadic".

    Returns:
        None

    Raises:
    - ValueError: If a required parameter is missing or empty.
    - TypeError: If metadata is provided but is not a dictionary.
    - FileNotFoundError: If input_path does not exist.
    """
    params = merge_with_defaults(user_params)
    validate_required(params)

    input_path = params["input_path"]
    output_dir = params["output_dir"]

    metadata = params["metadata"]
    exp_name = metadata["exp_name"]
    subject = metadata["subject"]

    if not os.path.exists(input_path):
        raise FileNotFoundError(f"Input directory not found: {input_path}")

    os.makedirs(output_dir, exist_ok=True)

    print("\nProcessing raw data...")
    print(f"  Input directory: {input_path}")
    print(f"  Output directory: {output_dir}")
    print(f"  Subject: {subject}")
    print(f"  Experiment: {exp_name}")
    print(f"  Types of judgments: {metadata['judgment_type']}\n")

    cfd.build_detailed_choice_mat(
        input_dir=input_path,
        output_dir=output_dir,
        exp_name=exp_name,
        subject=subject,
        metadata=metadata,
    )

    print("\nDone.")

`utils`

`anchor_coordinates`

Use the Gram Schmidt procedure to reorient a set of points such that one point is at the origin, another point lies on one axis (x1, 0, ..., 0) and the next point lies on a plane spanned by the two axes (x2, x3, 0, ..., 0), and so on. This constrains the minimization function that looks for the points of best fit, as rotations and translations are made redundant.

This can also take into effect the noise parameter as we do not constrain where on the first axis, the second point lies.

`anchor_points(points)`

points: 37 by d matrix Given a set of 37 points, find the rotation needed using Gram Schmidt and then transform the coordinates to the format needed - one point is zero and the rest are partly in a triangular form