Tuesday, February 22, 2011

Adaptivity Tasks

I will use this to keep track of what parts of the adaptivity code I have finished so far.

Refinement

Implementation
Set Initial Marks: Complete
Bisect Tetrahedron: Complete
Recursive Propagation: Complete
Progress Constraints: Not done; Requires theoretical development

Testing
Basic Testing: Complete
Physics Stub: Not done; Requires code from physics
Real Physics: Not done; Requires code from physics


Edge Contraction

Implementation
Build Patch:

-Check for Legality
--Progress Constraints/New Markings: Not done; Requires theoretical development
--Inversion: Complete
--Mesh Quality: Not done; Requires theoretical development
-Create New Patch: Complete

Update Front:
---Adjacency/Predecessor information: Complete
---Change Markings: Not done

Testing
Basic Testing: In progress
Physics Stub: Not done; Requires code from physics
Real Physics: Not done; Requires code from physics

Visualization
Possible for future; Not designed yet

Thursday, January 27, 2011

Qual Practice - Fall 2009 Part I (first two problems)

I am preparing for the qual now so I am going to write my answers to the previous qual questions here. That way my professors can see and critique my answers.

-----

Fall 2009 Part 1 - Problem 1:

CLAIM: The given problem can be reduced in polynomial time to the problem that is the same except that the paths are not required to have disjoint edges (only the nonterminal vertices are required to be disjoint.)

PROOF: Suppose we modify a given instance of the problem by stating that for each edge (A,B), we add a new non-terminal vertex C and replace the edge (A,B) with edges (A,C) and (C,B). Then for any path, each edge the path used in the original graph will correspond to a vertex that the path will use in the new graph. So then the condition of vertex disjointness will include the condition of edge disjointness.

(???)

----

Fall 2009 Part 1 - Problem 2

CLAIM 0: Given any collection of numbers A_i which can collectively be represented by X bits, their sum can be represented by at most X bits.

PROOF 0: Proof by induction on n, the number of numbers in the collection.
Base Case (n=1): Trivial, since the sum is just the number itself.
Inductive Case (from n=k, prove true for n=k+1): Let Y be the number of bits needed to represent A_1 through A_k. Then their sum also can be represented with Y bits.
If A_k+1 requires less than or equal to Y bits to represent, then it is at most 2^Y - 1, and since the sum of the others is also at most 2^Y - 1, the total sum is at most 2^(Y+1)-1 and can be represented in at most Y+1 bits. Since A_k+1 requires at least one bit to represent, X>=Y+1 and the statement holds
Suppose A_k+1 requires more than Y bits to represent. We know it requires X-Y bits. Then we have the sum of the others is also at most 2^(X-Y) - 1 (since that is greater than the bound of 2^Y-1) so the total sum is at most 2^(X-Y+1)-1 and can be represented in (X-Y+1) bits. Since Y must be at least 1, (X-Y+1) <=X and the statement holds.

CLAIM 1: If A can be represented by m bits and B can be represented by n bits, then AB can be represented by at most (m+n) bits.

PROOF 1: We must have abs(A) <= (2^m)-1, and abs(B) <= (2^n)-1, so we have abs(AB) = abs(A)*abs(B) <= (2^m)-1 * (2^n) -1 <= (2^(m+n))-1.

CLAIM 1.5: Given any collection of numbers A_1 through A_n which can collectively be represented by X bits, their product can be represented by at most X bits.

PROOF 1.5 (by induction on n)
Base Case(n=1): This is trivial, since the product of a single number is just the number itself.
Inductive Case (from n=k to n=k+1): Suppose that A_k+1 can be represented in Y bits. Then A_1 through A_k can collectively be represented in X-Y. So the product of A_1 through A_k can be represented in X-Y bits by the inductive hypothesis, and by Claim 1 the product of A_1 through A_k+1 can be represented in (X-Y)+Y = X bits.

CLAIM 2: An n*n matrix that can be represented in U bits has a determinant that can also be represented in U bits.

