Aesara fails to compile graph for model that runs on pymc3

Questions version agnostic
1

Hi folks, I hope you all are doing well!

I am running into an odd issue as I am trying to run a Gaussian Processes model in pymc v4 roughly following the steps outlined in @AlexAndorra presidential poll case study. Below is the actual code I am using.

months = dates_to_idx(df_train_reg["ds"]) # transform dates in float indicating months passed, shape = (110,)
coords = {"time": df_train_reg["ds"]}

with pm.Model(coords=coords) as gp_demand_model:
    
    baseline = pm.Normal("baseline", mu=0, sigma=0.11)
    no_homeless = pm.Normal("no_homeless", mu=0.2, sigma=0.1)
    covid = pm.Normal("covid", mu=-0.3, sigma=0.1)
    violence = pm.Normal("violence", mu=0.1, sigma=0.1)
    
    # time trends
    amplitude_trend = pm.HalfNormal("amplitude_trend", sigma=3)
    ls_trend = pm.Gamma("ls_trend", alpha=5, beta=2)
    cov_trend = amplitude_trend ** 2 * pm.gp.cov.Matern52(1, ls_trend)
    
    gp = pm.gp.Latent(cov_func=cov_trend)
    f_time = gp.prior("f_time", X=months[:, None])
    
    # data
    homeless_data = pm.Data("homeless_data", df_train_std["no_homeless"].values, dims="time")
    covid_data = pm.Data("covid_data", df_train_std["new_covid_cases"].values, dims="time")
    violence_data = pm.Data("violence_data", df_train_std["viol_events"].values, dims="time")
    
    dem_mu = pm.Deterministic(
        "dem_mu",
        var = (
            baseline 
            + f_time 
            + no_homeless * homeless_data 
            + covid * covid_data
            + violence * violence_data
        ),
        dims="time"
    )
    dem_sigma = pm.HalfNormal("dem_sigma", sigma=0.2)
    
    # Likelihood
    demand = pm.Normal(
        "demand", 
        mu=dem_mu, 
        sigma=dem_sigma, 
        observed=df_train_std["y"].values, 
        dims="time"
    )
    
    priors = pm.sample_prior_predictive()

Below is the graph of the model from graphviz:

image
image1060Ă—559 22.6 KB

If I comment out the line where I define f_time (and the relative line in the definition of dem_mu) the code runs without problem (note that here I am only generating prior predictive samples). When those lines are not commented I instead get the following error:

---------------------------------------------------------------------------
CompileError                              Traceback (most recent call last)
File ~\anaconda3\envs\test\lib\site-packages\aesara\link\vm.py:1246, in VMLinker.make_all(self, profiler, input_storage, output_storage, storage_map)
   1242 # no-recycling is done at each VM.__call__ So there is
   1243 # no need to cause duplicate c code by passing
   1244 # no_recycling here.
   1245 thunks.append(
-> 1246     node.op.make_thunk(node, storage_map, compute_map, [], impl=impl)
   1247 )
   1248 linker_make_thunk_time[node] = time.time() - thunk_start

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\c\op.py:131, in COp.make_thunk(self, node, storage_map, compute_map, no_recycling, impl)
    130 try:
--> 131     return self.make_c_thunk(node, storage_map, compute_map, no_recycling)
    132 except (NotImplementedError, MethodNotDefined):
    133     # We requested the c code, so don't catch the error.

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\c\op.py:96, in COp.make_c_thunk(self, node, storage_map, compute_map, no_recycling)
     95         raise NotImplementedError("float16")
---> 96 outputs = cl.make_thunk(
     97     input_storage=node_input_storage, output_storage=node_output_storage
     98 )
     99 thunk, node_input_filters, node_output_filters = outputs

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\c\basic.py:1198, in CLinker.make_thunk(self, input_storage, output_storage, storage_map)
   1197 init_tasks, tasks = self.get_init_tasks()
-> 1198 cthunk, module, in_storage, out_storage, error_storage = self.__compile__(
   1199     input_storage, output_storage, storage_map
   1200 )
   1202 res = _CThunk(cthunk, init_tasks, tasks, error_storage, module)

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\c\basic.py:1133, in CLinker.__compile__(self, input_storage, output_storage, storage_map)
   1132 output_storage = tuple(output_storage)
