personal/time-to-botec
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

don't download whole git dir for squiggle_c

Browse Source
This commit is contained in:
NunoSempere 2023-11-23 14:39:25 +00:00
parent 05b1b99561
commit 06e78cb079
65 changed files with 13 additions and 1002187 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
squiggle.c/makefile
View File

@ -4,6 +4,18 @@ build:
gcc -O3 samples.c ./squiggle_c/squiggle.c ./squiggle_c/squiggle_more.c -lm -fopenmp -o $(OUTPUT) gcc -O3 samples.c ./squiggle_c/squiggle.c ./squiggle_c/squiggle_more.c -lm -fopenmp -o $(OUTPUT)
install: install:
rm -r squiggle_c
wget https://git.nunosempere.com/personal/squiggle.c/raw/branch/master/squiggle.c
wget https://git.nunosempere.com/personal/squiggle.c/raw/branch/master/squiggle.h
wget https://git.nunosempere.com/personal/squiggle.c/raw/branch/master/squiggle_more.c
wget https://git.nunosempere.com/personal/squiggle.c/raw/branch/master/squiggle_more.h
mkdir temp
mv squiggle* temp
mv temp squiggle_c
wget https://git.nunosempere.com/personal/squiggle.c/raw/branch/master/examples/more/12_time_to_botec_parallel/example.c -O samples.c
sed -i 's|../../..|squiggle_c|' samples.c
install-git:
rm -r squiggle_c rm -r squiggle_c
git clone https://git.nunosempere.com/personal/squiggle.c git clone https://git.nunosempere.com/personal/squiggle.c
mv squiggle.c squiggle_c mv squiggle.c squiggle_c

7
squiggle.c/squiggle_c/LICENSE.txt
View File

@ -1,7 +0,0 @@
Copyright 2023 and ss., Nuño Sempere
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

385
squiggle.c/squiggle_c/README.md
View File