PROOF 2: Proof by induction on n.
Base Case (n=1): This is trivial, since the determinant of a 1x1 matrix is simply that matrix's single entry.
Inductive Case (from n=k, prove true for n=k+1): Consider a (k+1)x(k+1) matrix M. Let A be the number of bits required to express the entries in the first row; then U-A is the number of bits required to express the entries in the remainder of the matrix. Now consider the expansion of the determinant by minors using the first row of the matrix. We have:

(note: M(1,i) means the given entry of M, M{1,i} means the minor)

abs(det(M))
= abs(sum(i goes from 1 to k+1) of M(1,i) * (det(M{1,i}))
<= sum(i goes from 1 to k+1) of abs(M(1,i) * det(M{1,i}))
<= sum(i goes from 1 to k+1) of abs(M(1,i) * (2^(U-A)-1)) by the inductive hypothesis, since M{1,i} can be represented in at most U-A bits
= (2^(U-A)-1) * sum(i goes from 1 to k+1) of abs(M(1,i))
<= (2^(U-A)-1) * (2^A)-1 (by Claim 0, since all the M(1,i) can be represented in a total of A bits)
<= (2^U) - 1

So det(M) can be represented in U bits.

Since U is a polynomial in U, Part A of the problem holds.

For Part B, suppose that D is the total number of bits needed to store the denominators of the rational numbers, and thus U-D is the number of bits needed to store the numerators. Let k be the product of all the denominators. We know that k can be represented in at most D bits by Claim 1.5. Consider the matrix kM. We know kM is an integer matrix since we multiplied each value by a number that it is a multiple of its denominator. Also, the number of bits needed to store kM is at most (U-D)+((n^2)D) <= (U-D)+(UD) <= U^2. Thus by Part A the number of bits needed to store det(kM) is at most U^2. We know det(M) = det(kM)/(k^n), so the determinant of M can be stored with at most U^2 bits for the numerator plus no more than nD <= UD <= U^2 additional bits for the denominator, so the total number of bits is at most 2U^2, which is a polynomial in U.

For Part C:

Given an instance P = (A,c,b,K), consider the instance Q = (dA, dc, db, dK), where d is equal to the denominator of K, times all the denominators of the elements of A, times all the denominators of the elements of c, times all the denominators of the elements of b.

CLAIM 3: The instance Q can be represented in a number of bits polynomial in the number of bits required to represent P.

PROOF 3: Suppose P can be represented in A bits. Then since d is equal to the product of all the denominators, and since the denominators can be represented in less than or equal to A bits collectively, d can be represented in less than or equal to A bits. When you multiply each of the numbers d, you are adding at most A bits to the binary representation of each number, and there are (n^2+n+n+1) <= (4n^2) numbers. So the number of bits you end up with is at most A+(4n^2)A <= (5n^2)A <= (5A^2)A = 5A^3, which is polynomial in A.

CLAIM 4: All the numbers in Q are integers.

PROOF 4: Each number in Q is a rational number in P, multiplied by d, which is a multiple of its denominator.

(???)

--------

Fall 2009 Part 1 - Problem 3

Part A:

Proof by induction on n.

Base Case (n=1)

Tuesday, January 11, 2011

A new formulation

I have come up with a new formulation of the theorem in the last post that allows us to weaken the conditions necessary (and also make the proof a bit simpler). I also cleaned up the definitions a bit. Here is the new proof.

DEFINITIONS:

A "mark set" M is defined by a function M(f) that maps faces to edges in those faces. The domain of M, denoted by d(M), is a subset of the faces in the mesh.

A face f is "set by M" if f is in d(M).

An edge e is "marked by M" if there exists a face f such that M(f)=e. (In other words, e is in the range of M.)

M is "complete" if d(M) includes all the faces in the mesh.

Two mark sets M and N are "complementary" if the union of d(M) and d(N) is all the cells in the mesh, and the intersection of d(M) and d(N) is null.

The sum M+N of two mark sets whose domains are disjoint is defined so that:

(M+N)(f) = M(f), if f is in d(M),
(M+N)(f) = N(f), if f is in d(N),
and (M+N)(f) is undefined, if f is in neither d(M) nor d(N).

M is "consistent for C" (where C is a cell) if there exists an edge e in C such that M(f)=M(g)=e, where f and g are the two faces in C that contain e.

M is "valid" if for each cell C such that all four of the faces of that cell are set by M, M is consistent for C.

CLAIM:

Consider a mesh with two complementary mark sets defined on that mesh: M and K.

If the following condition is true:

(1) For each cell C, at most two of its edges are marked by K.

Then is is possible to find a set N whose domain is identical to that of M such that:

(2A) N+K is a valid mark set. (It is also complete by construction.)
(2B) For any face x that does not contain at least one edge marked by K, we have N(x)=M(x).

PROOF:

Consider the mark set N defined (over the same domain as M) as follows:

(3A) If the face f does not contain any edges that are marked by K, then N(f) := M(f).
(3B) If the face f does contain an edge that is marked by K, then N(f) := x, where x is one of those edges. (If there are more than one of those edges, we can choose arbitrarily.)

The following lemma is easy to see is true by construction:

(L) If the face f (where f is any face in the mesh) contains an edge that is marked by K, then (N+K)(f) equals one of those edges.

Also, by construction, N satisfies (2B) above. We must now show that N+K is valid.

Consider any cell C in the mesh. There are several possibilities:

(4A) C contains zero edges that are marked by K. Then all of C's faces are set by M, and (3A) is applied for all of C's faces, so the markings of C's faces in (N+K) are the same as the markings of C('s faces in M, which we know is consistent by validity of M.

(4B) C contains exactly one edge that is marked by K; call it e. Then by L, we have (N+K)(f)=(N+K)(g)=e, where f and g are the two faces in C that contain e. Thus we know N+K is consistent for C.

(4C) C contains exactly two edges that are marked by K. If these edges do not share a vertex, they thus are not both contained in the same face, and the argument in (4B) applies to each of the edges. If the edges do share a vertex, that means that both of the edges are contained in one face, and each edge is also contained in another face that it doesn't share with the other edge. Suppose that the two edges are e and f, and e is contained in faces x and y, and f is contained in faces x and z. By (L) we have (N+K)(y) = e, and (N+K)(z) = f. We have that (N+K)(x) can equal either e or f, and we don't know which, but whichever one it equals it will still satisfy the consistency condition.

----

The way we use this proceeds as before - put all the faces that got merged in K, and put the existing markings for the rest of the faces in M.

Sunday, January 9, 2011

A performance result

DEFINITION:

We will use "edges", "faces", and "cells" to mean 1-cells, 2-cells, and 3-cells respectively. (This terminology is appropriate here because of the unique role played by each dimension of cell in this problem.)

A set of marks M is defined by a function M(x) that maps faces to the edges that are marked in them. If a face is shared by two cells, it is considered just one face.

An edge e is "marked by M" if there exists at least one cell x such that M(x)=e.

CLAIM:

Suppose we have a mesh that is validly marked - call this set of marks M. Then suppose we take a subset L of the faces in that mesh, which we will call the "locked" faces. Define a function K(x) that maps each face in L to one of its edges. (Effectively, L and K(x) define a partial set of marks. This set does not need to coincide with M in any way. We define "marked by K" in the same manner that we defined "marked by M" above.)

Also, note that M doesn't have to be a complete set of marks; we don't have to specify values of M(f) for faces f that are in L (and M only has to satisfy the condition that each cell has an edge that is marked twice for those cells that have values of M(f) for all four faces; see below).

If the following conditions are true:

(1A) On each cell, at most one face is locked.
(1B) Any non-locked face shares an edge with at most one locked cell.

Then there exists a valid set of markings N of the mesh such that:

(2A) For any face x in L, N(x) = K(x).
(2B) For any face x that does not contain at least one edge that is marked by K, N(x) = M(x).

PROOF:

Consider the function N(x) defined as follows:

(3A) If x is in L, then N(x) = K(x).
(3B) If x is not in L but it contains an edge marked by K, then N(x) equals that edge. (Note that by condition 1B, it is impossible for x to contain two edges each marked by K.)
(3C) Otherwise, N(x) = M(x).

By construction, this function satisfies conditions 2A and 2B above. We now must show that this set of markings is valid.

Consider any cell C. There are several possibilities:

(4A) Suppose none of the edges in C are marked by K. Then (3C) is used for all faces in C, so the set of markings for that cell in N is the same as that in M, and we know M is valid for C by hypothesis.

(4B) Suppose there exists an edge e in C that is marked by K, but no cells in C are in L. Then there will be two faces in C, call them f and g, that contain e. By (3B), we have N(f)=N(g)=e, so C is validly marked since two of its faces share a marked edge.

(4C) Suppose there exists a face f in C that is in L. Let e=K(f). By (2A), we know that e=N(f). Also, C has another face g that also contains e. By (1A), g is not locked, so by (3B) we have e=N(g), which means f and g share a marked edge and thus the cell is validly marked.

-----

So, what does all this mean for our original problem? Well, let's say we have our original mesh, and then we contract an edge. Put all the faces that got "merged" by that edge contraction in L, and for the rest, keep the original markings and call it M. Now whatever we set the markings in L to, we can guarantee that we won't have to change the markings on any cell that isn't adjacent to one of the edges we just marked in L. This means that we have an effective upper bound on how many cells we will have to change, and also means that we have a simpler optimization problem since we just have to look at how to mark the cells in L, and let the algorithm above do the rest.

Of course this is contingent on the conditions 1A and 1B above being true, but those shouldn't be too hard to satisfy.

Saturday, January 1, 2011

More on the graph method

So I thought a little more about how the graph method (number 3 in the list in the last post) works when you are doing edge contractions. In an edge contraction, each pair of faces on opposite sides of the contracted edge gets merged into one face. That means that if the two faces had different edges marked, then you have to switch one of them so it matches the other edge. This can be simulated in the graph method by having an extra node for each of these pairs that has out-degree 2, one going to each of one of the possible switches. The trickier thing is dealing with the fact that there are multiple of these forced switches at the same time. Suppose that on one of the cells we have to do a given marking switch due to the edge contraction, then we propagate it back until eventually we hit another marking switch that we had to make due to the same edge contraction. This means that we have made both of those changes at one time.

When you put it all together this gives an algorithm for the problem:

1. When you create the mesh and initial markings, also create the graph. You also have to update the graph each time you refine, coarsen, or change the markings - but this only affects the cells where the markings were changed so it takes (on average) constant time.

2. When you contract an edge, make a list of all the pairs of opposing faces that don't have edge markings that correspond.

3. While that list still has elements in it:
3a. Do the search process starting from one of the elements in the list (chosen arbitrarily), searching for either (a) one of the other elements in the list or (b) a node with out-degree 0.
3b. Make the edge switches called for by the result of the search, updating the graph as necessary.
3c. Remove the chosen element form the list. Also if the search process ended at one of the other elements in the list, remove that element as well.

This by itself is probably not the best algorithm we can come up with because it is likely to change the same cells over and over, and also looks that each required change one by one rather than considering them all at the same time. But it is a simple algorithm that provides a baseline that we can use as a point of comparison to evaluate more sophisticated options. Probably the next thing I will do is implement this algorithm (not necessarily in the TentPitcher codebase, just a prototype implementation in Python or something) to see how well it works.

Thursday, December 30, 2010

Mesh Adaptivity

The research problem we are trying to solve right now is based on this paper. The paper describes a tetrahedral mesh adaptivity scheme based on "edge markings". Each face (2-cell) of a tetrahedron has one of its edges (1-cells) marked. A set of markings is valid if the following are true:

1. Whenever two tetrahedrons share a face, that face is marked the same way in each of the tetrahedrons. (That is, each face has a unique marking, even if it's in more than one tetrahedron.)

2. In each tetrahedron, there exists an edge that is marked on both of the faces that contain it.

The paper shows how to use these markings to do refinement. The problem we are trying to figure out is how to do coarsening. In particular, when you coarsen the mesh by contracting an edge to nothing, you merge pairs of opposite faces into one face, so the question is how to change the markings so that it is still valid. We want to do this in a way that involves changing the fewest number of markings.

There are a few different ways we can go about doing this:

1. The method described in the paper about how to set up the initial markings says that we can establish a total order over all the edges, then have each face mark the edge that is first in that order. This raises the question of whether or not this property - that there exists a total order over the edges such that the marked edge in each face is the first in that total order - is preserved over the refinement process. If so, then we have a simple algorithm for merging two faces when doing refinement: just choose one of the faces arbitrarily and have the new edges inherit that edge's position in the order, then remark them according to that order. This obviously will not affect any tetrahedron that doesn't contain one of the edges that was just merged, so we can keep the re-marking down to a small number of tetrahedrons.

When we generate a tetrahedron of type P, all the marked edges don't include the new vertex, so the property is preserved if we just put all the new edges last in the partial order. When we generate a tetrahedron of type A, it is a little more complicated - the property is still preserved but we have to put some of the edges between other edges in the order (I'll put a diagram up soon.) There is still a complication of when we generate a new tetrahedron where the new vertex was already generated (by refining a different tetrahedron that the refined edge was also a part of) because now we have to make sure that the same position for that edge in the order will satisfy the conditions on both of the faces I'm not sure yet if that's actually satisfied.

2. We could look at relations between this problem and the Boolean satisfiability problem (SAT). This problem has a very natural Boolean-satisfiability representation: the proposition "face F has edge E marked" is one of the variables for each pair (F,E), then the conditions that each face has exactly one edge marked as well as the conditions (1) and (2) above are easily expressed. If we wanted to try to minimize the number of changes we have to make then we could think of it as a sort of weighted maximum satisfiability problem where we try to also satisfy as many of the propositions "marking X wasn't changed" as possible. We could then use a MaxSAT solver like this one to figure out the best solution. This direction is probably not the best way of doing it because we have to convert the problem to a different form, solve it, and convert the answer back, plus the size of the problem will be the size of the whole mesh and we really want to only look at a small neighborhood of cells (for parallelization purposes if nothing else).

Another approach is to try to reduce the satisfiability problem to our problem; that would show that it is NP-complete. Of course that wouldn't give us an actual algorithm; it would just tell us we should stop looking for an exact polynomial-time algorithm. So this would be nice to have but isn't the first direction we should look.

3. Another way to think about the problem is to imagine each possible "switch" of an marking on one face from one edge to another edge as a node in a directed graph. Suppose that we switch a marking on a face of one tetrahedron to make the markings on that tetrahedron valid. That may make the markings on the adjacent tetrahedron (the one on the "other side" of that face) not valid because we also have to switch that marking on that tetrahedron. That may require switching yet another marking, and so on. Note that any position is only one "marking switch" away from a valid position, so we only have to do one marking switch at each step. In the directed graph, for each marking switch that necessitates another marking switch, put an edge to the possible choices of marking switches to make it valid, and for those marking switches that will make it valid on both sides and allow you to halt the process, mark those as "end nodes." Then the problem is just to find the shortest path from the change you are forced to make due to the merging to an end node, which is a simple matter of breadth-first search (or Dijkstra's algorithm, if we want to weight the marking switches by "how hard" they are to make). It will probably end up a little more complicated because an edge contraction affects lots of faces rather than just one.

I'll try to post up some more ideas in the next couple days.

Monday, December 20, 2010

Masters or Ph.D?

I just talked to Jeff Erickson last week, and he gave me some very interesting feedback. He said that I am very smart and willing to work very hard, and that I am able to get a lot done, provided that I am give a specific problem to work on. However, the one thing that I haven't shown myself to be abole to do yet is to make independent decisions about which directions to take my research and what are important research problems to work on, and that's what I will need to be able to do to get a Ph.D. He suggested that a better direction is to focus on getting a Masters' thesis first, and then at that point re-evaluating whether I want to continue toward a Ph.D. or just stop there.

The more I think about it, the more I think that a Masters may be the best option for me, for a few reasons:

1. First, one thing about getting a Ph.D. in computer science is that it doesn't actually involve writing that much code; or at least code isn't the primary product. As my advisor said, the main "product" you're producing as a Ph.D. student (and later as a professor) is ideas; code is just the way of implementing those ideas. But I think my specialty isn't necessarily really in developing theory, but more in the combination of math/theory and code. In particular, there are lots of people who are experts in the theory (or who are experts in a specific problem domain, like engineers) but aren't that good at actually writing code, and there are lots of people who can crank out lots of code without really understanding the theory, but I am very good at both understanding the theory and algorithms and translating it into code. So my main niche is to bridge that gap.

2. In particular, when I look at the things that really sound cool to me and "wow, I would really like to be working on this," it's generally not the theoretical stuff; it's more the stuff that involves taking a real-world problem, one not necessarily in computer science, and figuring out how to formalize and understand it in a mathematical context so you can write code to solve it. For example:

2a. At my work at the U.S. Census Bureau, we were studying disclosure avoidance in survey data - that is, how to disseminate data from local level surveys in ways that preserve the statistical patterns in the data but don't allow users to figure out who answered what to any particualr question. This question involves first defining specifically what it means for a respondent's privacy to be preserved (What if an adversary can make a good guess at the missing information but only some of the time? What if an adversary has access to a different database unrelated to the survey that they can combine with that information? Is there a way to prove things about what an adversary can figure out?), then to figure out what algorithms can implement that, then to code and test those algorithms.

2b. I was offered a job at Amazon.com after graduating from University of Maryland, though I chose not to take it so I could go to graduate school. One of the projects they were working on there was figuring out how to optimize their web pages to increase sales - so they had collected all this data about how the page layout and the time it takes to load the page affect users' behavior. Again, we have a practical problem, then there's a step of formalizing the problem (Exactly what function are you trying to optimize? Is it just number of sales, or are other things like which products get sold or how long it takes the user to go through the process also important? How can you determine which changes in output are caused by which factors? What data can be used to shed light on these questions?), figuring out what algorithms can be used to analyze the data, and then implementing and testing these algorithms.

2c. Another example of something I tried to do once was when I was learning how to produce foam tipped arrows for Belegarth (a game in which players dress up as medieval warriors and fight each other with padded replicas of medieval weapons). There are several factors that you need to consider when making the arrow, such as the type of foam to use, the weight of the arrow, how hard or soft the tip is, and so on. My first response was to try to figure out if there was a way to mathematically model what happens when the arrow strikes the target so as to better understand why the design is the way it is and if there is any way to improve it. I did some research and soon concluded that the formulas themselves wouldn't be overly complicated but that the data you would need to calibrate the model isn't readily available. But you can see this is another example of how I like to try to understand a problem by putting it into a mathematical framework.

3. Of course, most of the problems above aren't actually computer science research problems; they're engineering or business problems that have a math/computer science component. But I think that might be the kinds of things I want to work on: to work closely with the engineers or business people to understand what the computational needs are, then go and figure out what computer science techniques can meet those needs*. I talked to Bob Haber about this and he said that if I wanted to actually "do an application" for my Ph.D. thesis, I would essentially have to go back and take some undergraduate level engineering (or physics, or whatever) courses to learn all the field-specific stuff I would need to know to do research in that area. So if I want to do this kind of stuff, it would make much more sense to work at a lab or company than in academia.

So, the plan we have come up with is as follows. I'll work toward a Masters, and at the same time I'll also start looking for internships and jobs elsewhere. Los Alamos National Laboratories and Sandia National Laboratories do lots of scientific computing stuff so they are good possibilities. Then once I get a Masters, I will decide whether I want to stop there or whether I want to keep going.

*Although, I might be naively optimistic as to how much sophisticated algorithmic stuff is actually involved in these types of problems. I am reminded of a professor from University of Maryland (I think it was Bill Gasarch but I'm not sure) who, when he was teaching about job-scheduling algorithms, planned to bring his wife (who worked in a machine shop) in to talk about all the exciting new techniques they use to optimize their scheduling to improve production. He quickly cancelled this plan when she informed him that the actual algorithm they use is "whoever yells the loudest gets their job scheduled first."