This is a code used in a computer simulations discussed in paper with woking title **Assortative pairing alone leads to structured biota in organisms with cultural transmission**. Their results show, that for some forms of positional inheritance (inspired by old biometric theories), assortment alone can be responsible for the formation of groups with a large between-group and small within-group variance. No ecological specialization or deliberate selective imitation is necessary, the proportionality of "parent" and "offspring" variance does the heavy lifting here. Equivalents of "Urban subcultures" emerge with no additional assumptions. The preprint can be found at: https://doi.org/10.48550/arXiv.2202.13685 Folder **Single run** contains easily intelligable R code that allows the reader to **simulate one population and develop it with specified set of constants for a specified number of generations** (the functions to replice the simulations are stored in *simulation_combined_nonscaled_hdbscan.R* ). Functions that do the clustering, reduce the dimensionality of the simulated space to 3D and rotate it to mimimize mismatch between adjecent panels for nice visulaizations are included. It is possible to draw multi-panel tiffs ( *visualize_4panels_hdbscan.R* ) or pdf animations ( *visualize_animation_hdbscan.R* ). One example of each is included in the folder. It is best to start from *single_run_parade.Rproj* and run script *exemplar_runs_r.R*. This script recreates all demonstrations of single runs used in the article and accomapnying supplement. It sources all necessary scripts. Scripts *visual_functions_all.R* and *visualize_niceplot.R* store bits used in the two higher-order functions in *visualize_4panels_hdbscan.R* and *visualize_4panels_hdbscan.R*. (We used R studio to run all the scripts.) **To run multiple simulations on a grid of parameter values** as quickly as possible, use the scripts in folder **Parallel runs with Numpy**. It is best to first open the file *nw_example.py* first and check if everything works properly in Python (all .py scripts are writen in Python 3.7, we used Spyder to run the scripts). It is possible to test whether the infrastructure works directly in Python by un-commenting lines 52-67 (equivalent to running a single simulation in **Single run**). Having Anaconda installed is highly recommended (It installs library Numpy automatically). The only other package that needs to be forged separately is *dbscan*. Script *nw_example.py* sources and imports functions from the other two pzthon scripts that can be explored for details. These scripts can be loaded and functions can be called directly from R that is also used to parallelize the task on multiple cores (using *parallel* R package). It is, again, best to start from the *run_parallel.Rproj* project file. Then, the series of simulations can be managed from *grid_generate_parallel.R*. At the beginning of the script, *generate_function_core.R* is loaded onto each core - this includes loading R package *reticulate* that allows to call Python code from R. Make sure this package is installed on your computer. Because generating a lot of simulations can take long on a regular desktop, an exemplar dataset *results.all.txt* is included. Each simulation's results are stored in a single line of the dataset. This dataset can be used to draw all summary images included in the manuscript (such as *Figure_contour_01_clust.tif* and *Figure_image_01_clust.tif* included here as examples). To summarize and visualize the output, code in *grid_visualize_3by3_with_crosses.R* is used (this repository also presents its simple equivalent that does not draw crosses per exemplar simulation runs). It contains, among other things, function *doFloorC* that can be employed to visualize any composite image resulting from simulations on a multi-D grid of parameters. Each "floor" is 2D and theese floors can be easily arranged in higer order configurations - like the 3×3 panel arrangement - using the higher level function *drawFloors*. These functions can draw both "contour" and "image" plots. Since we had enough simulation runs across all parameter values, image plots that do not interpolate any information from the results were used. For quick approximation with lower number of simulation runs, contour plots may come useful. The script demonstrates how to draw both. Each folder **Parallel runs with Numpy** and **Single run** contains two folders with Alternatives - slight modifications of the algorithm discussed in the supplement of the paper. The set of bell curves that illustrate the difference between Galton-Pearson and PVDI inheritance model, can be produced with the script in *Figure1.R*. Numerical solutions of some boundary values of nu in PVDI model can be found in *Prob_between.R* script.