@ -1,385 +0,0 @@
# squiggle.c
squiggle.c is a [grug-brained](https://grugbrain.dev/) self-contained C99 library that provides functions for simple Monte Carlo estimation, based on [Squiggle](https://www.squiggle-language.com/).
## Why C?
- Because it is fast
- Because I enjoy it
- Because C is honest
- Because it will last long
- Because it can fit in my head
- Because if you can implement something in C, you can implement it anywhere else
- Because it can be made faster if need be
- e.g., with a multi-threading library like OpenMP,
- o by implementing faster but more complex algorithms
- or more simply, by inlining the sampling functions (adding an `inline` directive before their function declaration)
- **Because there are few abstractions between it and machine code** (C => assembly => machine code with gcc, or C => machine code, with tcc), leading to fewer errors beyond the programmer's control.
## Getting started
You can follow some example usage in the examples/ folder
1. In the [1st example](examples/01_one_sample/example.c), we define a small model, and draw one sample from it
2. In the [2nd example](examples/02_many_samples/example.c), we define a small model, and return many samples
3. In the [3rd example](examples/03_gcc_nested_function/example.c), we use a gcc extension—nested functions—to rewrite the code from point 2. in a more linear way.
4. In the [4th example](examples/04_sample_from_cdf_simple/example.c), we define some simple cdfs, and we draw samples from those cdfs. We see that this approach is slower than using the built-in samplers, e.g., the normal sampler.
5. In the [5th example](examples/05_sample_from_cdf_beta/example.c), we define the cdf for the beta distribution, and we draw samples from it.
6. In the [6th example](examples/06_gamma_beta/example.c), we take samples from simple gamma and beta distributions, using the samplers provided by this library.
7. In the [7th example](examples/07_ci_beta/example.c), we get the 90% confidence interval of a beta distribution
8. The [8th example](examples/08_nuclear_war/example.c) translates the models from Eli and Nuño from [Samotsvety Nuclear Risk Forecasts — March 2022](https://forum.nunosempere.com/posts/KRFXjCqqfGQAYirm5/samotsvety-nuclear-risk-forecasts-march-2022#Nu_o_Sempere) into squiggle.c, then creates a mixture from both, and returns the mean probability of death per month and the 90% confidence interval.
8. The [9th example](examples/09_burn_10kg_fat/example.c) estimates how many minutes per day I would have to jump rope in order to lose 10kg of fat in half a year.
## Commentary
### squiggle.c is short
[squiggle.c](squiggle.c) is less than 600 lines of C, with a core of <250 lines. The reader could just read it and grasp its contents.
### Core strategy
This library provides some basic building blocks. The recommended strategy is to:
1. Define sampler functions, which take a seed, and return 1 sample
2. Compose those sampler functions to define your estimation model
3. At the end, call the last sampler function many times to generate many samples from your model
### Cdf auxiliary functions
### Nested functions and compilation with tcc.
GCC has an extension which allows a program to define a function inside another function. This makes squiggle.c code more linear and nicer to read, at the cost of becoming dependent on GCC and hence sacrificing portability and increasing compilation times. Conversely, compiling with tcc (tiny c compiler) is almost instantaneous, but leads to longer execution times and doesn't allow for nested functions.
| GCC | tcc |
| --- | --- |
| slower compilation | faster compilation |
| allows nested functions | doesn't allow nested functions |
| faster execution | slower execution |
My recommendation would be to use tcc while drawing a small number of samples for fast iteration, and then using gcc for the final version with lots of samples, and possibly with nested functions for ease of reading by others.
### Guarantees and licensing
- I offer no guarantees about stability, correctness, performance, etc. I might, for instance, abandon the version in C and rewrite it in Zig, Nim or Rust.
- This project mostly exists for my own usage & for my own amusement.
- Caution! Think carefully before using this project for anything important
- If you wanted to pay me to provide some stability or correctness, guarantees, or to tweak this library for your own usage, or to teach you how to use it, you could do so [here](https://nunosempere.com/consulting). Although this theoretical possibility exists, I don't I don't anticipate that this would be a good idea on most cases.
This project is released under the MIT license, a permissive open-source license. You can see it in the LICENSE.txt file.
### Design choices
This code should aim to be correct, then simple, then fast.
- It should be correct. The user should be able to rely on it and not think about whether errors come from the library.
- Nonetheless, the user should understand the limitations of sampling-based methods. See the section on [Tests and the long tail of the lognormal](https://git.nunosempere.com/personal/squiggle.c#tests-and-the-long-tail-of-the-lognormal) for a discussion of how sampling is bad at capturing some aspects of distributions with long tails.
- It should be clear, conceptually simple. Simple for me to implement, simple for others to understand.
- It should be fast. But when speed conflicts with simplicity, choose simplicity. For example, there might be several possible algorithms to sample a distribution, each of which is faster over part of the domain. In that case, it's conceptually simpler to just pick one algorithm, and pay the—normally small—performance penalty. In any case, though, the code should still be *way faster* than Python.
Note that being terse, or avoiding verbosity, is a non-goal. This is in part because of the constraints that C imposes. But it also aids with clarity and conceptual simplicity, as the issue of correlated samples illustrates in the next section.
### Correlated samples
In the original [squiggle](https://www.squiggle-language.com/) language, there is some ambiguity about what this code means:
```js
a = 1 to 10
b = 2 * a
c = b/a
c
```
Likewise in [squigglepy](https://github.com/rethinkpriorities/squigglepy):
```python
import squigglepy as sq
import numpy as np
a = sq.to(1, 3)
b = 2 * a
c = b / a
c_samples = sq.sample(c, 10)
print(c_samples)
```
Should `c` be equal to `2`? or should it be equal to 2 times the expected distribution of the ratio of two independent draws from a (`2 * a/a`, as it were)?
In squiggle.c, this ambiguity doesn't exist, at the cost of much greater overhead & verbosity:
```c
// correlated samples
// gcc -O3 correlated.c squiggle.c -lm -o correlated
#include "squiggle.h"
#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
int main(){
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with a seed of 0
double a = sample_to(1, 10, seed);
double b = 2 * a;
double c = b / a;
printf("a: %f, b: %f, c: %f\n", a, b, c);
// a: 0.607162, b: 1.214325, c: 0.500000
free(seed);
}
```
vs
```c
// uncorrelated samples
// gcc -O3 uncorrelated.c ../../squiggle.c -lm -o uncorrelated
#include "squiggle.h"
#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
double draw_xyz(uint64_t* seed){
// function could also be placed inside main with gcc nested functions extension.
return sample_to(1, 20, seed);
}
int main(){
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with a seed of 0
double a = draw_xyz(seed);
double b = 2 * draw_xyz(seed);
double c = b / a;
printf("a: %f, b: %f, c: %f\n", a, b, c);
// a: 0.522484, b: 10.283501, c: 19.681936
free(seed)
}
```
### Tests and the long tail of the lognormal
Distribution functions can be tested with:
```sh
cd tests
make && make run
```
`make verify` is an alias that runs all the tests and just displays the ones that are failing.
These tests are somewhat rudimentary: they get between 1M and 10M samples from a given sampling function, and check that their mean and standard deviations correspond to what they should theoretically should be.
If you run `make run` (or `make verify`), you will see errors such as these:
```
[-] Mean test for normal(47211.047473, 682197.019012) NOT passed.
Mean of normal(47211.047473, 682197.019012): 46933.673278, vs expected mean: 47211.047473
delta: -277.374195, relative delta: -0.005910
[-] Std test for lognormal(4.584666, 2.180816) NOT passed.
Std of lognormal(4.584666, 2.180816): 11443.588861, vs expected std: 11342.434900
delta: 101.153961, relative delta: 0.008839
[-] Std test for to(13839.861856, 897828.354318) NOT passed.
Std of to(13839.861856, 897828.354318): 495123.630575, vs expected std: 498075.002499
delta: -2951.371925, relative delta: -0.005961
```
These tests I wouldn't worry about. Due to luck of the draw, their relative error is a bit over 0.005, or 0.5%, and so the test fails. But it would surprise me if that had some meaningful practical implication.
The errors that should raise some worry are:
```
[-] Mean test for lognormal(1.210013, 4.766882) NOT passed.
Mean of lognormal(1.210013, 4.766882): 342337.257677, vs expected mean: 288253.061628
delta: 54084.196049, relative delta: 0.157985
[-] Std test for lognormal(1.210013, 4.766882) NOT passed.
Std of lognormal(1.210013, 4.766882): 208107782.972184, vs expected std: 24776840217.604111
delta: -24568732434.631927, relative delta: -118.057730
[-] Mean test for lognormal(-0.195240, 4.883106) NOT passed.
Mean of lognormal(-0.195240, 4.883106): 87151.733198, vs expected mean: 123886.818303
delta: -36735.085104, relative delta: -0.421507
[-] Std test for lognormal(-0.195240, 4.883106) NOT passed.
Std of lognormal(-0.195240, 4.883106): 33837426.331671, vs expected std: 18657000192.914921
delta: -18623162766.583248, relative delta: -550.371727
[-] Mean test for lognormal(0.644931, 4.795860) NOT passed.
Mean of lognormal(0.644931, 4.795860): 125053.904456, vs expected mean: 188163.894101
delta: -63109.989645, relative delta: -0.504662
[-] Std test for lognormal(0.644931, 4.795860) NOT passed.
Std of lognormal(0.644931, 4.795860): 39976300.711166, vs expected std: 18577298706.170452
delta: -18537322405.459286, relative delta: -463.707799
```
What is happening in this case is that you are taking a normal, like `normal(-0.195240, 4.883106)`, and you are exponentiating it to arrive at a lognormal. But `normal(-0.195240, 4.883106)` is going to have some noninsignificant weight on, say, 18. But `exp(18) = 39976300`, and points like it are going to end up a nontrivial amount to the analytical mean and standard deviation, even though they have little probability mass.
The reader can also check that for more plausible real-world values, like those fitting a lognormal to a really wide 90% confidence interval from 10 to 10k, errors aren't eggregious:
```
[x] Mean test for to(10.000000, 10000.000000) PASSED
[-] Std test for to(10.000000, 10000.000000) NOT passed.
Std of to(10.000000, 10000.000000): 23578.091775, vs expected std: 25836.381819
delta: -2258.290043, relative delta: -0.095779
```
Overall, I would caution that if you really care about the very far tails of distributions, you might want to instead use tools which can do some of the analytical manipulations for you, like the original Squiggle, Simple Squiggle (both linked below), or even doing lognormal multiplication by hand, relying on the fact that two lognormals multiplied together result in another lognormal with known shape.
In fact, squiggle.c does have a few functions for algebraic manipulations of simple distributions at the end of squiggle.c. But these are pretty rudimentary, and I don't know whether I'll end up expanding or deleting them.
### Results of running clang-tidy
clang-tidy is a utility to detect common errors in C/C++. You can run it with:
```
make tidy
```
It emits one warning about something I already took care of, so by default I've suppressed it. I think this is good news in terms of making me more confident that this simple library is correct :).
### Division between core functions and squiggle_moreneous expansions
This library differentiates between core functions, which are pretty tightly scoped, and expansions and convenience functions, which are more meandering. Expansions are in `squiggle_more.c` and `squiggle_more.h`. To use them, take care to link them:
```
// In your C source file
#include "squiggle_more.h"
```
```
# When compiling:
gcc -std=c99 -Wall -O3 example.c squiggle.c squiggle_more.c -lm -o ./example
```
#### Extra: Cdf auxiliary functions
I provide some Take a cdf, and return a sample from the distribution produced by that cdf. This might make it easier to program models, at the cost of a 20x to 60x slowdown in the parts of the code that use it.
#### Extra: Error propagation vs exiting on error
The process of taking a cdf and returning a sample might fail, e.g., it's a Newton method which might fail to converge because of cdf artifacts. The cdf itself might also fail, e.g., if a distribution only accepts a range of parameters, but is fed parameters outside that range.
This library provides two approaches:
1. Print the line and function in which the error occured, then exit on error
2. In situations where there might be an error, return a struct containing either the correct value or an error message:
```C
struct box {
int empty;
double content;
char* error_msg;
};
```
The first approach produces terser programs but might not scale. The second approach seems like it could lead to more robust programmes, but is more verbose.
Behaviour on error can be toggled by the `EXIT_ON_ERROR` variable. This library also provides a convenient macro, `PROCESS_ERROR`, to make error handling in either case much terser—see the usage in example 4 in the examples/ folder.
Overall, I'd describe the error handling capabilities of this library as pretty rudimentary. For example, this program might fail in surprising ways if you ask for a lognormal with negative standard deviation, because I haven't added error checking for that case yet.
## Extra: confidence intervals
// to do
## Extra paralellism
// to do
## Related projects
- [Squiggle](https://www.squiggle-language.com/)
- [SquigglePy](https://github.com/rethinkpriorities/squigglepy)
- [Simple Squiggle](https://nunosempere.com/blog/2022/04/17/simple-squiggle/)
- [time to botec](https://github.com/NunoSempere/time-to-botec)
- [beta]()
## To do list
- [ ] Document paralellism
- [ ] Document confidence intervals
- [ ] Point out that, even though the C standard is ambiguous about this, this code assumes that doubles are 64 bit precision (otherwise the xorshift should be different).
- [ ] Document rudimentary algebra manipulations for normal/lognormal
- [ ] Think through whether to delete cdf => samples function
- [ ] Think through whether to:
- simplify and just abort on error
- complexify and use boxes for everything
- leave as is
- [ ] Systematize references
- [ ] Support all distribution functions in <https://www.squiggle-language.com/docs/Api/Dist>
- [ ] do so efficiently
- [ ] Add more functions to do algebra and get the 90% c.i. of normals, lognormals, betas, etc.
- Think through which of these make sense.
- [ ] Disambiguate sample_laplace--successes vs failures || successes vs total trials as two distinct and differently named functions
## Done
- [x] Add example for only one sample
- [x] Add example for many samples
- [ ] ~~Add a custom preprocessor to allow simple nested functions that don't rely on local scope?~~
- [x] Use gcc extension to define functions nested inside main.
- [x] Chain various `sample_mixture` functions
- [x] Add beta distribution
- See <https://stats.stackexchange.com/questions/502146/how-does-numpy-generate-samples-from-a-beta-distribution> for a faster method.
- [ ] ~~Use OpenMP for acceleration~~
- [x] Add function to get sample when given a cdf
- [x] Don't have a single header file.
- [x] Structure project a bit better
- [x] Simplify `PROCESS_ERROR` macro
- [x] Add README
- [x] Schema: a function which takes a sample and manipulates it,
- [x] and at the end, an array of samples.
- [x] Explain boxes
- [x] Explain nested functions
- [x] Explain exit on error
- [x] Explain individual examples
- [x] Rename functions to something more self-explanatory, e.g,. `sample_unit_normal`.
- [x] Add summarization functions: mean, std
- [x] Add sampling from a gamma distribution
- https://dl.acm.org/doi/pdf/10.1145/358407.358414
- [x] Explain correlated samples
- [ ] ~~Add tests in Stan?~~
- [x] Test summary statistics for each of the distributions.
- [x] For uniform
- [x] For normal
- [x] For lognormal
- [x] For lognormal (to syntax)
- [x] For beta distribution
- [x] Clarify gamma/standard gamma
- [x] Add efficient sampling from a beta distribution
- https://dl.acm.org/doi/10.1145/358407.358414
- https://link.springer.com/article/10.1007/bf02293108
- https://stats.stackexchange.com/questions/502146/how-does-numpy-generate-samples-from-a-beta-distribution
- https://github.com/numpy/numpy/blob/5cae51e794d69dd553104099305e9f92db237c53/numpy/random/src/distributions/distributions.c
- [x] Pontificate about lognormal tests
- [x] Give warning about sampling-based methods.
- [x] Have some more complicated & realistic example
- [x] Add summarization functions: 90% ci (or all c.i.?)
- [x] Link to the examples in the examples section.
- [x] Add a few functions for doing simple algebra on normals, and lognormals
- [x] Add prototypes
- [x] Use named structs
- [x] Add to header file
- [x] Provide example algebra
- [x] Add conversion between 90% ci and parameters.
- [x] Use that conversion in conjuction with small algebra.
- [x] Consider ergonomics of using ci instead of c_i
- [x] use named struct instead
- [x] demonstrate and document feeding a struct directly to a function; my_function((struct c_i){.low = 1, .high = 2});
- [ ] Consider desirability of defining shortcuts for those functions. Adds a level of magic, though.
- [ ] Test results
- [x] Move to own file? Or signpost in file? => signposted in file.
- [x] Write twitter thread: now [here](https://twitter.com/NunoSempere/status/1707041153210564959); retweets appreciated.
- [ ] ~~Think about whether to write a simple version of this for [uxn](https://100r.co/site/uxn.html), a minimalistic portable programming stack which, sadly, doesn't have doubles (64 bit floats)~~

BIN
squiggle.c/squiggle_c/examples/core/00_example_template/example
View File

Binary file not shown.

14
squiggle.c/squiggle_c/examples/core/00_example_template/example.c
View File

@ -1,14 +0,0 @@
#include "../../../squiggle.h"
#include <stdio.h>
#include <stdlib.h>
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
// ...
free(seed);
}

BIN
squiggle.c/squiggle_c/examples/core/01_one_sample/example
View File

Binary file not shown.

43
squiggle.c/squiggle_c/examples/core/01_one_sample/example.c
View File

@ -1,43 +0,0 @@
#include "../../../squiggle.h"
#include <stdio.h>
#include <stdlib.h>
// Estimate functions
double sample_0(uint64_t* seed)
{
return 0;
}
double sample_1(uint64_t* seed)
{
return 1;
}
double sample_few(uint64_t* seed)
{
return sample_to(1, 3, seed);
}
double sample_many(uint64_t* seed)
{
return sample_to(2, 10, seed);
}
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
double p_a = 0.8;
double p_b = 0.5;
double p_c = p_a * p_b;
int n_dists = 4;
double weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
double (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many };
double result_one = sample_mixture(samplers, weights, n_dists, seed);
printf("result_one: %f\n", result_one);
free(seed);
}

BIN
squiggle.c/squiggle_c/examples/core/02_time_to_botec/example
View File

Binary file not shown.

32
squiggle.c/squiggle_c/examples/core/02_time_to_botec/example.c
View File

@ -1,32 +0,0 @@
#include "../../../squiggle.h"
#include <stdio.h>
#include <stdlib.h>
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
double p_a = 0.8;
double p_b = 0.5;
double p_c = p_a * p_b;
double sample_0(uint64_t* seed){ return 0; }
double sample_1(uint64_t* seed) { return 1; }
double sample_few(uint64_t* seed) { return sample_to(1, 3, seed); }
double sample_many(uint64_t* seed) { return sample_to(2, 10, seed); }
int n_dists = 4;
double weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
double (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many };
int n_samples = 1000000;
double* result_many = (double*)malloc(n_samples * sizeof(double));
for (int i = 0; i < n_samples; i++) {
result_many[i] = sample_mixture(samplers, weights, n_dists, seed);
}
printf("Mean: %f\n", array_mean(result_many, n_samples));
free(seed);
}

BIN
squiggle.c/squiggle_c/examples/core/03_gcc_nested_function/example
View File

Binary file not shown.

38
squiggle.c/squiggle_c/examples/core/03_gcc_nested_function/example.c
View File

@ -1,38 +0,0 @@
#include "../../../squiggle.h"
#include <stdio.h>
#include <stdlib.h>
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
double p_a = 0.8;
double p_b = 0.5;
double p_c = p_a * p_b;
int n_dists = 4;
// These are nested functions. They will not compile without gcc.
double sample_0(uint64_t * seed) { return 0; }
double sample_1(uint64_t * seed) { return 1; }
double sample_few(uint64_t * seed) { return sample_to(1, 3, seed); }
double sample_many(uint64_t * seed) { return sample_to(2, 10, seed); }
double (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many };
double weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
int n_samples = 1000000;
double* result_many = (double*)malloc(n_samples * sizeof(double));
for (int i = 0; i < n_samples; i++) {
result_many[i] = sample_mixture(samplers, weights, n_dists, seed);
}
printf("result_many: [");
for (int i = 0; i < 100; i++) {
printf("%.2f, ", result_many[i]);
}
printf("]\n");
free(seed);
}

BIN
squiggle.c/squiggle_c/examples/core/04_gamma_beta/example
View File

Binary file not shown.

43
squiggle.c/squiggle_c/examples/core/04_gamma_beta/example.c
View File

@ -1,43 +0,0 @@
#include "../../../squiggle.h"
#include <stdio.h>
#include <stdlib.h>
// Estimate functions
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
int n = 1000 * 1000;
/*
for (int i = 0; i < n; i++) {
double gamma_0 = sample_gamma(0.0, seed);
// printf("sample_gamma(0.0): %f\n", gamma_0);
}
printf("\n");
*/
double* gamma_1_array = malloc(sizeof(double) * n);
for (int i = 0; i < n; i++) {
double gamma_1 = sample_gamma(1.0, seed);
// printf("sample_gamma(1.0): %f\n", gamma_1);
gamma_1_array[i] = gamma_1;
}
printf("gamma(1) summary statistics = mean: %f, std: %f\n", array_mean(gamma_1_array, n), array_std(gamma_1_array, n));
free(gamma_1_array);
printf("\n");
double* beta_1_2_array = malloc(sizeof(double) * n);
for (int i = 0; i < n; i++) {
double beta_1_2 = sample_beta(1, 2.0, seed);
// printf("sample_beta(1.0, 2.0): %f\n", beta_1_2);
beta_1_2_array[i] = beta_1_2;
}
printf("beta(1,2) summary statistics: mean: %f, std: %f\n", array_mean(beta_1_2_array, n), array_std(beta_1_2_array, n));
free(beta_1_2_array);
printf("\n");
free(seed);
}

BIN
squiggle.c/squiggle_c/examples/core/05_hundred_lognormals/example
View File

Binary file not shown.

17
squiggle.c/squiggle_c/examples/core/05_hundred_lognormals/example.c
View File

@ -1,17 +0,0 @@
#include "../../../squiggle.h"
#include <stdio.h>
#include <stdlib.h>
// Estimate functions
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
for (int i = 0; i < 100; i++) {
double sample = sample_lognormal(0, 10, seed);
printf("%f\n", sample);
}
free(seed);
}

1
squiggle.c/squiggle_c/examples/core/05_hundred_lognormals/run-sorted.sh
View File

@ -1 +0,0 @@
./example | sort -h

85
squiggle.c/squiggle_c/examples/core/makefile
View File

@ -1,85 +0,0 @@
# Interface:
# make all
# make format-all
# make run-all
# make one DIR=01_one_sample
# make format-one DIR=01_one_sample
# make run-one DIR=01_one_sample
# make time-linux-one DIR=01_one_sample
# make profile-one DIR=01_one_sample
# Compiler
CC=gcc
# CC=tcc # <= faster compilation
# Main file
SRC=example.c
OUTPUT=example
## Dependencies
SQUIGGLE=../../squiggle.c
MATH=-lm
DEPS=$(SQUIGGLE) $(MATH)
## Flags
DEBUG= #'-g'
STANDARD=-std=c99
WARNINGS=-Wall
OPTIMIZED=-O3 #-Ofast
## Formatter
STYLE_BLUEPRINT=webkit
FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT)
## make all
all:
$(CC) $(OPTIMIZED) $(DEBUG) 00_example_template/$(SRC) $(DEPS) -o 00_example_template/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 01_one_sample/$(SRC) $(DEPS) -o 01_one_sample/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 02_time_to_botec/$(SRC) $(DEPS) -o 02_time_to_botec/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 03_gcc_nested_function/$(SRC) $(DEPS) -o 03_gcc_nested_function/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 04_gamma_beta/$(SRC) $(DEPS) -o 04_gamma_beta/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 05_hundred_lognormals/$(SRC) $(DEPS) -o 05_hundred_lognormals/$(OUTPUT)
format-all:
$(FORMATTER) 00_example_template/$(SRC)
$(FORMATTER) 01_one_sample/$(SRC)
$(FORMATTER) 02_time_to_botec/$(SRC)
$(FORMATTER) 03_gcc_nested_function/$(SRC)
$(FORMATTER) 04_gamma_beta/$(SRC)
$(FORMATTER) 05_hundred_lognormals/$(SRC)
run-all:
00_example_template/$(OUTPUT)
01_one_sample/$(OUTPUT)
02_time_to_botec/$(OUTPUT)
03_gcc_nested_function/$(OUTPUT)
04_gamma_beta/$(OUTPUT)
05_hundred_lognormals/$(OUTPUT)
## make one DIR=01_one_sample
one: $(DIR)/$(SRC)
$(CC) $(OPTIMIZED) $(DEBUG) $(DIR)/$(SRC) $(DEPS) -o $(DIR)/$(OUTPUT)
## make format-one DIR=01_one_sample
format-one: $(DIR)/$(SRC)
$(FORMATTER) $(DIR)/$(SRC)
## make run-one DIR=01_one_sample
run-one: $(DIR)/$(OUTPUT)
$(DIR)/$(OUTPUT) && echo
## make time-linux-one DIR=01_one_sample
time-linux-one: $(DIR)/$(OUTPUT)
@echo "Requires /bin/time, found on GNU/Linux systems" && echo
@echo "Running 100x and taking avg time $(DIR)/$(OUTPUT)"
@t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do $(DIR)/$(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 1 thread: |" | sed 's|$$|ms|' && echo
## e.g., make profile-linux-one DIR=01_one_sample
profile-linux-one:
echo "Requires perf, which depends on the kernel version, and might be in linux-tools package or similar"
echo "Must be run as sudo"
$(CC) $(OPTIMIZED) $(DEBUG) $(DIR)/$(SRC) $(DEPS) -o $(DIR)/$(OUTPUT)
# $(CC) $(SRC) $(DEPS) -o $(OUTPUT)
sudo perf record $(DIR)/$(OUTPUT)
sudo perf report
rm perf.data

BIN
squiggle.c/squiggle_c/examples/more/00_example_template/example
View File

Binary file not shown.

15
squiggle.c/squiggle_c/examples/more/00_example_template/example.c
View File

@ -1,15 +0,0 @@
#include "../../../squiggle.h"
#include "../../../squiggle_more.h"
#include <stdio.h>
#include <stdlib.h>
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
// ...
free(seed);
}

BIN
squiggle.c/squiggle_c/examples/more/01_sample_from_cdf/example
View File

Binary file not shown.

102
squiggle.c/squiggle_c/examples/more/01_sample_from_cdf/example.c
View File

@ -1,102 +0,0 @@
#include "../../../squiggle.h"
#include "../../../squiggle_more.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define NUM_SAMPLES 1000000
// Example cdf
double cdf_uniform_0_1(double x)
{
if (x < 0) {
return 0;
} else if (x > 1) {
return 1;
} else {
return x;
}
}
double cdf_squared_0_1(double x)
{
if (x < 0) {
return 0;
} else if (x > 1) {
return 1;
} else {
return x * x;
}
}
double cdf_normal_0_1(double x)
{
double mean = 0;
double std = 1;
return 0.5 * (1 + erf((x - mean) / (std * sqrt(2)))); // erf from math.h
}
// Some testers
void test_inverse_cdf_double(char* cdf_name, double cdf_double(double))
{
struct box result = inverse_cdf_double(cdf_double, 0.5);
if (result.empty) {
printf("Inverse for %s not calculated\n", cdf_name);
exit(1);
} else {
printf("Inverse of %s at %f is: %f\n", cdf_name, 0.5, result.content);
}
}
void test_and_time_sampler_double(char* cdf_name, double cdf_double(double), uint64_t* seed)
{
printf("\nGetting some samples from %s:\n", cdf_name);
clock_t begin = clock();
for (int i = 0; i < NUM_SAMPLES; i++) {
struct box sample = sampler_cdf_double(cdf_double, seed);
if (sample.empty) {
printf("Error in sampler function for %s", cdf_name);
} else {
// printf("%f\n", sample.content);
}
}
clock_t end = clock();
double time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
printf("Time spent: %f\n", time_spent);
}
int main()
{
// Test inverse cdf double
test_inverse_cdf_double("cdf_uniform_0_1", cdf_uniform_0_1);
test_inverse_cdf_double("cdf_squared_0_1", cdf_squared_0_1);
test_inverse_cdf_double("cdf_normal_0_1", cdf_normal_0_1);
// Testing samplers
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
// Test double sampler
test_and_time_sampler_double("cdf_uniform_0_1", cdf_uniform_0_1, seed);
test_and_time_sampler_double("cdf_squared_0_1", cdf_squared_0_1, seed);
test_and_time_sampler_double("cdf_normal_0_1", cdf_normal_0_1, seed);
// Get some normal samples using a previous approach
printf("\nGetting some samples from sample_unit_normal\n");
clock_t begin_2 = clock();
double* normal_samples = malloc(NUM_SAMPLES * sizeof(double));
for (int i = 0; i < NUM_SAMPLES; i++) {
normal_samples[i] = sample_unit_normal(seed);
// printf("%f\n", normal_sample);
}
clock_t end_2 = clock();
double time_spent_2 = (double)(end_2 - begin_2) / CLOCKS_PER_SEC;
printf("Time spent: %f\n", time_spent_2);
free(seed);
return 0;
}

BIN
squiggle.c/squiggle_c/examples/more/02_sample_from_cdf_beta/example
View File

Binary file not shown.

168
squiggle.c/squiggle_c/examples/more/02_sample_from_cdf_beta/example.c
View File

@ -1,168 +0,0 @@
#include "../../../squiggle.h"
#include "../../../squiggle_more.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define NUM_SAMPLES 10000
#define STOP_BETA 1.0e-8
#define TINY_BETA 1.0e-30
// Incomplete beta function
struct box incbeta(double a, double b, double x)
{
// Descended from <https://github.com/codeplea/incbeta/blob/master/incbeta.c>,
// <https://codeplea.com/incomplete-beta-function-c>
// but modified to return a box struct and doubles instead of doubles.
// [ ] to do: add attribution in README
// Original code under this license:
/*
* zlib License
*
* Regularized Incomplete Beta Function
*
* Copyright (c) 2016, 2017 Lewis Van Winkle
* http://CodePlea.com
*
* This software is provided 'as-is', without any express or implied
* warranty. In no event will the authors be held liable for any damages
* arising from the use of this software.
*
* Permission is granted to anyone to use this software for any purpose,
* including commercial applications, and to alter it and redistribute it
* freely, subject to the following restrictions:
*
* 1. The origin of this software must not be misrepresented; you must not
* claim that you wrote the original software. If you use this software
* in a product, an acknowledgement in the product documentation would be
* appreciated but is not required.
* 2. Altered source versions must be plainly marked as such, and must not be
* misrepresented as being the original software.
* 3. This notice may not be removed or altered from any source distribution.
*/
if (x < 0.0 || x > 1.0) {
return PROCESS_ERROR("x out of bounds [0, 1], in function incbeta");
}
/*The continued fraction converges nicely for x < (a+1)/(a+b+2)*/
if (x > (a + 1.0) / (a + b + 2.0)) {
struct box symmetric_incbeta = incbeta(b, a, 1.0 - x);
if (symmetric_incbeta.empty) {
return symmetric_incbeta; // propagate error
} else {
struct box result = {
.empty = 0,
.content = 1 - symmetric_incbeta.content
};
return result;
}
}
/*Find the first part before the continued fraction.*/
const double lbeta_ab = lgamma(a) + lgamma(b) - lgamma(a + b);
const double front = exp(log(x) * a + log(1.0 - x) * b - lbeta_ab) / a;
/*Use Lentz's algorithm to evaluate the continued fraction.*/
double f = 1.0, c = 1.0, d = 0.0;
int i, m;
for (i = 0; i <= 200; ++i) {
m = i / 2;
double numerator;
if (i == 0) {
numerator = 1.0; /*First numerator is 1.0.*/
} else if (i % 2 == 0) {
numerator = (m * (b - m) * x) / ((a + 2.0 * m - 1.0) * (a + 2.0 * m)); /*Even term.*/
} else {
numerator = -((a + m) * (a + b + m) * x) / ((a + 2.0 * m) * (a + 2.0 * m + 1)); /*Odd term.*/
}
/*Do an iteration of Lentz's algorithm.*/
d = 1.0 + numerator * d;
if (fabs(d) < TINY_BETA)
d = TINY_BETA;
d = 1.0 / d;
c = 1.0 + numerator / c;
if (fabs(c) < TINY_BETA)
c = TINY_BETA;
const double cd = c * d;
f *= cd;
/*Check for stop.*/
if (fabs(1.0 - cd) < STOP_BETA) {
struct box result = {
.empty = 0,
.content = front * (f - 1.0)
};
return result;
}
}
return PROCESS_ERROR("More loops needed, did not converge, in function incbeta");
}
struct box cdf_beta(double x)
{
if (x < 0) {
struct box result = { .empty = 0, .content = 0 };
return result;
} else if (x > 1) {
struct box result = { .empty = 0, .content = 1 };
return result;
} else {
double successes = 1, failures = (2023 - 1945);
return incbeta(successes, failures, x);
}
}
// Some testers
void test_inverse_cdf_box(char* cdf_name, struct box cdf_box(double))
{
struct box result = inverse_cdf_box(cdf_box, 0.5);
if (result.empty) {
printf("Inverse for %s not calculated\n", cdf_name);
exit(1);
} else {
printf("Inverse of %s at %f is: %f\n", cdf_name, 0.5, result.content);
}
}
void test_and_time_sampler_box(char* cdf_name, struct box cdf_box(double), uint64_t* seed)
{
printf("\nGetting some samples from %s:\n", cdf_name);
clock_t begin = clock();
for (int i = 0; i < NUM_SAMPLES; i++) {
struct box sample = sampler_cdf_box(cdf_box, seed);
if (sample.empty) {
printf("Error in sampler function for %s", cdf_name);
} else {
// printf("%f\n", sample.content);
}
}
clock_t end = clock();
double time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
printf("Time spent: %f\n", time_spent);
}
int main()
{
// Test inverse cdf box
test_inverse_cdf_box("cdf_beta", cdf_beta);
// Test box sampler
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
test_and_time_sampler_box("cdf_beta", cdf_beta, seed);
// Ok, this is slower than python!!
// Partly this is because I am using a more general algorithm,
// which applies to any cdf
// But I am also using absurdly precise convergence conditions.
// This could be optimized.
free(seed);
return 0;
}

BIN
squiggle.c/squiggle_c/examples/more/03_ci_beta/example
View File

Binary file not shown.

22
squiggle.c/squiggle_c/examples/more/03_ci_beta/example.c
View File

@ -1,22 +0,0 @@
#include "../../../squiggle.h"
#include "../../../squiggle_more.h"
#include <stdio.h>
#include <stdlib.h>
// Estimate functions
double beta_1_2_sampler(uint64_t* seed)
{
return sample_beta(1, 2.0, seed);
}
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
ci beta_1_2_ci_90 = get_90_confidence_interval(beta_1_2_sampler, seed);
printf("90%% confidence interval of beta(1,2) is [%f, %f]\n", beta_1_2_ci_90.low, beta_1_2_ci_90.high);
free(seed);
}

BIN
squiggle.c/squiggle_c/examples/more/04_nuclear_war/example
View File

Binary file not shown.

68
squiggle.c/squiggle_c/examples/more/04_nuclear_war/example.c
View File

@ -1,68 +0,0 @@
#include "../../../squiggle.h"
#include "../../../squiggle_more.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
double probability_of_dying_nuno(uint64_t* seed)
{
double first_year_russian_nuclear_weapons = 1953;
double current_year = 2022;
double laplace_probability_nuclear_exchange_year = sample_beta(1, current_year - first_year_russian_nuclear_weapons + 1, seed);
double laplace_probability_nuclear_exchange_month = 1 - pow(1 - laplace_probability_nuclear_exchange_year, (1.0 / 12.0));
double london_hit_conditional_on_russia_nuclear_weapon_usage = sample_beta(7.67, 69.65, seed);
// I.e., a beta distribution with a range of 0.05 to 0.16 into: https://nunosempere.com/blog/2023/03/15/fit-beta/
// 0.05 were my estimate and Samotsvety's estimate in March 2022, respectively:
// https://forum.effectivealtruism.org/posts/KRFXjCqqfGQAYirm5/samotsvety-nuclear-risk-forecasts-march-2022#Nu_o_Sempere
double informed_actor_not_able_to_escape = sample_beta(3.26212166586967, 3.26228162008564, seed);
// 0.2 to 0.8, i.e., 20% to 80%, again using the previous tool
double proportion_which_die_if_bomb_drops_in_london = sample_beta(10.00, 2.45, seed); // 60% to 95%
double probability_of_dying = laplace_probability_nuclear_exchange_month * london_hit_conditional_on_russia_nuclear_weapon_usage * informed_actor_not_able_to_escape * proportion_which_die_if_bomb_drops_in_london;
return probability_of_dying;
}
double probability_of_dying_eli(uint64_t* seed)
{
double russia_nato_nuclear_exchange_in_next_month = sample_beta(1.30, 1182.99, seed); // .0001 to .003
double london_hit_conditional = sample_beta(3.47, 8.97, seed); // 0.1 to 0.5
double informed_actors_not_able_to_escape = sample_beta(2.73, 5.67, seed); // .1 to .6
double proportion_which_die_if_bomb_drops_in_london = sample_beta(3.00, 1.46, seed); // 0.3 to 0.95;
double probability_of_dying = russia_nato_nuclear_exchange_in_next_month * london_hit_conditional * informed_actors_not_able_to_escape * proportion_which_die_if_bomb_drops_in_london;
return probability_of_dying;
}
double mixture(uint64_t* seed)
{
double (*samplers[])(uint64_t*) = { probability_of_dying_nuno, probability_of_dying_eli };
double weights[] = { 0.5, 0.5 };
return sample_mixture(samplers, weights, 2, seed);
}
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
int n = 1000 * 1000;
double* mixture_result = malloc(sizeof(double) * n);
for (int i = 0; i < n; i++) {
mixture_result[i] = mixture(seed);
}
printf("mixture_result: [ ");
for (int i = 0; i < 9; i++) {
printf("%.6f, ", mixture_result[i]);
}
printf("... ]\n");
ci ci_90 = get_90_confidence_interval(mixture, seed);
printf("mean: %f\n", array_mean(mixture_result, n));
printf("90%% confidence interval: [%f, %f]\n", ci_90.low, ci_90.high);
free(seed);
}

BIN
squiggle.c/squiggle_c/examples/more/04_nuclear_war/scratchpad/example
View File

Binary file not shown.

20
squiggle.c/squiggle_c/examples/more/04_nuclear_war/scratchpad/example.c
View File

@ -1,20 +0,0 @@
#include "../../../squiggle.h"
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
double firstYearRussianNuclearWeapons = 1953;
double currentYear = 2023;
for(int i=0; i<10; i++){
double laplace_beta = sample_beta(currentYear-firstYearRussianNuclearWeapons + 1, 1, seed);
printf("%f\n", laplace_beta);
}
free(seed);
}

53
squiggle.c/squiggle_c/examples/more/04_nuclear_war/scratchpad/makefile
View File

@ -1,53 +0,0 @@
# Interface:
# make
# make build
# make format
# make run
# Compiler
CC=gcc
# CC=tcc # <= faster compilation
# Main file
SRC=example.c ../../../squiggle.c
OUTPUT=example
## Dependencies
MATH=-lm
## Flags
DEBUG= #'-g'
STANDARD=-std=c99
WARNINGS=-Wall
OPTIMIZED=-O3 #-Ofast
# OPENMP=-fopenmp
## Formatter
STYLE_BLUEPRINT=webkit
FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT)
## make build
build: $(SRC)
$(CC) $(OPTIMIZED) $(DEBUG) $(SRC) $(MATH) -o $(OUTPUT)
format: $(SRC)
$(FORMATTER) $(SRC)
run: $(SRC) $(OUTPUT)
OMP_NUM_THREADS=1 ./$(OUTPUT) && echo
time-linux:
@echo "Requires /bin/time, found on GNU/Linux systems" && echo
@echo "Running 100x and taking avg time $(OUTPUT)"
@t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do $(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 1 thread: |" | sed 's|$$|ms|' && echo
## Profiling
profile-linux:
echo "Requires perf, which depends on the kernel version, and might be in linux-tools package or similar"
echo "Must be run as sudo"
$(CC) $(SRC) $(MATH) -o $(OUTPUT)
sudo perf record ./$(OUTPUT)
sudo perf report
rm perf.data

BIN
squiggle.c/squiggle_c/examples/more/05_burn_10kg_fat/example
View File

Binary file not shown.

48
squiggle.c/squiggle_c/examples/more/05_burn_10kg_fat/example.c
View File

@ -1,48 +0,0 @@
#include "../../../squiggle.h"
#include "../../../squiggle_more.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
double sample_minutes_per_day_jumping_rope_needed_to_burn_10kg(uint64_t* seed)
{
double kcal_jumping_rope_minute = sample_to(15, 20, seed);
double kcal_jumping_rope_hour = kcal_jumping_rope_minute * 60;
double kcal_in_kg_of_fat = 7700;
double num_kg_of_fat_to_lose = 10;
double hours_jumping_rope_needed = kcal_in_kg_of_fat * num_kg_of_fat_to_lose / kcal_jumping_rope_hour;
double days_until_end_of_year = 152; // as of 2023-08-01
double hours_per_day = hours_jumping_rope_needed / days_until_end_of_year;
double minutes_per_day = hours_per_day * 60;
return minutes_per_day;
}
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
int n = 1000 * 1000;
double* result = malloc(sizeof(double) * n);
for (int i = 0; i < n; i++) {
result[i] = sample_minutes_per_day_jumping_rope_needed_to_burn_10kg(seed);
}
printf("## How many minutes per day do I have to jump rope to lose 10kg of fat by the end of the year?\n");
printf("Mean: %f\n", array_mean(result, n));
printf("A few samples: [ ");
for (int i = 0; i < 9; i++) {
printf("%.6f, ", result[i]);
}
printf("... ]\n");
ci ci_90 = get_90_confidence_interval(sample_minutes_per_day_jumping_rope_needed_to_burn_10kg, seed);
printf("90%% confidence interval: [%f, %f]\n", ci_90.low, ci_90.high);
free(seed);
}

BIN
squiggle.c/squiggle_c/examples/more/06_nuclear_recovery/example
View File

Binary file not shown.

85
squiggle.c/squiggle_c/examples/more/06_nuclear_recovery/example.c
View File

@ -1,85 +0,0 @@
#include "../../../squiggle.h"
#include "../../../squiggle_more.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
double yearly_probability_nuclear_collapse(double year, uint64_t* seed)
{
double successes = 0;
double failures = (year - 1960);
return sample_laplace(successes, failures, seed);
// ^ can change to (successes + 1)/(trials + 2)
// to get a probability,
// rather than sampling from a distribution over probabilities.
}
double yearly_probability_nuclear_collapse_2023(uint64_t* seed)
{
return yearly_probability_nuclear_collapse(2023, seed);
}
double yearly_probability_nuclear_collapse_after_recovery(double year, double rebuilding_period_length_years, uint64_t* seed)
{
// assumption: nuclear
double successes = 1.0;
double failures = (year - rebuilding_period_length_years - 1960 - 1);
return sample_laplace(successes, failures, seed);
}
double yearly_probability_nuclear_collapse_after_recovery_example(uint64_t* seed)
{
double year = 2070;
double rebuilding_period_length_years = 30;
// So, there was a nuclear collapse in 2040,
// then a recovery period of 30 years
// and it's now 2070
return yearly_probability_nuclear_collapse_after_recovery(year, rebuilding_period_length_years, seed);
}
double yearly_probability_nuclear_collapse_after_recovery_antiinductive(uint64_t* seed)
{
return yearly_probability_nuclear_collapse(2023, seed) / 2;
}
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
int num_samples = 1000000;
// Before a first nuclear collapse
printf("## Before the first nuclear collapse\n");
ci ci_90_2023 = get_90_confidence_interval(yearly_probability_nuclear_collapse_2023, seed);
printf("90%% confidence interval: [%f, %f]\n", ci_90_2023.low, ci_90_2023.high);
double* yearly_probability_nuclear_collapse_2023_samples = malloc(sizeof(double) * num_samples);
for (int i = 0; i < num_samples; i++) {
yearly_probability_nuclear_collapse_2023_samples[i] = yearly_probability_nuclear_collapse_2023(seed);
}
printf("mean: %f\n", array_mean(yearly_probability_nuclear_collapse_2023_samples, num_samples));
// After the first nuclear collapse
printf("\n## After the first nuclear collapse\n");
ci ci_90_2070 = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_example, seed);
printf("90%% confidence interval: [%f, %f]\n", ci_90_2070.low, ci_90_2070.high);
double* yearly_probability_nuclear_collapse_after_recovery_samples = malloc(sizeof(double) * num_samples);
for (int i = 0; i < num_samples; i++) {
yearly_probability_nuclear_collapse_after_recovery_samples[i] = yearly_probability_nuclear_collapse_after_recovery_example(seed);
}
printf("mean: %f\n", array_mean(yearly_probability_nuclear_collapse_after_recovery_samples, num_samples));
// After the first nuclear collapse (antiinductive)
printf("\n## After the first nuclear collapse (antiinductive)\n");
ci ci_90_antiinductive = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_antiinductive, seed);
printf("90%% confidence interval: [%f, %f]\n", ci_90_antiinductive.low, ci_90_antiinductive.high);
double* yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples = malloc(sizeof(double) * num_samples);
for (int i = 0; i < num_samples; i++) {
yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples[i] = yearly_probability_nuclear_collapse_after_recovery_antiinductive(seed);
}
printf("mean: %f\n", array_mean(yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples, num_samples));
free(seed);
}

