New sample calculation with proof

[rajeee]

Pull Request Description

(Proof refers to the assertions in the code that compares new samples to old samples)

Adds python version of run_sampling.rb
Simplifies sample calculation (for modest speed boost) in run_sampling.rb.

Checklist

Not all may apply:

• Tests (and test files) have been updated
• Documentation has been updated
• Changelog has been updated
• openstudio tasks.rb update_measures has been run
• No unexpected regression test changes on CI (checked comparison artifacts)

[rajeee]

New version.
Once the new version is accepted to be okay, we can remove the old version. Keeping both for now to make it possible to compare.

[rajeee]

For most cases, the new version is equal to old version; When they are not, making sure that the difference can be justified.

[shorowit]

I sampled project_national w/ 350k datapoints on Windows Subsystem for Linux.

• run_sampling.rb (v1): 4m6s
• run_sampling.rb (v2): 3m32s
• run_sampling.py: 16m28s

All three produced CSV files of roughly the same size, so that's a good sign.

Logs:
run_sampling_rb_v1.txt
run_sampling_rb_v2.txt
run_sampling_py.txt

[rajeee]

@shorowit Thanks!
I am really surprised with the runtimes though. In my benchmarking, python version was running much faster! I ran it directly on windows though.

[rajeee]

Here is what I got on Mac and on Eagle.
My M1 Mac
Python: 33.4s
Ruby: 1m53s

Eagle Login Node:
Python: 1m29s
Ruby: 3m18s

run_sampling.py.m1mac.txt
run_sampling.py.eagle_login_node.txt
run_sampling.rb.m1mac.txt
run_sampling.rb.eagle_login_node.txt

WSL seems to be doing something really bad for python!

[shorowit]

Here are my results on Windows.

• run_sampling.rb (v1): 8m12s
• run_sampling.rb (v2): 7m36s
• run_sampling.py: 1m1s

Very impressive! Definitely helps to make use of multiple cpu cores (my machine has 16). Not sure what's going on w/ my WSL results above.

[aspeake]

@rajeee @shorowit thanks guys. Sounds like this could be really helpful to @vtnate and myself. We will run some checks and most likely use it for our custom sampler in the buildstock-projections library. I'll keep you updated on the progress of the down-selection feature.

[rajeee]

There is a certain limitation of current sampling approach.

Consider this simplified case of just two TSVs.
Bedrooms.TSV

Option=1	Option=2	Option=3	Option=4	Option=5	sampling_probability
0.2	0.2	0.2	0.2	0.2	1

Ceiling Fan.TSV

Dependency=Bedrooms	Option=None	Option=Standard Efficiency	Option=Premium Efficiency	sampling_probability
1	0.4	0.4	0.2	0.2
2	0.4	0.4	0.2	0.2
3	0.4	0.4	0.2	0.2
4	0.4	0.4	0.2	0.2
5	0.4	0.4	0.2	0.2

For n=10, one would expect two samples each for bedrooms = 1 through 5.
One would also expect, 4 samples for Ceiling Fan = None, 4 samples for Ceiling Fan = Standard Efficiency and 2 samples for Ceiling Fan=Premium Efficiency.

However, run_sampling.rb on these TSVs would generate these samples on buildstock.csv:

Building	Bedrooms	Ceiling Fan
4	1	None
9	1	Standard Efficiency
2	2	None
5	2	Standard Efficiency
3	3	None
8	3	Standard Efficiency
1	4	None
10	4	Standard Efficiency
6	5	Standard Efficiency
7	5	None

There are 5 samples each for None and Standard Efficiency and no samples for Premium Efficiency.

This limitation is easily explained by the way current sampling works. It first samples Bedrooms that has no dependency, so it assigns two samples for each number of bedrooms. Then for each of these two samples, it tries to distribute the Fan Efficiency. Since there are 3 fan efficiencies, but only two samples, it picks the two option with the highest probability. While this ensures that the distribution follows the TSV specification as close as possible for this particular bedroom, on doing so for all bedrooms, we miss out on the fact that the probability for Option=Premium, small as it may be, should add up and we should eventually get some samples for it.