-> 1133 thunk, module = self.cthunk_factory(
   1134     error_storage,
   1135     input_storage,
   1136     output_storage,
   1137     storage_map,
   1138 )
   1139 return (
   1140     thunk,
   1141     module,
   (...)
   1150     error_storage,
   1151 )

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\c\basic.py:1629, in CLinker.cthunk_factory(self, error_storage, in_storage, out_storage, storage_map)
   1628         node.op.prepare_node(node, storage_map, None, "c")
-> 1629     module = get_module_cache().module_from_key(key=key, lnk=self)
   1631 vars = self.inputs + self.outputs + self.orphans

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\c\cmodule.py:1228, in ModuleCache.module_from_key(self, key, lnk)
   1227 location = dlimport_workdir(self.dirname)
-> 1228 module = lnk.compile_cmodule(location)
   1229 name = module.__file__

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\c\basic.py:1538, in CLinker.compile_cmodule(self, location)
   1537     _logger.debug(f"LOCATION {location}")
-> 1538     module = c_compiler.compile_str(
   1539         module_name=mod.code_hash,
   1540         src_code=src_code,
   1541         location=location,
   1542         include_dirs=self.header_dirs(),
   1543         lib_dirs=self.lib_dirs(),
   1544         libs=libs,
   1545         preargs=preargs,
   1546     )
   1547 except Exception as e:

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\c\cmodule.py:2639, in GCC_compiler.compile_str(module_name, src_code, location, include_dirs, lib_dirs, libs, preargs, py_module, hide_symbols)
   2635     # We replace '\n' by '. ' in the error message because when Python
   2636     # prints the exception, having '\n' in the text makes it more
   2637     # difficult to read.
   2638     # compile_stderr = compile_stderr.replace("\n", ". ")
-> 2639     raise CompileError(
   2640         f"Compilation failed (return status={status}):\n{' '.join(cmd)}\n{compile_stderr}"
   2641     )
   2642 elif config.cmodule__compilation_warning and compile_stderr:
   2643     # Print errors just below the command line.

CompileError: Compilation failed (return status=1):
"C:\Users\andri\anaconda3\envs\test_pymc3\Library\mingw-w64\bin\g++.exe" -shared -g -O3 -fno-math-errno -Wno-unused-label -Wno-unused-variable -Wno-write-strings -Wno-c++11-narrowing -fno-exceptions -fno-unwind-tables -fno-asynchronous-unwind-tables -march=broadwell -mmmx -mno-3dnow -msse -msse2 -msse3 -mssse3 -mno-sse4a -mcx16 -msahf -mmovbe -maes -mno-sha -mpclmul -mpopcnt -mabm -mno-lwp -mfma -mno-fma4 -mno-xop -mbmi -mbmi2 -mno-tbm -mavx -mavx2 -msse4.2 -msse4.1 -mlzcnt -mno-rtm -mno-hle -mrdrnd -mf16c -mfsgsbase -mrdseed -mprfchw -madx -mfxsr -mxsave -mxsaveopt -mno-avx512f -mno-avx512er -mno-avx512cd -mno-avx512pf -mno-prefetchwt1 -mclflushopt -mxsavec -mxsaves -mno-avx512dq -mno-avx512bw -mno-avx512vl -mno-avx512ifma -mno-avx512vbmi -mno-clwb -mno-pcommit -mno-mwaitx --param l1-cache-size=32 --param l1-cache-line-size=64 --param l2-cache-size=12288 -mtune=generic -DNPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION -m64 -DMS_WIN64 -I"C:\Users\andri\anaconda3\envs\test\lib\site-packages\numpy\core\include" -I"C:\Users\andri\anaconda3\envs\test\include" -I"C:\Users\andri\anaconda3\envs\test\lib\site-packages\aesara\link\c\c_code" -L"C:\Users\andri\anaconda3\envs\test\libs" -L"C:\Users\andri\anaconda3\envs\test" -o "C:\Users\andri\AppData\Local\Aesara\compiledir_Windows-10-10.0.22000-SP0-Intel64_Family_6_Model_158_Stepping_10_GenuineIntel-3.8.13-64\tmp11_4bwf_\m47260be5189a2297f95a5722358fab6ed80907bc9b3b7adb3279eddf44b57064.pyd" "C:\Users\andri\AppData\Local\Aesara\compiledir_Windows-10-10.0.22000-SP0-Intel64_Family_6_Model_158_Stepping_10_GenuineIntel-3.8.13-64\tmp11_4bwf_\mod.cpp" -lblas "C:\Users\andri\anaconda3\envs\test\python38.dll"
C:/Users/andri/anaconda3/envs/test_pymc3/Library/mingw-w64/bin/../lib/gcc/x86_64-w64-mingw32/5.3.0/../../../../x86_64-w64-mingw32/bin/ld.exe: cannot find -lblas
collect2.exe: error: ld returned 1 exit status


