Algorithm to sample high-dimensional parameter space of expensive cost function I am looking for recommendations for algorithms (apparently this is the right SO site for that) that efficiently scan the high-dimensional parameter space of a cost function that is very expensive to evaluate.

By expensive to evaluate I mean that one can only do O(1) evaluations per second on a modern computer.
By efficiently I mean that the algorithm should learn to avoid regions where there are no solutions: the evaluation of the cost function may fail for certain sets of parameters.
Still under the aspect of efficiency, the algorithm should yield a set of samples that represent a (ideally unbiased) coverage of all solutions in the whole parameter space within a reasonable number of samples.
The parameters are continuous numbers within predefined bounds.
The cost function is &quot;mostly well-behaved&quot;.

I am not so much after importance sampling, so I guess I am looking for an improved random sampling that is adaptive in the sense that it learns which parts of parameter space to avoid.
Background information:
I am looking into a specific physics model where from ~20 input parameters a set of outputs (masses and decay probabilities) are computed. This is the bottleneck CPU-wise. My goal is to identify allowed ranges for the input parameters for which the outputs fulfil certain constraints that are implemented in terms of a cost function. Not all possible combinations of input parameters give meaningful physics, therefore a large (but acceptably so) fraction of calculations will fail.
At the moment I am successfully using parallel simulated annealing to find solutions (input parameters) that come very close to the known minimum of the cost function. So finding a (nearly) optimal solution is &quot;easy&quot;.
What I would now like to improve is the coverage: sampling from the whole input parameter space (e.g. below some threshold of the cost function).
This would allow projections onto 1-D or 2-D subspaces of the input-parameter space to be done to see where valid models lie.
In this example

the yellow bin holds most of the samples because the cost function depends on this particular output such that it ought to be close to 125. The interesting part though is the dependence of the output on the input. One can guess from the plot that lower values of the input would also give the desired output, and below ~ 1000 we see some structure, where the input apparently is too low to achieve the desired output. However, the statistics is very low. Therefore, I am looking for another algorithm which yields a sampling that covers input-parameter space more evenly (subject to the constraint that invalid input parameter sets should be avoided).
I don't know how to sample from the space of nearly-optimal solutions, but I'll share one approach you could consider.
It appears that you want to estimate the density of good solutions, at different points in the input space.  It appears you are gridding up the parameter space and estimating the density in each cell of the grid.  Here is a possible alternate method to do that.
Suppose we have a cell $C \subseteq \mathbb{R}^{20}$ in the parameter space, which might take the form $C = [\ell_1,u_1] \times \cdots \times [\ell_{20},u_{20}]$, i.e., $\ell_i$ is the lower bound on the $i$th coordinate and $u_i$ the upper bound.   Let $V(C)$ be the volume of cell $C$.  We'll estimate the density of good solutions in $C$.  Pick a volume $v$ that is a fraction of $V(C)$, let $\alpha=v^{1/20}$, and then repeat the following 100 times:

Pick a random point $w \in C$.  Form the smaller cell $C_x = [w_1-\alpha/2,w_1+\alpha/2] \times \cdots \times [w_{20}-\alpha/2,w_{20}+\alpha/2]$.  Use your optimization method to search whether there is a good solution in the cell $C_x$.  (For many optimization methods, it is easy to enforce constraints that set an upper and lower bound on each coordinate in the parameter space.)

Now based on the number of times that you were able to find a good solution, you can form an estimate of the density of good solutions in $C$.  If this count is close to zero, double $v$ and repeat.  If this count is close to 100, halve $v$ and repeat.  Repeat until you get a count that is not too close to 0 or 100, and the result should give you a pretty good estimate of the density in $C$.
You can now iterate this, once per cell $C$, and use that to produce the kind of plot you are hoping for.
Here I've adopted standard methods for #SAT, i.e., for estimating the number of solutions to a SAT instance, to your specific problem.