BIN
squiggle.c/squiggle_c/examples/more/07_algebra/example
View File

Binary file not shown.

26
squiggle.c/squiggle_c/examples/more/07_algebra/example.c
View File

@ -1,26 +0,0 @@
#include "../../../squiggle.h"
#include "../../../squiggle_more.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
normal_params n1 = { .mean = 1.0, .std = 3.0 };
normal_params n2 = { .mean = 2.0, .std = 4.0 };
normal_params sn = algebra_sum_normals(n1, n2);
printf("The sum of Normal(%f, %f) and Normal(%f, %f) is Normal(%f, %f)\n",
n1.mean, n1.std, n2.mean, n2.std, sn.mean, sn.std);
lognormal_params ln1 = { .logmean = 1.0, .logstd = 3.0 };
lognormal_params ln2 = { .logmean = 2.0, .logstd = 4.0 };
lognormal_params sln = algebra_product_lognormals(ln1, ln2);
printf("The product of Lognormal(%f, %f) and Lognormal(%f, %f) is Lognormal(%f, %f)\n",
ln1.logmean, ln1.logstd, ln2.logmean, ln2.logstd, sln.logmean, sln.logstd);
free(seed);
}

BIN
squiggle.c/squiggle_c/examples/more/08_algebra_and_conversion/example
View File