During handling of the above exception, another exception occurred:

CompileError                              Traceback (most recent call last)
Input In [64], in <cell line: 4>()
     37 # Likelihood
     38 demand = pm.Normal(
     39     "demand", 
     40     mu=dem_mu, 
   (...)
     43     dims="time"
     44 )
---> 46 priors = pm.sample_prior_predictive()

File ~\anaconda3\envs\test\lib\site-packages\pymc\sampling.py:2242, in sample_prior_predictive(samples, model, var_names, random_seed, return_inferencedata, idata_kwargs, compile_kwargs)
   2239 compile_kwargs.setdefault("allow_input_downcast", True)
   2240 compile_kwargs.setdefault("accept_inplace", True)
-> 2242 sampler_fn = compile_forward_sampling_function(
   2243     vars_to_sample,
   2244     vars_in_trace=[],
   2245     basic_rvs=model.basic_RVs,
   2246     givens_dict=None,
   2247     random_seed=random_seed,
   2248     **compile_kwargs,
   2249 )
   2251 values = zip(*(sampler_fn() for i in range(samples)))
   2253 data = {k: np.stack(v) for k, v in zip(names, values)}

File ~\anaconda3\envs\test\lib\site-packages\pymc\sampling.py:1741, in compile_forward_sampling_function(outputs, vars_in_trace, basic_rvs, givens_dict, **kwargs)
   1730 # Populate the givens list
   1731 givens = [
   1732     (
   1733         node,
   (...)
   1738     for node, value in givens_dict.items()
   1739 ]
-> 1741 return compile_pymc(inputs, fg.outputs, givens=givens, on_unused_input="ignore", **kwargs)

File ~\anaconda3\envs\test\lib\site-packages\pymc\aesaraf.py:1034, in compile_pymc(inputs, outputs, random_seed, mode, **kwargs)
   1032 opt_qry = mode.provided_optimizer.including("random_make_inplace", check_parameter_opt)
   1033 mode = Mode(linker=mode.linker, optimizer=opt_qry)
-> 1034 aesara_function = aesara.function(
   1035     inputs,
   1036     outputs,
   1037     updates={**rng_updates, **kwargs.pop("updates", {})},
   1038     mode=mode,
   1039     **kwargs,
   1040 )
   1041 return aesara_function

File ~\anaconda3\envs\test\lib\site-packages\aesara\compile\function\__init__.py:317, in function(inputs, outputs, mode, updates, givens, no_default_updates, accept_inplace, name, rebuild_strict, allow_input_downcast, profile, on_unused_input)
    311     fn = orig_function(
    312         inputs, outputs, mode=mode, accept_inplace=accept_inplace, name=name
    313     )
    314 else:
    315     # note: pfunc will also call orig_function -- orig_function is
    316     #      a choke point that all compilation must pass through
--> 317     fn = pfunc(
    318         params=inputs,
    319         outputs=outputs,
    320         mode=mode,
    321         updates=updates,
    322         givens=givens,
    323         no_default_updates=no_default_updates,
    324         accept_inplace=accept_inplace,
    325         name=name,
    326         rebuild_strict=rebuild_strict,
    327         allow_input_downcast=allow_input_downcast,
    328         on_unused_input=on_unused_input,
    329         profile=profile,
    330         output_keys=output_keys,
    331     )
    332 return fn

File ~\anaconda3\envs\test\lib\site-packages\aesara\compile\function\pfunc.py:374, in pfunc(params, outputs, mode, updates, givens, no_default_updates, accept_inplace, name, rebuild_strict, allow_input_downcast, profile, on_unused_input, output_keys, fgraph)
    360     profile = ProfileStats(message=profile)
    362 inputs, cloned_outputs = construct_pfunc_ins_and_outs(
    363     params,
    364     outputs,
   (...)
    371     fgraph=fgraph,
    372 )
--> 374 return orig_function(
    375     inputs,
    376     cloned_outputs,
    377     mode,
    378     accept_inplace=accept_inplace,
    379     name=name,
    380     profile=profile,
    381     on_unused_input=on_unused_input,
    382     output_keys=output_keys,
    383     fgraph=fgraph,
    384 )

File ~\anaconda3\envs\test\lib\site-packages\aesara\compile\function\types.py:1763, in orig_function(inputs, outputs, mode, accept_inplace, name, profile, on_unused_input, output_keys, fgraph)
   1751     m = Maker(
   1752         inputs,
   1753         outputs,
   (...)
   1760         fgraph=fgraph,
   1761     )
   1762     with config.change_flags(compute_test_value="off"):
-> 1763         fn = m.create(defaults)
   1764 finally:
   1765     t2 = time.time()

File ~\anaconda3\envs\test\lib\site-packages\aesara\compile\function\types.py:1656, in FunctionMaker.create(self, input_storage, trustme, storage_map)
   1653 start_import_time = aesara.link.c.cmodule.import_time
   1655 with config.change_flags(traceback__limit=config.traceback__compile_limit):
-> 1656     _fn, _i, _o = self.linker.make_thunk(
   1657         input_storage=input_storage_lists, storage_map=storage_map
   1658     )
   1660 end_linker = time.time()
   1662 linker_time = end_linker - start_linker

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\basic.py:254, in LocalLinker.make_thunk(self, input_storage, output_storage, storage_map, **kwargs)
    247 def make_thunk(
    248     self,
    249     input_storage: Optional["InputStorageType"] = None,
   (...)
    252     **kwargs,
    253 ) -> Tuple["BasicThunkType", "InputStorageType", "OutputStorageType"]:
--> 254     return self.make_all(
    255         input_storage=input_storage,
    256         output_storage=output_storage,
    257         storage_map=storage_map,
    258     )[:3]

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\vm.py:1255, in VMLinker.make_all(self, profiler, input_storage, output_storage, storage_map)
   1253             thunks[-1].lazy = False
   1254     except Exception:
-> 1255         raise_with_op(fgraph, node)
   1257 t1 = time.time()
   1259 if self.profile:

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\utils.py:534, in raise_with_op(fgraph, node, thunk, exc_info, storage_map)
    529     warnings.warn(
    530         f"{exc_type} error does not allow us to add an extra error message"
    531     )
    532     # Some exception need extra parameter in inputs. So forget the
    533     # extra long error message in that case.
--> 534 raise exc_value.with_traceback(exc_trace)

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\vm.py:1246, in VMLinker.make_all(self, profiler, input_storage, output_storage, storage_map)
   1241 thunk_start = time.time()
   1242 # no-recycling is done at each VM.__call__ So there is
   1243 # no need to cause duplicate c code by passing
   1244 # no_recycling here.
   1245 thunks.append(
-> 1246     node.op.make_thunk(node, storage_map, compute_map, [], impl=impl)
   1247 )
   1248 linker_make_thunk_time[node] = time.time() - thunk_start
   1249 if not hasattr(thunks[-1], "lazy"):
   1250     # We don't want all ops maker to think about lazy Ops.
   1251     # So if they didn't specify that its lazy or not, it isn't.
   1252     # If this member isn't present, it will crash later.

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\c\op.py:131, in COp.make_thunk(self, node, storage_map, compute_map, no_recycling, impl)
    127 self.prepare_node(
    128     node, storage_map=storage_map, compute_map=compute_map, impl="c"
    129 )
    130 try:
--> 131     return self.make_c_thunk(node, storage_map, compute_map, no_recycling)
    132 except (NotImplementedError, MethodNotDefined):
    133     # We requested the c code, so don't catch the error.
    134     if impl == "c":

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\c\op.py:96, in COp.make_c_thunk(self, node, storage_map, compute_map, no_recycling)
     94         print(f"Disabling C code for {self} due to unsupported float16")
     95         raise NotImplementedError("float16")
---> 96 outputs = cl.make_thunk(
     97     input_storage=node_input_storage, output_storage=node_output_storage
     98 )
     99 thunk, node_input_filters, node_output_filters = outputs
    101 @is_cthunk_wrapper_type
    102 def rval():

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\c\basic.py:1198, in CLinker.make_thunk(self, input_storage, output_storage, storage_map)
   1170 """
   1171 Compiles this linker's fgraph and returns a function to perform the
   1172 computations, as well as lists of storage cells for both the inputs
   (...)
   1195   first_output = ostor[0].data
   1196 """
   1197 init_tasks, tasks = self.get_init_tasks()
-> 1198 cthunk, module, in_storage, out_storage, error_storage = self.__compile__(
   1199     input_storage, output_storage, storage_map
   1200 )
   1202 res = _CThunk(cthunk, init_tasks, tasks, error_storage, module)
   1203 res.nodes = self.node_order

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\c\basic.py:1133, in CLinker.__compile__(self, input_storage, output_storage, storage_map)
   1131 input_storage = tuple(input_storage)
   1132 output_storage = tuple(output_storage)
-> 1133 thunk, module = self.cthunk_factory(
   1134     error_storage,
   1135     input_storage,
   1136     output_storage,
   1137     storage_map,
   1138 )
   1139 return (
   1140     thunk,
   1141     module,
   (...)
   1150     error_storage,
   1151 )

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\c\basic.py:1629, in CLinker.cthunk_factory(self, error_storage, in_storage, out_storage, storage_map)
   1627     for node in self.node_order:
   1628         node.op.prepare_node(node, storage_map, None, "c")
-> 1629     module = get_module_cache().module_from_key(key=key, lnk=self)
   1631 vars = self.inputs + self.outputs + self.orphans
   1632 # List of indices that should be ignored when passing the arguments
   1633 # (basically, everything that the previous call to uniq eliminated)

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\c\cmodule.py:1228, in ModuleCache.module_from_key(self, key, lnk)
   1226 try:
   1227     location = dlimport_workdir(self.dirname)
-> 1228     module = lnk.compile_cmodule(location)
   1229     name = module.__file__
   1230     assert name.startswith(location)

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\c\basic.py:1538, in CLinker.compile_cmodule(self, location)
   1536 try:
   1537     _logger.debug(f"LOCATION {location}")
-> 1538     module = c_compiler.compile_str(
   1539         module_name=mod.code_hash,
   1540         src_code=src_code,
   1541         location=location,
   1542         include_dirs=self.header_dirs(),
   1543         lib_dirs=self.lib_dirs(),
   1544         libs=libs,
   1545         preargs=preargs,
   1546     )
   1547 except Exception as e:
   1548     e.args += (str(self.fgraph),)

File ~\anaconda3\envs\test\lib\site-packages\aesara\link\c\cmodule.py:2639, in GCC_compiler.compile_str(module_name, src_code, location, include_dirs, lib_dirs, libs, preargs, py_module, hide_symbols)
   2631                 print(
   2632                     "Check if package python-dev or python-devel is installed."
   2633                 )
   2635     # We replace '\n' by '. ' in the error message because when Python
   2636     # prints the exception, having '\n' in the text makes it more
   2637     # difficult to read.
   2638     # compile_stderr = compile_stderr.replace("\n", ". ")
-> 2639     raise CompileError(
   2640         f"Compilation failed (return status={status}):\n{' '.join(cmd)}\n{compile_stderr}"
   2641     )
   2642 elif config.cmodule__compilation_warning and compile_stderr:
   2643     # Print errors just below the command line.
   2644     print(compile_stderr)

CompileError: Compilation failed (return status=1):
"C:\Users\andri\anaconda3\envs\test_pymc3\Library\mingw-w64\bin\g++.exe" -shared -g -O3 -fno-math-errno -Wno-unused-label -Wno-unused-variable -Wno-write-strings -Wno-c++11-narrowing -fno-exceptions -fno-unwind-tables -fno-asynchronous-unwind-tables -march=broadwell -mmmx -mno-3dnow -msse -msse2 -msse3 -mssse3 -mno-sse4a -mcx16 -msahf -mmovbe -maes -mno-sha -mpclmul -mpopcnt -mabm -mno-lwp -mfma -mno-fma4 -mno-xop -mbmi -mbmi2 -mno-tbm -mavx -mavx2 -msse4.2 -msse4.1 -mlzcnt -mno-rtm -mno-hle -mrdrnd -mf16c -mfsgsbase -mrdseed -mprfchw -madx -mfxsr -mxsave -mxsaveopt -mno-avx512f -mno-avx512er -mno-avx512cd -mno-avx512pf -mno-prefetchwt1 -mclflushopt -mxsavec -mxsaves -mno-avx512dq -mno-avx512bw -mno-avx512vl -mno-avx512ifma -mno-avx512vbmi -mno-clwb -mno-pcommit -mno-mwaitx --param l1-cache-size=32 --param l1-cache-line-size=64 --param l2-cache-size=12288 -mtune=generic -DNPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION -m64 -DMS_WIN64 -I"C:\Users\andri\anaconda3\envs\test\lib\site-packages\numpy\core\include" -I"C:\Users\andri\anaconda3\envs\test\include" -I"C:\Users\andri\anaconda3\envs\test\lib\site-packages\aesara\link\c\c_code" -L"C:\Users\andri\anaconda3\envs\test\libs" -L"C:\Users\andri\anaconda3\envs\test" -o "C:\Users\andri\AppData\Local\Aesara\compiledir_Windows-10-10.0.22000-SP0-Intel64_Family_6_Model_158_Stepping_10_GenuineIntel-3.8.13-64\tmp11_4bwf_\m47260be5189a2297f95a5722358fab6ed80907bc9b3b7adb3279eddf44b57064.pyd" "C:\Users\andri\AppData\Local\Aesara\compiledir_Windows-10-10.0.22000-SP0-Intel64_Family_6_Model_158_Stepping_10_GenuineIntel-3.8.13-64\tmp11_4bwf_\mod.cpp" -lblas "C:\Users\andri\anaconda3\envs\test\python38.dll"
C:/Users/andri/anaconda3/envs/test_pymc3/Library/mingw-w64/bin/../lib/gcc/x86_64-w64-mingw32/5.3.0/../../../../x86_64-w64-mingw32/bin/ld.exe: cannot find -lblas
collect2.exe: error: ld returned 1 exit status

Apply node that caused the error: Dot22Scalar(Elemwise{true_div,no_inplace}.0, InplaceDimShuffle{1,0}.0, TensorConstant{-2.0})
Toposort index: 19
Inputs types: [TensorType(float64, (None, 1)), TensorType(float64, (1, None)), TensorType(float64, ())]

HINT: Use a linker other than the C linker to print the inputs' shapes and strides.
HINT: Re-running with most Aesara optimizations disabled could provide a back-trace showing when this node was created. This can be done by setting the Aesara flag 'optimizer=fast_compile'. If that does not work, Aesara optimizations can be disabled with 'optimizer=None'.
HINT: Use the Aesara flag `exception_verbosity=high` for a debug print-out and storage map footprint of this Apply node.

I believe that this is triggered by a shape mismatch but I cannot figure out where it comes from. I remember having this problem in the past but it had something to do specifically with shared data variables. It must also be connected to the new way that v4 handles shapes because I can run this exact code using v3. However, I should mention that for v3 I need to restrict the cores to 1 or it throws the “usual” error cannot pickle a fortran object. So I guess that the two things must be somehow related.

Does anyone have an idea of what could be going wrong?
Thanks!

2

Your installation looks broken, specifically blas is missing. I would advise to install pymc via conda-forge as per the official instructions if that’s a possibility: Installation — PyMC 4.1.4 documentation

3

Hi @ricardoV94, thanks for your response (and also for the very nice Aesara, Aeppl,… intro you ran a few weeks ago btw, it was very helpful:). I’m gonna take advantage of this chance sorry, can I ask if there is a way to “create” a new distribution in the way you discussed in the video for the case of log Betas but for a truncated distribution, e.g. truncated StudentT?).

