There has been much interest in analyzing genome-scale DNA sequence data to infer population histories, but inference methods developed hitherto are limited in model complexity and computational scalability. Here we present an efficient, flexible statistical method, diCal2, that can utilize wholegenome sequence data from multiple populations to infer complex demographic models involving population size changes, population splits, admixture, and migration. Applying our method to data from Australian, East Asian, European, and Papuan populations, we find that the population ancestral to Australians and Papuans started separating from East Asians and Europeans about 100,000 years ago, and that the separation of East Asians and Europeans started about 50,000 years ago, with pervasive gene flow between all pairs of populations.1. An arbitrary number of populations specified by the user.2. An arbitrary pattern of population splits and mergers. 3. More general population size changes (e.g., piecewise-exponential). 4. Arbitrary migration patterns with time-varying continuous migration rates or pulse admixture events. 5. An arbitrary poly-allelic mutation model at each site (including di-or tetra-allelic).