Binary file not shown.

33
squiggle.c/squiggle_c/examples/more/08_algebra_and_conversion/example.c
View File

@ -1,33 +0,0 @@
#include "../../../squiggle.h"
#include "../../../squiggle_more.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
// Convert to 90% confidence interval form and back
lognormal_params ln1 = { .logmean = 1.0, .logstd = 3.0 };
ci ln1_ci = convert_lognormal_params_to_ci(ln1);
printf("The 90%% confidence interval of Lognormal(%f, %f) is [%f, %f]\n",
ln1.logmean, ln1.logstd,
ln1_ci.low, ln1_ci.high);
lognormal_params ln1_params2 = convert_ci_to_lognormal_params(ln1_ci);
printf("The lognormal which has 90%% confidence interval [%f, %f] is Lognormal(%f, %f)\n",
ln1_ci.low, ln1_ci.high,
ln1_params2.logmean, ln1_params2.logstd);
lognormal_params ln2 = convert_ci_to_lognormal_params((ci) { .low = 1, .high = 10 });
lognormal_params ln3 = convert_ci_to_lognormal_params((ci) { .low = 5, .high = 50 });
lognormal_params sln = algebra_product_lognormals(ln2, ln3);
ci sln_ci = convert_lognormal_params_to_ci(sln);
printf("Result of some lognormal products: to(%f, %f)\n", sln_ci.low, sln_ci.high);
free(seed);
}

