Hi Pymc Folks,

I’m working on simulating audio where multiple bats may be emitting sounds - and trying to see which subset of possible scenarios matches my observed audio data best (e.g. with 3 individuals A,B,C we can have 7 scenarios/combinations where at least one of them emits a sound) . I’m doing this using ABC as I need to understand the match to my audio data itself.

My simulation function takes two parameters per draw - an integer which decides the scenario to be simulated. This integer is used to subset a binary ‘scenario array’. The other input parameter is a random draw of potential positions from where sound was emitted.

Before I even began to write the audio simulation portion - I ran into inconsistencies on how the drawn integer was treated - it seems like pymc is treating it like an np.ndarray when I print its `type()`

, but throws an IndexError when I try to use the integer to subset the scenario array.

I understand from looking around the docs and here that `pytensor`

and `numpy`

objects don’t mix well. Does anyone have any ideas on 1) why this odd type inconsistency occurs (the output type is `numpy.ndarray`

- but it is not really this kind of object?), and 2) if there is anyway to convert a pytensor into a numpy array so that I can perform all operations within the simulation only with numpy objects.

System specs: Windows 10, Python 3.11, Pymc 5.5.0. Same error also seen in Python 3.11, Pymc 5.5.0 and Ubuntu 20.04 too.

```
# -*- coding: utf-8 -*-
"""
Created on Tue Jul 4 18:22:40 2023
@author: theja
"""
from itertools import combinations
import numpy as np
import pymc as pm
import os
os.environ['PYTHONIOENCODING']='utf-8'
def generate_all_possible_scenarios(nbats):
'''
Generates all possible 2^nbats - 1 possible combinations where
at least 1 bat is calling
'''
all_scenarios = []
scenario_number = 0
for subset in range(1, nbats+1):
subscenario_combis = combinations(range(1,nbats+1),subset)
for each in subscenario_combis:
binary_scenario = np.zeros(nbats)
batinds = np.array((each))-1
binary_scenario[batinds] = 1
all_scenarios.append(binary_scenario )
scenario_number += 1
all_scenarios = np.array(all_scenarios)
return all_scenarios
nbats = 5
# generate various possible scenarios individual called (1=call, 0=silent)
allscenarios_asnp = generate_all_possible_scenarios(nbats)
# set probability of all scenarios to be the same
num_scenarios = allscenarios_asnp.shape[0]
p_categories = np.tile(1/num_scenarios, num_scenarios)
p_categories = p_categories/p_categories.sum() # 're' normalise everything to be sure it adds to 1
# x,y,z ranges of possible emission points for each individual per row
# ordered as xmin, xmax, ymin, ymax, zmin, zmax
xyz_ranges = np.ones((nbats,6))
for batnum in range(nbats):
xyz_ranges[batnum,:] = batnum + np.linspace(0,1,6)
def generate_fake_cc(rng, scenario_number, xyz_emissions, size=None):
'''
An incomplete but functional example to show issues with subsetting
Parameters
----------
rng :
scenario_number : int? or pytensor?
The row number to choose from the binary scenario matrix
xyz_emissions : np.array? or pytensor?
The set of drawn potential emission points
size :
'''
print(type(xyz_emissions), xyz_emissions.shape,
type(scenario_number), scenario_number)
# PROBLEM POINT when uncommented
# calling_bats = allscenarios_asnp[scenario_number,:]
# dummy output just to get the function to run
return np.array([5, 2]).reshape(-1,2)
if __name__ == "__main__":
with pm.Model() as mo:
# choose which scenario is being tested
current_scenario = pm.Categorical('current_scenario',
p_categories,
shape=(1,))
# generate the hypothesised x,y,z for the emitting bats
x_hyp = pm.Uniform('x_hyp', lower=xyz_ranges[:,0],
upper=xyz_ranges[:,1], shape=(nbats,))
y_hyp = pm.Uniform('y_hyp', lower=xyz_ranges[:,2],
upper=xyz_ranges[:,3],shape=(nbats,))
z_hyp = pm.Uniform('z_hyp', lower=xyz_ranges[:,4],
upper=xyz_ranges[:,5],shape=(nbats,))
xyz_stack = pm.math.stack([x_hyp, y_hyp, z_hyp], axis=1)
print('xyz_stack:',xyz_stack.shape)
# #scenario_data[1] = xyz_stack
# # recreate the simulated audio and cross-correlations
sim = pm.Simulator('sim', generate_fake_cc, params=(current_scenario, xyz_stack),
observed=np.array([3,3]).reshape(-1,2),
)
with mo:
idata = pm.sample_smc(draws=100, chains=2)
idata.extend(pm.sample_posterior_predictive(idata))
```