I did install pymc following the instructions (and re-created the environment a couple of times to make sure) and it works just fine with other models, which is the main reason of my confusion. This -blas is missing error has come up before and I could never pin it down (once I think it went away with removing the assignment of the data to shared Data variables) but it does not seem to be connected to the installation unless it could be that only certain model structures could be affected? In this case the error appears only if I add the f_time = gp.prior("f_time", X=months[:, None]) line.

1 Like
4

Is there a possibility to provide a fully reproducible example to check if it happens on my machine as well? If you can narrow down to the smallest model that fails with that error?

@bwengals do you see any problem with the GP specification?

5

I plan to implement truncated distributions pretty soon! It’s not yet available but you can hack something with the logcdf functions already if you need

6

Nope it looks right to me. I agree though it looks like an environment issue, ld.exe: cannot find -lblas. I think you aren’t getting this error when you comment out the GP because it’s the only part of your model that uses blas.

7

Ah this is great thanks!
Do you know of any example where this hack is done just as a reference?

8

Thank you both, ok I will try to rebuild the environment once more to see if I can sort this out. Otherwise I will put together something I can share to see if you have the same issue. Thanks again for now!

9

Hi all,

I rebuilt the environment (using exclusively conda create -c conda-forge -n pymc_env "pymc>=4" but still I have the same problem.
The only thing I can think of is that this is the result of an incorrect path search on Windows machines (see the slashes direction in the line cannot find -lblas). This seems weird but I don’t know what else could be going wrong.
I put together a minimal example that maybe you could try to run (possibly on a Windows machine)? Incidentally, while I was doing this I remembered what was the problem I had with using shared data variables. When I set mutable=True, the same issue with cannot find -lblas emerges (sometimes only, I have not been able to pin down why this does not occur every time but it has something to do with the order in which I comment out other lines), which points to the fact that indeed blas jumps in for some things and as long as I don’t need it everything else runs fine.

So, in short, you should be able to run the example below just fine but, when you un-comment the f_time = ... and the + f_time lines, the error I posted above should pop up.
Thanks again for your help!

import numpy as np
import pymc as pm

seed = sum(map(ord, "Test gaussian processes with pymc v4"))
rng = np.random.default_rng(seed=seed)

time = np.linspace(1, 100, 100)
covariate = rng.normal(0, 1, 100)
target_obs = (
    pm.GaussianRandomWalk.dist(0, 0.2, steps=99, init_dist=pm.Normal.dist(0, 0.1)).eval()
    + (0.5 * covariate)
)

coords = {"time": time}

with pm.Model(coords=coords):
    
    baseline = pm.Normal("baseline", mu=target_obs.mean(), sigma=0.11)
    covariate_param = pm.Normal("no_homeless", mu=0.2, sigma=0.1)
    
    # time trends
    amplitude_trend = pm.HalfNormal("amplitude_trend", sigma=3)
    ls_trend = pm.Gamma("ls_trend", alpha=5, beta=2)
    cov_trend = amplitude_trend ** 2 * pm.gp.cov.Matern52(1, ls_trend)
    
    gp = pm.gp.Latent(cov_func=cov_trend)
    # f_time = gp.prior("f_time", X=time[:, None])
    
    # data
    covariate_data = pm.Data("covariate_data", covariate, dims="time", mutable=True)
    
    target_mu = pm.Deterministic(
        "target_mu",
        var = (
            baseline 
            # + f_time 
            + covariate_param * covariate_data 
        ),
        dims="time"
    )
    target_sigma = pm.HalfNormal("dem_sigma", sigma=0.2)
    
    # Likelihood
    target = pm.Normal(
        "target", 
        mu=target_mu, 
        sigma=target_sigma, 
        observed=target_obs, 
        dims="time"
    )
    
    priors = pm.sample_prior_predictive()
1 Like
10

Hi @ricardoV94 and @bwengals!
Apologies for bugging you, I was just wondering if you had a chance to run the piece of code above? I’m curious to see if there is something I can do to get this running on my Windows machine.

Thanks!

1 Like
11

I’ll give it a try when I am in front of a machine.

In the meantime you can try it in Google Colab. Just run !pip install ”pymc>4” in the first cell

1 Like
12

Will do, thanks a lot Ricardo!

13

Hi, I’m following this discussion with interest, because I encountered the same error message and after googling ended up here. I was able to reproduce the minimum code example from Paolo with the same error message.

Thanks,

14

Apologies for the delay. I have finally opened that PR Implement Truncated distributions by ricardoV94 · Pull Request #6113 · pymc-devs/pymc · GitHub

Anyway you can implement a truncated distribution by adding an interval transform so it doesn’t sample out of the truncation bounds and a Potential with the needed normalization constant. If it’s only for the likelihood you don’t need the transform.

There are some examples here: Google Colab which came from Help with pm.Potential - #12 by sspickle

15

Hey Ricardo!
Thanks a lot for this update, this will be super helpful!

On another note, unfortunately, I’m still having the same issues with lblas with v4 on my local machine. I haven’t managed to figure out what’s causing them but it sounds like it’s not only me with this problem.

16

If the discourse forum was of no help, feel free to open an issue on Github, perhaps the devs there might be more helpful.

17

Will do, thanks again Ricardo!