BIN
squiggle.c/squiggle_c/examples/more/09_ergonomic_algebra/example
View File

Binary file not shown.

26
squiggle.c/squiggle_c/examples/more/09_ergonomic_algebra/example.c
View File

@ -1,26 +0,0 @@
#include "../../../squiggle.h"
#include "../../../squiggle_more.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define ln lognormal_params
#define to(...) convert_ci_to_lognormal_params((ci)__VA_ARGS__)
#define from(...) convert_lognormal_params_to_ci((ln)__VA_ARGS__)
#define times(a, b) algebra_product_lognormals(a, b)
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
ln a = to({ .low = 1, .high = 10 });
ln b = to({ .low = 5, .high = 500 });
ln c = times(a, b);
printf("Result: to(%f, %f)\n", from(c).low, from(c).high);
printf("One sample from it is: %f\n", sample_lognormal(c.logmean, c.logstd, seed));
free(seed);
}

BIN
squiggle.c/squiggle_c/examples/more/10_twitter_thread_example/example
View File

Binary file not shown.

47
squiggle.c/squiggle_c/examples/more/10_twitter_thread_example/example.c
View File

@ -1,47 +0,0 @@
#include "../../../squiggle.h"
#include "../../../squiggle_more.h"
#include <stdio.h>
#include <stdlib.h>
double sample_0(uint64_t* seed)
{
return 0;
}
double sample_1(uint64_t* seed)
{
return 1;
}
double sample_normal_mean_1_std_2(uint64_t* seed)
{
return sample_normal(1, 2, seed);
}
double sample_1_to_3(uint64_t* seed)
{
return sample_to(1, 3, seed);
}
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
int n_dists = 4;
double weights[] = { 1, 2, 3, 4 };
double (*samplers[])(uint64_t*) = {
sample_0,
sample_1,
sample_normal_mean_1_std_2,
sample_1_to_3
};
int n_samples = 10;
for (int i = 0; i < n_samples; i++) {
printf("Sample #%d: %f\n", i, sample_mixture(samplers, weights, n_dists, seed));
}
free(seed);
}

