diff --git a/lilypond/attempt2/genlp.py b/lilypond/attempt2/genlp.py index 8470ef4..7ba0b33 100755 --- a/lilypond/attempt2/genlp.py +++ b/lilypond/attempt2/genlp.py @@ -20,8 +20,6 @@ def __init__(self, x, y, t, v): self.y = y self.t = t self.v = v - self.r = r - self.i = i def __repr__(self): return "%f,%f,%f,%f,0,0" % (self.x, self.y, self.t, self.v) @@ -33,7 +31,7 @@ def main(): lambda t, y: dimX[0], lambda t, y: dimX[1]) # '_' assigned numerical error # single grain simulation is boring - N = np.random.poisson(maxLam) + N = np.random.poisson(lam) generate_grains(N, dimX, dimY, dimT, lamDensity, lam) def generate_grains(N, dimX, dimY, dimT, lamDensity, maxLam = None): @@ -44,9 +42,9 @@ def generate_grains(N, dimX, dimY, dimT, lamDensity, maxLam = None): n = 0 while n < N: potentialGrain = (np.random.uniform(*dimX), np.random.uniform(*dimY), np.random.uniform(*dimT)) - if np.random.uniform(0,1) <= lamd(*potentialGrain)/maxlam: + if np.random.uniform(0,1) <= lamDensity(*potentialGrain)/maxLam: # accept point - print(Grain(*tmp, np.random.uniform()/100)) + print(Grain(*potentialGrain, np.random.uniform()/100)) n += 1 if __name__ == "__main__": diff --git a/lilypond/attempt2/grain.c b/lilypond/attempt2/grain.c new file mode 100644 index 0000000..f0373ac --- /dev/null +++ b/lilypond/attempt2/grain.c @@ -0,0 +1,355 @@ +/** + * I don't like how coupled this is all getting... + */ + +#include "grain.h" + +/** + * Contains all information pertaining to a single grain + * x,y are spatial coordinates + * t is the birth time + * v is growth rate + * r is maximum radius. Set to -1 if unknown and 0 if suffocated. + * i is the Grain it collides with if r is positive + * lst is a linkedlist of all pairs it's a part of + * spair is a special pair since it's the best pair involving a stationary + */ +struct Grain { + size_t id; + double x, y, t; + double v; + double r; + ssize_t i; + struct GrainPairListNode *lst; + struct GrainPair *spair; +}; + +/** + * Contains all information about a collision between two grains + * score contains the cached value of the time/distance/general metric + * between two grains. Should be updated every timestep. + * gnode is the global node attached to this GrainPair. + */ +struct GrainPair { + struct Grain* g1; + struct Grain* g2; + double score; + struct GrainPairListNode *glst; +} + +/** + * Linkedlist containing all the grain pairs. Very useful. + * l and r are left and right + */ +struct GrainPairListNode { + struct GrainPair *gp; + struct GrainPairListNode *l; + struct GrainPairListNode *r; +}; + +/** + * Performs free_gpln on all nodes to the left of lst + */ +void free_gpln_l(struct GrainPairListNode *lst) { + if (lst) { /* NULL */ + return; + } + + struct GrainPairListNode *l = lst->l; + + free_gpln(lst); + + free_gpln_l(l); +} + +/** + * Performs free_gpln on all nodes to the right of lst + */ +void free_gpln_r(struct GrainPairListNode *lst) { + if (lst) { /* NULL */ + return; + } + + struct GrainPairListNode *r = lst->r; + + free_gpln(lst); + + free_gpln_r(r); +} + +/** + * Creates (calloc) space for a new GrainPair + * Returns a new GrainPair on success + * NULL on failure + */ +struct Grain *new_g(size_t id, double x, double y, double t, + double v) { + struct Grain *ret = calloc(1, sizeof(ret)); + + ret->id = id; + ret->x = x; + ret->y = y; + ret->t = t; + ret->v = v; + + ret->r = -1; + ret->i = -1; + ret->lst = NULL; + + return ret; +} + +/** + * Frees and cleans up the Grain. + * Will also free all GrainPairs ensuring they properly + * remove their GrainPairListNodes from their parent lists. + */ +void free_g(struct Grain *g) { + for (struct GrainPairListNode *node = g->lst; + free(g); +} + +/** + * Goes through its list and removes the pairs with two stationary in them. + * + */ +void set_g(struct Grain *g, double r, size_t i) { + +} + +size_t get_g_id(struct Grain *g) { + +} + +double get_g_x(struct Grain *g) { + +} + +double get_g_y(struct Grain *g) { + +} + +double get_g_t(struct Grain *g) { + +} + +double get_g_v(struct Grain *g) { + +} + +/** + * r is maximum radius. Set to -1 if unknown and 0 if suffocated. + */ +double get_g_r(struct Grain *g) { + +} + +/** + * i is the Grain it collides with if r is positive + * otherwise return -1 + */ +ssize_t get_g_i(struct Grain *g) { + +} + +/** + * Creates (calloc) space for a new GrainPair + * Returns a new GrainPair on success + * NULL on failure + * + * [IMPL] This function should append itself to each grains linked + * list. We can do this because these pairs are expected to be unique + * unlike the GrainPairListNodes + */ +struct GrainPair *new_gp(struct Grain *g1, struct Grain *g2) { + +} + +/** + * Frees and cleans up the GrainPair. + * Will also free all GrainPairListNodes ensuring they're properly + * removed from their lists. + * Doesn't get rid of the Grains. + */ +void free_gp(struct GrainPair *gp) { + +} + +/** + * Gets the first grain of the grain pair + */ +double get_gp_1(struct GrainPair *gp) { + +} + +/** + * Gets the second grain of the grain pair + */ +double get_gp_2(struct GrainPair *gp) { + +} + +/** + * Checks if both grains in the pair are already realised + */ +bool is_gp_done(struct GrainPair *gp) { + +} + +/** + * Sets the glst to what the GrainPair considers is the GlobalListNode it resides in + */ +void set_gp_glst(struct Grainpair *gp, struct GrainPairListNode *glst) { + +} + +/** + * Computes the score of the GrainPair but doesn't store it + */ +double calc_gp_score(struct GrainPair *gp) { + +} + +/** + * Gets the score of the grain pair + */ +double get_gp_score(struct GrainPair *gp) { + +} + +/** + * Sets the score of the grain pair + */ +void set_gp_score(struct GrainPair *gp, double newscore) { + +} + +/** + * Create (calloc) space for a new GrainPairListNode + * Returns a new GrainPairListNode on success + * NULL on failure + */ +struct GrainPairListNode *new_gpln(struct GrainPair *gp) { + +} + +/** + * Frees and cleans up a node. + * Doesn't get rid of the GrainPair. + */ +void free_gpln(struct GrainPairListNode *node) { + if (!node) { /* NULL */ + return NULL; + } + + free(node); +} + +/** + * Gets the left node of the GrainPairListNode + * Returns NULL if there isn't one + */ +struct GrainPairListNode *get_gpln_l(struct GrainPairListNode *node) { + if (!node) { /* NULL */ + return NULL; + } + + return node->r; +} + +/** + * Gets the right node of the GrainPairListNode + * Returns NULL if there isn't one + */ +struct GrainPairListNode *get_gpln_r(struct GrainPairListNode *node) { + if (!node) { /* NULL */ + return NULL; + } + + return node->r; +} + +/** + * Gets the GrainPair from the GrainPairListNode + * Returns NULL if there isn't one + */ +struct GrainPair *get_gpln_gp(struct GrainPairListNode *node) { + if (!node) { /* NULL */ + return NULL; + } + + return node->gp; +} + +/** + * Inserts the node to the left of the list + */ +void add_gpln_l(struct GrainPairListNode *lst, struct GrainPairListNode *node) { + if (!(lst && node)) { /* NULL */ + return; + } + + struct GrainPairListNode oldl = get_gpln_l(lst); + + lst->l = node; + node->r = lst; + + if (oldl) { /* NOTNULL */ + oldl->r = node; + node->l = oldl; + } +} + +/** + * Inserts the node to the right of the list + */ +void add_gpln_r(struct GrainPairListNode *lst, struct GrainPairListNode *node) { + if (!(lst && node)) { /* NULL */ + return; + } + + struct GrainPairListNode oldr = get_gpln_r(lst); + + lst->r = node; + node->l = lst; + + if (oldr) { /* NOTNULL */ + oldr->l = node; + node->r = oldr; + } +} + +/** + * Removes the current GrainPairListNode from its list + */ +void rm_gpln(struct GrainPairListNode *node) { + if (!node) { /* NULL */ + return; + } + + struct GrainPairListNode *l = get_gpln_l(node); + struct GrainPairListNode *r = get_gpln_r(node); + + node->l = NULL; + node->r = NULL; + + if (l) { /* NOTNULL */ + l->r = r; + } + if (r) { /* NOTNULL */ + r->l = l; + } +} + +/** + * Performs free_gpln on all nodes of lst including self + */ +void free_gpln_a(struct GrainPairListNode *lst) { + if (!lst) { /* NULL */ + return; + } + + free_gpln_l(lst->l); + free_gpln_r(lst->r); + free_gpln(lst); +} diff --git a/lilypond/attempt2/grain.h b/lilypond/attempt2/grain.h new file mode 100644 index 0000000..ed8afe2 --- /dev/null +++ b/lilypond/attempt2/grain.h @@ -0,0 +1,160 @@ +#ifndef GRAIN_H +#define GRAIN_H + +#include +#include + +/* Definitions in grain.c */ +struct Grain; +struct GrainPair; +struct GrainPairListNode; + +/** + * Creates (calloc) space for a new GrainPair + * Returns a new GrainPair on success + * NULL on failure + */ +struct Grain *new_g(size_t id, double x, double y, double t, + double v); +/** + * Frees and cleans up the Grain. + * Will also free all GrainPairs ensuring they properly + * remove their GrainPairListNodes from their parent lists. + */ +void free_g(struct Grain *g); + +/** + * Once we've solved for r and i we do a lot of behind-the-scenes + * updating. + * + * [IMPL] This function actually does a LOT of updating. + * Removes the spair + * Searches through its lst in search of other grains + * Removes it's own entry from other grains and replaces spair + * if it's better otherwise nah. + */ +void set_g(struct Grain *g, double r, size_t i); + +size_t get_g_id(struct Grain *g); +double get_g_x(struct Grain *g); +double get_g_y(struct Grain *g); +double get_g_t(struct Grain *g); +double get_g_v(struct Grain *g); + +/** + * r is maximum radius. Is -1 if unknown and 0 if suffocated. + */ +double get_g_r(struct Grain *g); + +/** + * i is the Grain it collides with if r is positive + * otherwise return -1 + */ +ssize_t get_g_i(struct Grain *g); + +/** + * Creates (calloc) space for a new GrainPair + * Returns a new GrainPair on success + * NULL on failure + * + * [IMPL] This function should append itself to each grains linked + * list. We can do this because these pairs are expected to be unique + * unlike the GrainPairListNodes + */ +struct GrainPair *new_gp(struct Grain *g1, struct Grain *g2); + +/** + * Frees and cleans up the GrainPair. + * Will also free all GrainPairListNodes ensuring they're properly + * removed from their lists. + * Doesn't get rid of the Grains. + */ +void free_gp(struct GrainPair *gp); + +/** + * Gets the first grain of the grain pair + */ +double get_gp_1(struct GrainPair *gp); + +/** + * Gets the second grain of the grain pair + */ +double get_gp_2(struct GrainPair *gp); + +/** + * Checks if both grains in the pair are already realised + */ +bool is_gp_done(struct GrainPair *gp); + +/** + * Sets the glst to what the GrainPair considers is the GlobalListNode it resides in + */ +void set_gp_glst(struct Grainpair *gp, struct GrainPairListNode *glst); + +/** + * Computes the score of the GrainPair but doesn't store it + */ +double calc_gp_score(struct GrainPair *gp); + +/** + * Gets the score of the grain pair + */ +double get_gp_score(struct GrainPair *gp); + +/** + * Sets the score of the grain pair + */ +void set_gp_score(struct GrainPair *gp, double newscore); + +/** + * Create (calloc) space for a new GrainPairListNode + * Returns a new GrainPairListNode on success + * NULL on failure + */ +struct GrainPairListNode *new_gpln(struct GrainPair *gp); + +/** + * Frees and cleans up a node. + * Doesn't get rid of the GrainPair. + */ +void free_gpln(struct GrainPairListNode *node); + +/** + * Gets the left node of the GrainPairListNode + * Returns NULL if there isn't one + */ +struct GrainPairListNode *get_gpln_l(struct GrainPairListNode *node); + +/** + * Gets the right node of the GrainPairListNode + * Returns NULL if there isn't one + */ +struct GrainPairListNode *get_gpln_r(struct GrainPairListNode *node); + +/** + * Gets the GrainPair from the GrainPairListNode + * Returns NULL if there isn't one + */ +struct GrainPair *get_gpln_gp(struct GrainPairListNode *node); + +/** + * Inserts the node to the left of the list + */ +void add_gpln_l(struct GrainPairListNode *lst, struct GrainPairListNode *node); + +/** + * Inserts the node to the right of the list + */ +void add_gpln_r(struct GrainPairListNode *lst, struct GrainPairListNode *node); + +/** + * Removes the current GrainPairListNode from its list + */ +void rm_gpln(struct GrainPairListNode *node); + +/** + * Performs free_gpln on all nodes of lst including self + */ +void free_gpln_a(struct GrainPairListNode *lst); + +#endif diff --git a/lilypond/attempt2/solvelp b/lilypond/attempt2/solvelp index 6f3623b..b21572f 100755 Binary files a/lilypond/attempt2/solvelp and b/lilypond/attempt2/solvelp differ diff --git a/lilypond/attempt2/solvelp.c b/lilypond/attempt2/solvelp.c index 62fde6a..c1a41cd 100644 --- a/lilypond/attempt2/solvelp.c +++ b/lilypond/attempt2/solvelp.c @@ -1,4 +1,5 @@ #include +#include #include #include #include @@ -7,9 +8,7 @@ #include #include -/* Using the gnu science libraries */ -#include -#include +#include "grain.h" void print_usage(const char *progname) { fprintf(stderr, "Usage: %s [-i inputfile] [-o outputfile]\n", progname); @@ -74,7 +73,6 @@ void parse_args(int argc, char **argv) { int main(int argc, char **argv) { parse_args(argc, argv); - printf("hello\n"); return 0; } diff --git a/lilypond/attempt2/visualiser.py b/lilypond/attempt2/visualiser.py index 3c0022e..5235991 100755 --- a/lilypond/attempt2/visualiser.py +++ b/lilypond/attempt2/visualiser.py @@ -3,9 +3,15 @@ import csv import math +import argparse +import fileinput + from matplotlib import pyplot as plt from matplotlib import animation +parser = argparse.ArgumentParser(description="Takes an input of *.lp format and generates a visual") +parser.add_argument("inputfile", help="if specified reads a *.lp formatted file otherwise standard in") + class Circle: def __init__(self, ID, x, y, t, v, r, i): self.id = ID @@ -17,38 +23,33 @@ def __init__(self, ID, x, y, t, v, r, i): self.shape = plt.Circle(self.xy, 0) def __repr__(self): - return "Circle:{id: %d, x: %f, y: %f, t: %f, v: %f, r: %f, i: %d}" % (self.id, *self.xy, self.t, self.v, self.r, self.i) + return "Circle:{id: %d, x: %f, y: %f, t: %f, v: %f, r: %f, i: %d}" \ + % (self.id, *self.xy, self.t, self.v, self.r, self.i) def main(): - if len(sys.argv) != 2: - print("Usage: %s filename" % sys.argv[0]) - sys.exit(1) - - fileName = sys.argv[1] - with open(fileName) as csvFile: - reader = csv.reader(row for row in csvFile if not row.startswith("#")) - - circles = [] - dimX, dimY, dimT = (math.inf, -math.inf), (math.inf, -math.inf), (0, 0) - n = 0 # this is the index of the next circle - try: # parsing - for row in reader: - x, y, t, v, r, i = row # see example.lp for these values - # perform a very forgiving parse of the file - x, y, t, v, r, i = float(x), float(y), abs(float(t)), abs(float(v)), abs(float(r)), abs(int(i)) - dimX = (min(dimX[0], x), max(dimX[1], x)) - dimY = (min(dimY[0], y), max(dimY[1], y)) - dimT = (0, math.ceil((max(dimT[1], t + r/v)))) # this is determined by the birth and how long it takes to reach max radius - circles.append(Circle(n, x, y, t, v, r, i)) - n += 1 - except ValueError as e: - raise ValueError("Parsing error has occurred in %s at %s" % (fileName, row[0])) - - # check that all i < n - for circle in (circle for circle in circles if circle.i >= n): - raise IndexError("circle %d collides with circle %d which does not exist" % (circle.id, circle.i)) - - visualise(circles, dimX, dimY, dimT) + reader = csv.reader(row for row in fileinput.input() if not row.startswith("#")) + + circles = [] + dimX, dimY, dimT = (math.inf, -math.inf), (math.inf, -math.inf), (0, 0) + n = 0 # this is the index of the next circle + try: # parsing + for row in reader: + x, y, t, v, r, i = row # see example.lp for these values + # perform a very forgiving parse of the file + x, y, t, v, r, i = float(x), float(y), abs(float(t)), abs(float(v)), abs(float(r)), abs(int(i)) + dimX = (min(dimX[0], x), max(dimX[1], x)) + dimY = (min(dimY[0], y), max(dimY[1], y)) + dimT = (0, math.ceil((max(dimT[1], t + r/v)))) # this is determined by the birth and how long it takes to reach max radius + circles.append(Circle(n, x, y, t, v, r, i)) + n += 1 + except ValueError as e: + raise ValueError("Parsing error has occurred in %s at %s" % (fileinput.filename(), row[0])) + + # check that all i < n + for circle in (circle for circle in circles if circle.i >= n): + raise IndexError("circle %d collides with circle %d which does not exist" % (circle.id, circle.i)) + + visualise(circles, dimX, dimY, dimT) # displays the list of circles as a matplotlib figure diff --git a/lilypond/journal/2019-11-29.md b/lilypond/journal/2019-11-29.md new file mode 100644 index 0000000..e89ce1b --- /dev/null +++ b/lilypond/journal/2019-11-29.md @@ -0,0 +1,65 @@ +# Journal entry for 29 November 2019 +## What have I done so far. +* I've implemented the Poisson Lilypond model in Python but it's currently very slow +and uses numerical integration for the lambda. +There's a single remaining bug but for most cases it works just fine. + * The model has been extended so that various growth rates are possible and various birth times too. + * Implemented inhomogeneous intensity through thinning + +* I've read some of the literature on the model. + * I'm overwhelmed with how dense these papers are. + The proofs are going well over my heads but I'm understanding some of + the statements they're making. + * Percolation and the Limit Theory for the Poisson Lilypond Model (Last & Penrose, 2013) + * Dissected as much as I could + * On the absence of percolation in a line-segment based lilypond model (Hirsch, 2014) + * Have not read much of this paper + +* I've also gotten some books on percolation which may or may not be useful. + * Continuum percolation (Mester & Roy, 1996) + * Had loads of information on Poisson Boolean models. + * Percolation (Grimmett, 1999) + * Had mostly discrete models don't think there's anything useful there. + +## What do I intend to do the coming week. +* Rewrite the program in C for a small boost in speed. Also document it well. + +* I would like to visualise the graph made from the collisions of grains. + +* Also would like to be able to compare the models. A technique called coupling + you force two independent models to become dependent so you may easily compare them. + Run the simulation on the points normally then one where you ignore birth time etc. + +* Eventually figure out what I want to analyse. + +## Long term +* Look into extending it to 3D (which shouldn't be that hard) + +* Alternative shapes and distance metrics (manhattan distance) + +## Lingering questions. +* I'm currently unsure of where to take the project next. + * Trying to prove something new about the model seems mighty difficult. + * Reduce the complexity of the model and focus on a small part. That's what + Last and Penrose (2010) did they percolated with a possible jumping distance. + Doing so gave them a critical jumping distance which almost always percolates. + * Focus on random birth times + * Ebert and Last (2015) already proved the model exists and + is unique and does not percolate and a central limit. + They also showed how it relates to Apollonian packing. + * Focus on varying growth rates + * How should they be distributed? What would be natural? + * Focus on non-homogeneous intensity + * What lambda density should be used? + * Exploring clustering + * What would I even prove? + * Percolation + * Largest grain size + * Largest connected component (by what metric: volume, number of grains, dimensions) + * Density/saturation coverage/packing fraction + * Proving a central limit theorem + * Finding ways to improve the algorithm because it's currently running at n^4 time + +* An issue with the current model is that it's bounded by the dimensions we define. + In reality this model should have an infinite number of circles. Is there any way + we can approximate this behaviour?