Binary population synthesis (BPS) is an essential tool for extracting information about massive binary evolution from gravitational-wave (GW) detections of compact object mergers. It has been successfully used to constrain the most likely permutations of evolution models among hundreds of alternatives, while initial condition models, in contrast, have not yet received the same level of attention. Here, we introduce BOSSA, a detailed initial sampling code including a set of 192 initial condition permutations for BPS that capture both ”invariant” and ”varying” models, the latter accounting for a possible metallicity- and star formation rate (SFR)-dependence of the initial mass function (IMF); as well as correlations between the initial primary mass, orbital period, mass ratio and eccentricity of binaries. We include 24 metallicity-specific cosmic star formation history (cSFH) models and propose two alternate models for the mass-dependent binary fraction. We build a detailed pipeline for time-evolving BPS, such that each binary has well-defined initial conditions, and we are able to distinguish the contributions from populations of different ages. We discuss the meaning of the IMF for binaries and introduce a refined initial sampling procedure for component masses. We also discuss the treatment of higher-order multiple systems when normalizing a binary sample. In particular, we argue for how a consistent interpretation of the IMF implies that this is not the distribution from which any set of component masses should be independently drawn, and show how the individual IMF of primaries and companions is expected to deviate from the full IMF.