BIN
squiggle.c/squiggle_c/examples/more/11_billion_lognormals_paralell/example
View File

Binary file not shown.

29
squiggle.c/squiggle_c/examples/more/11_billion_lognormals_paralell/example.c
View File

@ -1,29 +0,0 @@
#include "../../../squiggle.h"
#include "../../../squiggle_more.h"
#include <stdio.h>
#include <stdlib.h>
// Estimate functions
int main()
{
// set randomness seed
// uint64_t* seed = malloc(sizeof(uint64_t));
// *seed = 1000; // xorshift can't start with 0
// ^ not necessary, because parallel_sampler takes care of the seed.
int n_samples = 1000 * 1000 * 1000;
int n_threads = 16;
double sampler(uint64_t* seed){
return sample_lognormal(0, 10, seed);
}
double* results = malloc(n_samples * sizeof(double));
parallel_sampler(sampler, results, n_threads, n_samples);
double avg = array_sum(results, n_samples)/n_samples;
printf("Average of 1B lognormal(0,10): %f", avg);
free(results);
// free(seed);
// ^ not necessary, because parallel_sampler takes care of the seed.
}

BIN
squiggle.c/squiggle_c/examples/more/12_time_to_botec_parallel/example
View File

Binary file not shown.

29
squiggle.c/squiggle_c/examples/more/12_time_to_botec_parallel/example.c
View File

@ -1,29 +0,0 @@
#include "../../../squiggle.h"
#include "../../../squiggle_more.h"
#include <stdio.h>
#include <stdlib.h>
int main()
{
double p_a = 0.8;
double p_b = 0.5;
double p_c = p_a * p_b;
double sample_0(uint64_t* seed){ return 0; }
double sample_1(uint64_t* seed) { return 1; }
double sample_few(uint64_t* seed) { return sample_to(1, 3, seed); }
double sample_many(uint64_t* seed) { return sample_to(2, 10, seed); }
int n_dists = 4;
double weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
double (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many };
double sampler_result(uint64_t* seed) {
return sample_mixture(samplers, weights, n_dists, seed);
}
int n_samples = 1000 * 1000, n_threads = 16;
double* results = malloc(n_samples * sizeof(double));
parallel_sampler(sampler_result, results, n_threads, n_samples);
printf("Avg: %f\n", array_sum(results, n_samples)/n_samples);
free(results);
}

108
squiggle.c/squiggle_c/examples/more/makefile
View File

@ -1,108 +0,0 @@
# Interface:
# make all
# make format-all
# make run-all
# make one DIR=06_nuclear_recovery
# make format-one DIR=06_nuclear_recovery
# make run-one DIR=06_nuclear_recovery
# make time-linux-one DIR=06_nuclear_recovery
# make profile-one DIR=06_nuclear_recovery
# Compiler
CC=gcc
# CC=tcc # <= faster compilation
# Main file
SRC=example.c
OUTPUT=example
## Dependencies
SQUIGGLE=../../squiggle.c
SQUIGGLE_MORE=../../squiggle_more.c
MATH=-lm
OPENMP=-fopenmp
DEPS=$(SQUIGGLE) $(SQUIGGLE_MORE) $(MATH) $(OPENMP)
## Flags
DEBUG= #'-g'
STANDARD=-std=c99
WARNINGS=-Wall
OPTIMIZED=-O3 #-Ofast
## Formatter
STYLE_BLUEPRINT=webkit
FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT)
## make all
all:
$(CC) $(OPTIMIZED) $(DEBUG) 00_example_template/$(SRC) $(DEPS) -o 00_example_template/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 01_sample_from_cdf/$(SRC) $(DEPS) -o 01_sample_from_cdf/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 02_sample_from_cdf_beta/$(SRC) $(DEPS) -o 02_sample_from_cdf_beta/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 03_ci_beta/$(SRC) $(DEPS) -o 03_ci_beta/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 04_nuclear_war/$(SRC) $(DEPS) -o 04_nuclear_war/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 05_burn_10kg_fat/$(SRC) $(DEPS) -o 05_burn_10kg_fat/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 06_nuclear_recovery/$(SRC) $(DEPS) -o 06_nuclear_recovery/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 07_algebra/$(SRC) $(DEPS) -o 07_algebra/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 08_algebra_and_conversion/$(SRC) $(DEPS) -o 08_algebra_and_conversion/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 09_ergonomic_algebra/$(SRC) $(DEPS) -o 09_ergonomic_algebra/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 10_twitter_thread_example/$(SRC) $(DEPS) -o 10_twitter_thread_example/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 11_billion_lognormals_paralell/$(SRC) $(DEPS) -o 11_billion_lognormals_paralell/$(OUTPUT)
$(CC) $(OPTIMIZED) $(DEBUG) 12_time_to_botec_parallel/$(SRC) $(DEPS) -o 12_time_to_botec_parallel/$(OUTPUT)
format-all:
$(FORMATTER) 00_example_template/$(SRC)
$(FORMATTER) 01_sample_from_cdf/$(SRC)
$(FORMATTER) 02_sample_from_cdf_beta/$(SRC)
$(FORMATTER) 03_ci_beta/$(SRC)
$(FORMATTER) 04_nuclear_war/$(SRC)
$(FORMATTER) 05_burn_10kg_fat/$(SRC)
$(FORMATTER) 06_nuclear_recovery/$(SRC)
$(FORMATTER) 07_algebra/$(SRC)
$(FORMATTER) 08_algebra_and_conversion/$(SRC)
$(FORMATTER) 09_ergonomic_algebra/$(SRC)
$(FORMATTER) 10_twitter_thread_example/$(SRC)
$(FORMATTER) 11_billion_lognormals_paralell/$(SRC)
$(FORMATTER) 12_time_to_botec_parallel/$(SRC)
run-all:
00_example_template/$(OUTPUT)
01_sample_from_cdf/$(OUTPUT)
02_sample_from_cdf_beta/$(OUTPUT)
03_ci_beta/$(OUTPUT)
04_nuclear_war/$(OUTPUT)
05_burn_10kg_fat/$(OUTPUT)
06_nuclear_recovery/$(OUTPUT)
07_algebra/$(OUTPUT)
08_algebra_and_conversion/$(OUTPUT)
09_ergonomic_algebra/$(OUTPUT)
10_twitter_thread_example/$(OUTPUT)
11_billion_lognormals_paralell/$(OUTPUT)
12_time_to_botec_parallel/$(OUTPUT)
## make one DIR=06_nuclear_recovery
one: $(DIR)/$(SRC)
$(CC) $(OPTIMIZED) $(DEBUG) $(DIR)/$(SRC) $(DEPS) -o $(DIR)/$(OUTPUT)
## make format-one DIR=06_nuclear_recovery
format-one: $(DIR)/$(SRC)
$(FORMATTER) $(DIR)/$(SRC)
## make run-one DIR=06_nuclear_recovery
run-one: $(DIR)/$(OUTPUT)
$(DIR)/$(OUTPUT) && echo
## make time-linux-one DIR=06_nuclear_recovery
time-linux-one: $(DIR)/$(OUTPUT)
@echo "Requires /bin/time, found on GNU/Linux systems" && echo
@echo "Running 100x and taking avg time $(DIR)/$(OUTPUT)"
@t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do $(DIR)/$(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time: |" | sed 's|$$|ms|' && echo
## e.g., make profile-linux-one DIR=06_nuclear_recovery
profile-linux-one:
echo "Requires perf, which depends on the kernel version, and might be in linux-tools package or similar"
echo "Must be run as sudo"
$(CC) $(OPTIMIZED) $(DEBUG) $(DIR)/$(SRC) $(DEPS) -o $(DIR)/$(OUTPUT)
# $(CC) $(SRC) $(DEPS) -o $(OUTPUT)
sudo perf record $(DIR)/$(OUTPUT)
sudo perf report
rm perf.data

21
squiggle.c/squiggle_c/makefile
View File

@ -1,21 +0,0 @@
MAKEFLAGS += --no-print-directory
## Formatter
STYLE_BLUEPRINT=webkit
FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT)
build-examples:
cd examples/core && make all
cd examples/more && make all
format-examples:
cd examples/core && make format-all
cd examples/more && make format-all
format: squiggle.c squiggle.h
$(FORMATTER) squiggle.c
$(FORMATTER) squiggle.h
lint:
clang-tidy squiggle.c -- -lm
clang-tidy extra.c -- -lm