Although this issue is demonstrated for a case where the number of samples is very small, this phenomena can occur even for the case of very large n when sampling within a narrow slice created by a large number of dependencies and their options.

I am thinking of a slight modification to the sampling approach where we accumulate the probability of discarded options, and use that to ensure that low probability options also get some samples.

[shorowit]

In my opinion, this is a feature, not a limitation/bug. Indeed the sampling algorithm was intentionally developed this way.

In your example, if you did end up end up with 2 Premium Efficiency ceiling fans, they would be randomly assigned to buildings. If you were to then look at ceiling fan energy use as a function of number of bedrooms, the trend would have random spikes based on this random assignment and a user would likely be confused. Or, a user might grab the two 2-bedroom building models and use them for an analysis assuming that these models reflect the highest probability building characteristics, but they would not.

The bottom line is that there is no good sampling result when you have too few samples for a given distribution. And it's important that ResStock maintain enough samples given its distributions. In the absence of having enough samples, your proposed approach may be what some users want and the current approach may be what other users want. The current algorithm reflects our consensus from many years ago as to which use cases are more important.

@afontani may have thoughts here, I think that he and I had some discussions about all of this many years ago.

[rajeee]

@shorowit
Thanks for the commentary.

I agree it's trade-off between two things - making global distribution more accurate vs making local distribution more accurate.

My feelings was that we have more confidence on global distribution (what % of homes nationally have premium fan) than on local distribution (of homes that have 2 bedrooms, what % have premium fan), so it might be more suitable to prioritize that. But, yeah, my presumption maybe wrong.

Also, I might be raising sort of a false flag. I did notice some parameters had slightly (very tiny in fact; couple of samples in a 100K total samples), different global distribution than what one would expect based on sampling_probability calculation. I was investigating what could have caused that. I am working on a test to validate the generated buildstock.csv against the TSVs and was investigating this as a part of that. I will see if we even have any substantial discrepancy for any parameters; maybe we don't.

Anyway, this was an interesting exercise on understanding the behavior of sampling.

[shorowit]

@rajeee I completely agree with your assessment. And to be clear, I'm not arguing that the current behavior is correct just because it's always been this way. I look forward to the results of your investigations and any improvements you can make! Developing a test around all of this is a definite good first step.

[afontani]

@rajeee: Thanks for thinking through this. I also look forward to any results or comparisons at larger sample sizes to see if there are any issues.

I agree with @shorowit about this not actually being a bug but a feature. I am not positive, but I think this behavior helps the algorithm converge to the distributions at $O(n^{-1})$ vs $O(n^{-1/2})$. I think the point of quota-based sampling is to sample the center of the distribution so that the estimate of the mean quantity of interest has a reasonable estimate at smaller sample sizes. Ultimately, I would want any modification to the algorithm to preserve the $O(n^{-1})$ convergence rate.

I think we will have issues with small sample sizes regardless of the algorithm. In this case, I would prefer the algorithm as it is currently written. But for larger sample sizes, like for a national run, 500k samples, I am open to modifying if there are use cases that properly motivate the change. I do worry about how people slice our data from a large run and if the user is getting a more representative set of samples or a set of edge cases.

[rajeee]

@shorowit @afontani Thanks for the comments. As per your feedback, I decided to keep the current approach since it has better accuracy in regions where we have low samples.

I continued the development of the python sampler here: NREL/buildstockbatch#308
In addition to sampling, it now also has a verifier that can verify buildstock.csv for correctness. It also calculates "sampling_errors" (for lack of a better word) for each TSV. You can learn more in that PR link.

I verified the samples generated by all three versions (ruby_v1, ruby_v2 and python) to be correct.

I did 5 sampling runs at n=350K between ruby_v1, ruby_v2 and python sampler using different random seed each time and compared the sampling_error (errors for both dependency groups, and errors for options). There is no qualitative difference between them. This is for 350K samples. So, we have up-to 10% option distribution error and 18% group distribution error for some of the large TSVs.

ZOOMED

[rajeee]

Error distribution convergence:

Convergence of the average errors:

Log-log plot asymptomatically converges to slope of -1.

[shorowit]

Moved to NREL/buildstockbatch#308