BIN
squiggle.c/squiggle_c/references/marsaglia-gamma.pdf
View File

Binary file not shown.

27
squiggle.c/squiggle_c/scratchpad/core.c
View File

@ -1,27 +0,0 @@
uint64_t xorshift64(uint64_t* seed)
{
// Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs"
// <https://en.wikipedia.org/wiki/Xorshift>
uint64_t x = *seed;
x ^= x << 13;
x ^= x >> 7;
x ^= x << 17;
return *seed = x;
}
double sample_unit_uniform(uint64_t* seed)
{
// samples uniform from [0,1] interval.
return ((double)xorshift64(seed)) / ((double)UINT64_MAX);
}
double sample_unit_normal(uint64_t* seed)
{
// // See: <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>
double u1 = sample_unit_uniform(seed);
double u2 = sample_unit_uniform(seed);
double z = sqrtf(-2.0 * log(u1)) * sin(2 * PI * u2);
return z;
}

56
squiggle.c/squiggle_c/scratchpad/makefile
View File

@ -1,56 +0,0 @@
# Interface:
# make
# make build
# make format
# make run
# Compiler
CC=gcc
# CC=tcc # <= faster compilation
# Main file
SRC=scratchpad.c ../squiggle.c
OUTPUT=scratchpad
## Dependencies
MATH=-lm
## Flags
DEBUG= #'-g'
STANDARD=-std=c99
WARNINGS=-Wall
OPTIMIZED=-O3 #-Ofast
# OPENMP=-fopenmp
## Formatter
STYLE_BLUEPRINT=webkit
FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT)
## make build
build: $(SRC)
$(CC) $(OPTIMIZED) $(DEBUG) $(SRC) $(MATH) -o $(OUTPUT)
format: $(SRC)
$(FORMATTER) $(SRC)
run: $(SRC) $(OUTPUT)
./$(OUTPUT)
verify: $(SRC) $(OUTPUT)
./$(OUTPUT) | grep "NOT passed" -A 2 --group-separator='' || true
time-linux:
@echo "Requires /bin/time, found on GNU/Linux systems" && echo
@echo "Running 100x and taking avg time $(OUTPUT)"
@t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do $(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 1 thread: |" | sed 's|$$|ms|' && echo
## Profiling
profile-linux:
echo "Requires perf, which depends on the kernel version, and might be in linux-tools package or similar"
echo "Must be run as sudo"
$(CC) $(SRC) $(MATH) -o $(OUTPUT)
sudo perf record ./$(OUTPUT)
sudo perf report
rm perf.data

7
squiggle.c/squiggle_c/scratchpad/plotting/.Rhistory
View File

@ -1,7 +0,0 @@
library(ggplot2)
data <- read.csv("samples.txt", header = FALSE)
data <- as.data.frame(data)
ggplot(data = data, aes(x = V1)) +
geom_bar()
ggplot(data = data, aes(x = V1)) +
geom_freqpoly()

3
squiggle.c/squiggle_c/scratchpad/plotting/makefile
View File

@ -1,3 +0,0 @@
build:
format:

12
squiggle.c/squiggle_c/scratchpad/plotting/notes.md
View File

@ -1,12 +0,0 @@
options:
- use gnuplot
- requires aggregation
- do aggregation in its own gnuplot language
- use something like awk
- use some other language, like python or R, specifically to do the aggregation
- do the aggregation in C
- use ggplot
- requires setting up bins correctly
- also can't pipe to an r script directly, sadly enough.
- use python for plotting?
- just present the confidence intervals in C?

0
squiggle.c/squiggle_c/scratchpad/plotting/r/.Rhistory
View File

7
squiggle.c/squiggle_c/scratchpad/plotting/r/plot.R
View File

@ -1,7 +0,0 @@
library(ggplot2)
data <- read.csv("samples.txt", header = FALSE)
data <- as.data.frame(data)
ggplot(data = data, aes(x = V1)) +
geom_freqpoly()

1000000
squiggle.c/squiggle_c/scratchpad/plotting/r/samples.txt
View File

File diff suppressed because it is too large Load Diff

BIN
squiggle.c/squiggle_c/scratchpad/plotting/src/example
View File

Binary file not shown.

50
squiggle.c/squiggle_c/scratchpad/plotting/src/example.c
View File

@ -1,50 +0,0 @@
#include "../../squiggle.h"
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
// Estimate functions
double sample_0(uint64_t* seed)
{
return 0;
}
double sample_1(uint64_t* seed)
{
return 1;
}
double sample_few(uint64_t* seed)
{
return sample_to(1, 3, seed);
}
double sample_many(uint64_t* seed)
{
return sample_to(2, 10, seed);
}
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
double p_a = 0.8;
double p_b = 0.5;
double p_c = p_a * p_b;
int n_dists = 4;
double weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
double (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many };
int n_samples = 1000000;
double* result_many = (double*)malloc(n_samples * sizeof(double));
for (int i = 0; i < n_samples; i++) {
result_many[i] = sample_mixture(samplers, weights, n_dists, seed);
printf("%f\n", result_many[i]);
}
// printf("Mean: %f\n", array_mean(result_many, n_samples));
free(seed);
}

53
squiggle.c/squiggle_c/scratchpad/plotting/src/makefile
View File

@ -1,53 +0,0 @@
# Interface:
# make
# make build
# make format
# make run
# Compiler
CC=gcc
# CC=tcc # <= faster compilation
# Main file
SRC=example.c ../../squiggle.c
OUTPUT=example
## Dependencies
MATH=-lm
## Flags
DEBUG= #'-g'
STANDARD=-std=c99
WARNINGS=-Wall
OPTIMIZED=-O3 #-Ofast
# OPENMP=-fopenmp
## Formatter
STYLE_BLUEPRINT=webkit
FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT)
## make build
build: $(SRC)
$(CC) $(OPTIMIZED) $(DEBUG) $(SRC) $(MATH) -o $(OUTPUT)
format: $(SRC)
$(FORMATTER) $(SRC)
run: $(SRC) $(OUTPUT)
./$(OUTPUT) && echo
time-linux:
@echo "Requires /bin/time, found on GNU/Linux systems" && echo
@echo "Running 100x and taking avg time $(OUTPUT)"
@t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do ./$(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 1 thread: |" | sed 's|$$|ms|' && echo
## Profiling
profile-linux:
echo "Requires perf, which depends on the kernel version, and might be in linux-tools package or similar"
echo "Must be run as sudo"
$(CC) $(SRC) $(MATH) -o $(OUTPUT)
sudo perf record ./$(OUTPUT)
sudo perf report
rm perf.data

BIN
squiggle.c/squiggle_c/scratchpad/scratchpad
View File

Binary file not shown.

22
squiggle.c/squiggle_c/scratchpad/scratchpad.c
View File

@ -1,22 +0,0 @@
#include "../squiggle.h"
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with a seed of 0
/*
for (int i = 0; i < 100; i++) {
double draw = sample_unit_uniform(seed);
printf("%f\n", draw);
}*/
// Test division
printf("\n%d\n", 10 % 3);
free(seed);
}

56
squiggle.c/squiggle_c/test/makefile
View File

@ -1,56 +0,0 @@
# Interface:
# make
# make build
# make format
# make run
# Compiler
CC=gcc
# CC=tcc # <= faster compilation
# Main file
SRC=test.c ../squiggle.c
OUTPUT=test
## Dependencies
MATH=-lm
## Flags
DEBUG= #'-g'
STANDARD=-std=c99
WARNINGS=-Wall
OPTIMIZED=-O3 #-Ofast
# OPENMP=-fopenmp
## Formatter
STYLE_BLUEPRINT=webkit
FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT)
## make build
build: $(SRC)
$(CC) $(OPTIMIZED) $(DEBUG) $(SRC) $(MATH) -o $(OUTPUT)
format: $(SRC)
$(FORMATTER) $(SRC)
run: $(SRC) $(OUTPUT)
./$(OUTPUT)
verify: $(SRC) $(OUTPUT)
./$(OUTPUT) | grep "NOT passed" -A 2 --group-separator='' || true
time-linux:
@echo "Requires /bin/time, found on GNU/Linux systems" && echo
@echo "Running 100x and taking avg time $(OUTPUT)"
@t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do $(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 1 thread: |" | sed 's|$$|ms|' && echo
## Profiling
profile-linux:
echo "Requires perf, which depends on the kernel version, and might be in linux-tools package or similar"
echo "Must be run as sudo"
$(CC) $(SRC) $(MATH) -o $(OUTPUT)
sudo perf record ./$(OUTPUT)
sudo perf report
rm perf.data

BIN
squiggle.c/squiggle_c/test/test
View File

Binary file not shown.

328
squiggle.c/squiggle_c/test/test.c
View File

@ -1,328 +0,0 @@
#include "../squiggle.h"
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#define TOLERANCE 5.0 / 1000.0
#define MAX_NAME_LENGTH 500
// Structs
struct array_expectations {
double* array;
int n;
char* name;
double expected_mean;
double expected_std;
double tolerance;
};
void test_array_expectations(struct array_expectations e)
{
double mean = array_mean(e.array, e.n);
double delta_mean = mean - e.expected_mean;
double std = array_std(e.array, e.n);
double delta_std = std - e.expected_std;
if ((fabs(delta_mean) / fabs(mean) > e.tolerance) && (fabs(delta_mean) > e.tolerance)) {
printf("[-] Mean test for %s NOT passed.\n", e.name);
printf("Mean of %s: %f, vs expected mean: %f\n", e.name, mean, e.expected_mean);
printf("delta: %f, relative delta: %f\n", delta_mean, delta_mean / fabs(mean));
} else {
printf("[x] Mean test for %s PASSED\n", e.name);
}
if ((fabs(delta_std) / fabs(std) > e.tolerance) && (fabs(delta_std) > e.tolerance)) {
printf("[-] Std test for %s NOT passed.\n", e.name);
printf("Std of %s: %f, vs expected std: %f\n", e.name, std, e.expected_std);
printf("delta: %f, relative delta: %f\n", delta_std, delta_std / fabs(std));
} else {
printf("[x] Std test for %s PASSED\n", e.name);
}
printf("\n");
}
// Test unit uniform
void test_unit_uniform(uint64_t* seed)
{
int n = 1000 * 1000;
double* unit_uniform_array = malloc(sizeof(double) * n);
for (int i = 0; i < n; i++) {
unit_uniform_array[i] = sample_unit_uniform(seed);
}
struct array_expectations expectations = {
.array = unit_uniform_array,
.n = n,
.name = "unit uniform",
.expected_mean = 0.5,
.expected_std = sqrt(1.0 / 12.0),
.tolerance = TOLERANCE,
};
test_array_expectations(expectations);
free(unit_uniform_array);
}
// Test uniforms
void test_uniform(double start, double end, uint64_t* seed)
{
int n = 1000 * 1000;
double* uniform_array = malloc(sizeof(double) * n);
for (int i = 0; i < n; i++) {
uniform_array[i] = sample_uniform(start, end, seed);
}
char* name = malloc(MAX_NAME_LENGTH * sizeof(char));
snprintf(name, MAX_NAME_LENGTH, "[%f, %f] uniform", start, end);
struct array_expectations expectations = {
.array = uniform_array,
.n = n,
.name = name,
.expected_mean = (start + end) / 2,
.expected_std = sqrt(1.0 / 12.0) * fabs(end - start),
.tolerance = fabs(end - start) * TOLERANCE,
};
test_array_expectations(expectations);
free(name);
free(uniform_array);
}
// Test unit normal
void test_unit_normal(uint64_t* seed)
{
int n = 1000 * 1000;
double* unit_normal_array = malloc(sizeof(double) * n);
for (int i = 0; i < n; i++) {
unit_normal_array[i] = sample_unit_normal(seed);
}
struct array_expectations expectations = {
.array = unit_normal_array,
.n = n,
.name = "unit normal",
.expected_mean = 0,
.expected_std = 1,
.tolerance = TOLERANCE,
};
test_array_expectations(expectations);
free(unit_normal_array);
}
// Test normal
void test_normal(double mean, double std, uint64_t* seed)
{
int n = 10 * 1000 * 1000;
double* normal_array = malloc(sizeof(double) * n);
for (int i = 0; i < n; i++) {
normal_array[i] = sample_normal(mean, std, seed);
}
char* name = malloc(MAX_NAME_LENGTH * sizeof(char));
snprintf(name, MAX_NAME_LENGTH, "normal(%f, %f)", mean, std);
struct array_expectations expectations = {
.array = normal_array,
.n = n,
.name = name,
.expected_mean = mean,
.expected_std = std,
.tolerance = TOLERANCE,
};
test_array_expectations(expectations);
free(name);
free(normal_array);
}
// Test lognormal
void test_lognormal(double logmean, double logstd, uint64_t* seed)
{
int n = 10 * 1000 * 1000;
double* lognormal_array = malloc(sizeof(double) * n);
for (int i = 0; i < n; i++) {
lognormal_array[i] = sample_lognormal(logmean, logstd, seed);
}
char* name = malloc(MAX_NAME_LENGTH * sizeof(char));
snprintf(name, MAX_NAME_LENGTH, "lognormal(%f, %f)", logmean, logstd);
struct array_expectations expectations = {
.array = lognormal_array,
.n = n,
.name = name,
.expected_mean = exp(logmean + pow(logstd, 2) / 2),
.expected_std = sqrt((exp(pow(logstd, 2)) - 1) * exp(2 * logmean + pow(logstd, 2))),
.tolerance = TOLERANCE,
};
test_array_expectations(expectations);
free(name);
free(lognormal_array);
}
// Test lognormal to
void test_to(double low, double high, uint64_t* seed)
{
int n = 10 * 1000 * 1000;
double* lognormal_array = malloc(sizeof(double) * n);
for (int i = 0; i < n; i++) {
lognormal_array[i] = sample_to(low, high, seed);
}
char* name = malloc(MAX_NAME_LENGTH * sizeof(char));
snprintf(name, MAX_NAME_LENGTH, "to(%f, %f)", low, high);
const double NORMAL95CONFIDENCE = 1.6448536269514722;
double loglow = logf(low);
double loghigh = logf(high);
double logmean = (loglow + loghigh) / 2;
double logstd = (loghigh - loglow) / (2.0 * NORMAL95CONFIDENCE);
struct array_expectations expectations = {
.array = lognormal_array,
.n = n,
.name = name,
.expected_mean = exp(logmean + pow(logstd, 2) / 2),
.expected_std = sqrt((exp(pow(logstd, 2)) - 1) * exp(2 * logmean + pow(logstd, 2))),
.tolerance = TOLERANCE,
};
test_array_expectations(expectations);
free(name);
free(lognormal_array);
}
// Test beta
void test_beta(double a, double b, uint64_t* seed)
{
int n = 10 * 1000 * 1000;
double* beta_array = malloc(sizeof(double) * n);
for (int i = 0; i < n; i++) {
beta_array[i] = sample_beta(a, b, seed);
}
char* name = malloc(MAX_NAME_LENGTH * sizeof(char));
snprintf(name, MAX_NAME_LENGTH, "beta(%f, %f)", a, b);
struct array_expectations expectations = {
.array = beta_array,
.n = n,
.name = name,
.expected_mean = a / (a + b),
.expected_std = sqrt((a * b) / (pow(a + b, 2) * (a + b + 1))),
.tolerance = TOLERANCE,
};
test_array_expectations(expectations);
free(name);
}
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with a seed of 0
printf("Testing unit uniform\n");
test_unit_uniform(seed);
printf("Testing small uniforms\n");
for (int i = 0; i < 100; i++) {
double start = sample_uniform(-10, 10, seed);
double end = sample_uniform(-10, 10, seed);
if (end > start) {
test_uniform(start, end, seed);
}
}
printf("Testing wide uniforms\n");
for (int i = 0; i < 100; i++) {
double start = sample_uniform(-1000 * 1000, 1000 * 1000, seed);
double end = sample_uniform(-1000 * 1000, 1000 * 1000, seed);
if (end > start) {
test_uniform(start, end, seed);
}
}
printf("Testing unit normal\n");
test_unit_normal(seed);
printf("Testing small normals\n");
for (int i = 0; i < 100; i++) {
double mean = sample_uniform(-10, 10, seed);
double std = sample_uniform(0, 10, seed);
if (std > 0) {
test_normal(mean, std, seed);
}
}
printf("Testing larger normals\n");
for (int i = 0; i < 100; i++) {
double mean = sample_uniform(-1000 * 1000, 1000 * 1000, seed);
double std = sample_uniform(0, 1000 * 1000, seed);
if (std > 0) {
test_normal(mean, std, seed);
}
}
printf("Testing smaller lognormals\n");
for (int i = 0; i < 100; i++) {
double mean = sample_uniform(-1, 1, seed);
double std = sample_uniform(0, 1, seed);
if (std > 0) {
test_lognormal(mean, std, seed);
}
}
printf("Testing larger lognormals\n");
for (int i = 0; i < 100; i++) {
double mean = sample_uniform(-1, 5, seed);
double std = sample_uniform(0, 5, seed);
if (std > 0) {
test_lognormal(mean, std, seed);
}
}
printf("Testing lognormals — sample_to(low, high) syntax\n");
for (int i = 0; i < 100; i++) {
double low = sample_uniform(0, 1000 * 1000, seed);
double high = sample_uniform(0, 1000 * 1000, seed);
if (low < high) {
test_to(low, high, seed);
}
}
// Bonus example
test_to(10, 10 * 1000, seed);
printf("Testing beta distribution\n");
for (int i = 0; i < 100; i++) {
double a = sample_uniform(0, 1000, seed);
double b = sample_uniform(0, 1000, seed);
if ((a > 0) && (b > 0)) {
test_beta(a, b, seed);
}
}
printf("Testing larger beta distributions\n");
for (int i = 0; i < 100; i++) {
double a = sample_uniform(0, 1000 * 1000, seed);
double b = sample_uniform(0, 1000 * 1000, seed);
if ((a > 0) && (b > 0)) {
test_beta(a, b, seed);
}
}
free(seed);
}