Sunday, November 29, 2015

Shortest path touching a circle


This came up a week or so ago in a discussion with a lab mate. The governing principle is well-known from physics, but I was too lazy to find a derivation online. Later on, I had an idea and wanted to give it a try, so here goes.

We're interested in the shortest path starting at point $a$ to point $b$ with the requirement that it intersects a circle of radius $r$ centered at point $c$, and does so before stopping at $b$. When the circle intersects the line segment $\overline{ab}$, the segment coincides with the shortest path. Hence, we are interested in paths that reflect off the circle, just like a ray of light does when it hits a shiny surface. I'll be referring to the figure to the right.

Eventually, we're interested in a statement involving the angle of incidence $\theta$. However, I found it easier to work instead with the central angle $\theta'$. We might as well parameterize the path by the point of incidence $x$ and write the length as:

\[ f(x) = |\overline{ax}| + |\overline{xb}|. \]

Now, what do we know about each of $\overline{ax}$ and $\overline{xb}$? Assume $d_a$ and $d_b$ are the distances from $c$ to $a$ and $b$, respectively. Let's look at $\triangle{cxa}$. By the law of cosines, we have that:

\[ |\overline{ax}|^2 = d_a^2 + r^2 - 2  d_a  r  \cos{\theta'}. \]
Likewise for $\triangle{cxb}$, where $\phi' = \angle{acb}$, we get:
\[ |\overline{xb}|^2 = d_b^2 + r^2 - 2  d_b  r  \cos{(\phi' - \theta')}. \]

It's quite disappointing that we don't obtain nice expressions for $|\overline{ax}|$ and $|\overline{xb}|$ to plug into $f(x)$. But, since we're only interested in minimizing $f$, we might hope that finding the derivative will get us around that. Let's start by rewriting $f$ as:

\[ f(x) = \sqrt{|\overline{ax}|^2} + \sqrt{|\overline{xb}|^2}. \]

We can then write the derivative as:

\[ \frac{\partial f(x)}{\partial \theta'} = \frac{1}{2\sqrt{|\overline{ax}|^2} }\frac{\partial |\overline{ax}|^2}{\partial \theta'} + \frac{1}{2\sqrt{|\overline{xb}|^2}.}\frac{\partial |\overline{xb}|^2}{\partial \theta'}. \]

Plugging in $\frac{\partial |\overline{ax}|^2}{\partial \theta'} = 2d_ar\sin{\theta'}$ and $\frac{\partial |\overline{xb}|^2}{\partial \theta'} = -2d_br\sin{(\phi'-\theta')}$ we get:

\[ \frac{\partial f(x)}{\partial \theta'} = \frac{d_ar\sin{\theta'}}{|\overline{ax}|} - \frac{d_br\sin{(\phi' - \theta')}}{|\overline{xb}|}. \]

Setting $\frac{\partial f(x)}{\partial \theta'} = 0$, cancelling $r$ and rearranging, we find that:

\[ \frac{d_a\sin{\theta'}}{|\overline{ax}|} = \frac{d_b\sin{(\phi' - \theta')}}{|\overline{xb}|}. \]

It would really help to realize that $d_a\sin{\theta'} = |\overline{aa'}|$ and $d_b\sin{(\phi'-\theta')} = |\overline{bb'}|$, where $a'$ and $b'$ are the projections of $a$ and $b$, respectively, on $\overrightarrow{cx}$. Making this substitution we get:

\[ \frac{|\overline{aa'}|}{|\overline{ax}|} = \frac{|\overline{bb'}|}{|\overline{xb}|}. \]

Letting $\phi = \angle axb$, this is the same as:

\[ \sin{\theta} = \sin{(\phi - \theta)}. \]

Or,

\[ \theta = \phi - \theta. \]

Which says that the angle of incidence should be equal to the angle of reflection, for the path to be optimal. Which is already known.

Update 01/14/16: There remains the issue of determining the point $x$ that minimizes the shortest path defined by $f$. I found it easier to work with angles $\alpha = \angle{cax}$ and $\beta = \angle{cbx}$. Applying the law of sines in $\triangle{cxa}$ and $\triangle{cxb}$ we find that:

\[ \sin{(\pi - \theta)} = \sin{\theta} = \frac{d_a}{r}\sin{\alpha}, \quad \sin{(\phi - \theta)} = \frac{d_b}{r}\sin{\beta}.\]

Since the point $x$ we are interested in makes $\theta = \phi - \theta$, we get:

\[ d_a\sin{\alpha} = d_b\sin{\beta}. \]

In a way, this equation only describes the angle of incidence from $a$ or $b$ to points on the circle. We still need a second equation to ensure these two angles correspond to a single point of incidence. This can be achieved by examining $\triangle{abx}$. Letting $\alpha' = \angle{cab}$ and $\beta' = \angle{cba}$ we can write:

\[ (\alpha' - \alpha) + (\beta' - \beta) + \phi = \pi. \]

Rewriting $\phi = 2 \theta$ and rearranging we get:

\[ \beta = (\alpha' + \beta' - \pi) - \alpha + 2\theta. \]

Noting that $\theta = \arcsin{(\frac{d_a}{r}\sin{\alpha})}$, we can now eliminate $\beta$ using the first equation to get an equation for $\alpha$ only, and the same can be written for $\beta$:

\[ \alpha = (\alpha' + \beta' - \pi) + 2 \arcsin{(\frac{d_a}{r}\sin{\alpha})} - \arcsin{(\frac{d_a}{d_b}\sin{\alpha)}}. \]

Other than a numerical approach, it is plausible to assume $\alpha$ is small especially when $d_a \gg r$. This might justify linearizing both the sines and arcsines to get something like:

\[ \alpha \approx \frac{\alpha' + \beta' - \pi}{1 + \frac{d_a}{d_b} - 2 \frac{d_a}{r}}. \]

Sunday, March 23, 2014

PS to PDF: convert and crop

I've done this task once before and forgot how. It took me a while to do it again, so I'm posting this here in case somebody else needs it. It might not be the best way to do it so if you can make it better, I'd love to hear from you.

I'm using pstill_dist and pdfcrop with two subdirectories for all input and output files.
#!/bin/bash
for FILE in $(find ./input -name '*.ps'); do
    FNAME="${FILE##*/}"
    FNAME="${FNAME%.*}"
    ./pstill "./input/${FNAME}.ps" -o "./output/${FNAME}.pdf"
    pdfcrop --clip "./output/${FNAME}.pdf" "./output/${FNAME}.pdf"
done

Sunday, May 8, 2011

Sorting Facebook feeds in chronological order

College students spend lots of time on Facebook and virtually use it to share and organize their shcedules and course materials. Obviously, things get out of control and you end up with a very long and entangled feed that you might find yourself in a real need to explore or even worse, scan for something important.

So, one day a student of mine posted (!), wondering if there were any group settings that would arrange posts in reverse chronological order, presumably to save some time. That post caught my attention and eventually, I decided to prepare a bookmarklet for this exact purpose. It was pretty much straightforward but unfortunately, it would only work for groups. Someone may find it worthwhile to adapt it to the home/profile feeds, though.

Update: It seems the student got interested in this sort of JS magic and she actually took it a step further to work for the home feed as well.
// This snippet is licensed under a Creative Commons Attribution 3.0 License.
// Authors: Ahmed Abdelkader and Nancy Iskander
javascript:{
    var flist = document.getElementById('home_stream');
    if (flist == null)
    flist = document.getElementById('pagelet_group_mall').children[0].children[0];
    var far = new Array();
    var dates = {};
    while (flist.children.length) {
        var abbrs = flist.children[0].getElementsByTagName('abbr');
        if (abbrs.length) {
            dates[flist.children[0].id] = abbrs[0].getAttribute('data-utime');
            far.push(flist.children[0]);
        }
        flist.removeChild(flist.children[0]);
    }
    far.sort(function(a, b){return dates[b.id] - dates[a.id];});
    for (var i = 0; far.length - i; ++i)
        flist.appendChild(far[i]);
    void(0);
}
To create the bookmarklet, we just need to compress that and put it in a hyperlink. Just drag this link: Feed-Repairo!, into your browser's Bookmarks Toolbar to add the bookmark and you can easily rename it if you please. That's it! Whenever you need to sort out your group feeds, all you need to do is click the bookmark. Let me know if it works for you or you encoutner any bugs. Enjoy!

Update 8/23/2013: Drag this link if you want to sort ascendingly Repairo-Feed!
Update 3/23/2014: flist is not at the root container anymore; so we added another .children[0] and all is fine. (Many thanks to Jeff E)

Wednesday, January 26, 2011

A lightweight anonymous visitor for Java collections

I got myself into a situation where I really wished that Java collections accepted a visitor interface of some kind so I could define one on the fly and get my job done. Alas, I decided to implement the utility myself and here's the result.

public interface IVisitor<T> {
void visit(T item);
}

public class CollectionUtils {
public static<T> void applyVisitor(Collection<T> collection, IVisitor<T> visitor) {
for (T item : collection)
visitor.visit(item);
}
}
A typical usage example would be:
public static void main(String[] args) {
CollectionUtils.applyVisitor(Arrays.asList(1, 2, 3, 4, 5), new IVisitor<Integer>() {
public void visit(Integer x) {
System.out.println(x);
}
});
}

Sunday, October 24, 2010

Efficient enumeration of all integers with a given pop count

Population count is the number of '1' bits in the binary representation of an integer.

The basic trick we employ is to get the next higher number with the same number of 1-bits, which we borrow from the Hacker's Delight (whole book). Below is a Java implementation of all 4 version of the snoob function.

public static int snoob1(int x) {
int r = x + (x & -x);
return r | ((x ^ r) >>> (2 + Integer.numberOfTrailingZeros(x)));
}

public static int snoob2(int x) {
int s = x & -x, r = x + s;
return r | ((x ^ r) >>> (33 - Integer.numberOfLeadingZeros(s)));
}

public static int snoob3(int x) {
int r = x + (x & -x);
return r | ((1 << (Integer.bitCount(x ^ r) - 2)) - 1);
}

public static int snoob4(int x) {
int s = x & -x, r = x + s;
return r | (((x ^ r) >> 2) / s);
}
We utlize this formula to implement the method below, where the width parameter is added for convenience.
public static void printAllIntegersOfWidthAndPop(int w, int n) {
if (n <= 0 || n > w || w > 32) return;
int x = (1 << n) - 1, g = x << (w - n), c = 0;
while (true) {
String s = Integer.toBinaryString(x);
while(s.length() < w) s = "0" + s;
System.out.println(++c + "\t" + s);
if (x == g) break;
else x = snoob*(x);
}
}
It is easy to figure out there will be Choose(w, n) such numbers. We know where the sequence starts and since we also know where it stops, we don't need to compute this number.

If you'd like to learn how to find the number of leading/trailing zeros or the pop count or learn about more cool bit twiddling hacks you can consult the book above (which the java.lang.Integer implementation of these methods is based on) or this magnificent web page.

Saturday, July 17, 2010

Generating combinations in lexicographical order using Java

Below is a Java port of the algorithm by James McCaffrey as presented in the MSDN article Generating the mth Lexicographical Element of a Mathematical Combination.

Please do take into consideration that this implementation only uses int as Java does not allow array indexing with long. This means that the valid input range is more restricted. Watch out for overflows!

/**
* Based on the Combinadic Algorithm explained by James McCaffrey
* in the MSDN article titled: "Generating the mth Lexicographical
* Element of a Mathematical Combination"
* <http://msdn.microsoft.com/en-us/library/aa289166(VS.71).aspx>
*
* @author Ahmed Abdelkader
* Licensed under Creative Commons Attribution 3.0
* <http://creativecommons.org/licenses/by/3.0/us/>
*/
public class Combinations {
/** returns the mth lexicographic element of combination C(n,k) **/
public static int[] element(int n, int k, int m) {
int[] ans = new int[k];

int a = n;
int b = k;
int x = (choose(n, k) - 1) - m; // x is the "dual" of m

for (int i = 0; i < k; ++i) {
a = largestV(a, b, x); // largest value v, where v < a and vCb < x
x = x - choose(a, b);
b = b - 1;
ans[i] = (n - 1) - a;
}

return ans;
}

/** returns the largest value v where v < a and Choose(v,b) <= x **/
public static int largestV(int a, int b, int x) {
int v = a - 1;

while (choose(v, b) > x)
--v;

return v;
}

/** returns nCk - watch out for overflows **/
public static int choose(int n, int k) {
if (n < 0 || k < 0)
return -1;
if (n < k)
return 0;
if (n == k || k == 0)
return 1;

int delta, iMax;

if (k < n - k) {
delta = n - k;
iMax = k;
} else {
delta = k;
iMax = n - k;
}

int ans = delta + 1;

for (int i = 2; i <= iMax; ++i) {
ans = (ans * (delta + i)) / i;
}

return ans;
}
}
The code below produced the output that follows:
public static void main(String[] args) {
int n = 5, k = 3;
int total = choose(n, k);
for(int i = 0; i < total; i ++) {
for(int x : element(n, k, i))
System.out.print(x + " ");
System.out.println();
}
}

0 1 2
0 1 3
0 1 4
0 2 3
0 2 4
0 3 4
1 2 3
1 2 4
1 3 4
2 3 4

Monday, June 21, 2010

G/G/1 Queue Model Simulation in Java

In this post, we present an object-oriented design and implementation of an event-driven G/G/1 queue model simulation using Java. The lecture notes on Computer System Analysis by Raj Jain were very helpful and are highly recommended. In particular, we would like to reference the introductory lectures on Simulation Modeling and Queueing Theory.

Please note that this implementation is provided as a guide/starting point and is not, by any means, complete. We preferred a simple design while keeping in mind where you may wish to extend it.

Update: Download the Eclipse project here.

An event-driven simulation operates by generating and processing events through time. For our queue model, we have two types of events: arrivals and departures or service completion. It is also necessary to determine which events should happen first.

public class Event implements Comparable<Event> {
protected double time;
protected int code;

...

public int compareTo(Event e) {
return Double.compare(time, e.time);
}
}
While the design will not make any assumptions about the interarrival or service time distributions, we still need a suitable way to represent them and a couple of concrete implementations for testing.
public abstract class Distribution {
public abstract double generateRV();
}

public class UniformDistribution extends Distribution {
double a, b;

...

public double generateRV() {
return a + rand.nextDouble() * (b - a);
}
}

public class ExponentialDistribution extends Distribution {
double lambda;

...

// Generating exponential variates.
public double generateRV() {
return -1/lambda * Math.log(rand.nextDouble());
}
}

public class NormalDistribution extends Distribution {
double mu, segma;

...

public double generateRV() {
return mu + rand.nextGaussian() * segma;
}
}
Now, we should be ready to implement the queue class. We preferred to keep the queue simple by moving data collection into a separate arbitrary observer.
/**
* @author Ahmed Abdelkader
* Licensed under Creative Commons Attribution 3.0
* <http://creativecommons.org/licenses/by/3.0/us/>
*/
public class GG1Q {
protected Distribution arrivalDistribution;
protected Distribution serviceDistribution;
protected double t;
protected PriorityQueue<Event> eventQueue = new PriorityQueue<Event>();
protected int customers;
protected QueueObserver observer;

...

public void run(double duration) {
eventQueue.add(new Event(t + arrivalDistribution.generateRV(), Event.ARRIVAL));
while(t < duration) {
Event e = eventQueue.poll();
t = e.time;
switch(e.code) {
case Event.ARRIVAL:
customers++;
eventQueue.add(new Event(t + arrivalDistribution.generateRV(), Event.ARRIVAL));
if(customers == 1)
eventQueue.add(new Event(t + serviceDistribution.generateRV(), Event.SERVICE));
break;
case Event.SERVICE:
customers--;
if(customers > 0)
eventQueue.add(new Event(t + serviceDistribution.generateRV(), Event.SERVICE));
break;
}
if(observer != null)
observer.stateChanged(this, e);
}
}

...
}
We define the observer interface and implement two sample observers: one for estimating the queue length PDF and the other for the waiting time CDF. A composite observer allows more than one observer to collect data of a single run.
public interface QueueObserver {
void stateChanged(GG1Q q, Event e);
void printStats();
}

public class CompositeQueueObserver implements QueueObserver {
protected ArrayList<QueueObserver> observers = new ArrayList<QueueObserver>();

public void addObserver(QueueObserver observer) {
observers.add(observer);
}

public void stateChanged(GG1Q q, Event e) {
for(QueueObserver observer : observers)
observer.stateChanged(q, e);
}

public void printStats() {
for(QueueObserver observer : observers)
observer.printStats();
}
}

public class QueueLengthQueueObserver implements QueueObserver {
protected TreeMap<Integer, Integer> lengthFrequencies = new TreeMap<Integer, Integer>();
protected int total;

public void stateChanged(GG1Q q, Event e) {
if(e.getCode() != Event.ARRIVAL) return;
Integer length = lengthFrequencies.get(q.getLength());
if(length == null) length = 0;
lengthFrequencies.put(q.getLength(), length + 1);
total++;
}

public void printStats() {
System.out.println(getClass().getSimpleName());
for(int length : lengthFrequencies.keySet())
System.out.println(length + "\t" + 1.0*lengthFrequencies.get(length)/total);
}
}

public class WaitingTimeQueueObserver implements QueueObserver {
protected ArrayList<Double> waitingTimes = new ArrayList<Double>();
protected int arrivalIndex, serviceIndex;
protected HashMap<Integer, Double> arrivalTimes = new HashMap<Integer, Double>();

public void stateChanged(GG1Q q, Event e) {
switch(e.getCode()) {
case Event.ARRIVAL:
arrivalTimes.put(arrivalIndex++, e.getTime());
break;
case Event.SERVICE:
waitingTimes.add(e.getTime() - arrivalTimes.get(serviceIndex++));
break;
}
}

public void printStats() {
System.out.println(getClass().getSimpleName());
if(waitingTimes.size() == 0) return;
double acc = 0;
Collections.sort(waitingTimes);
for(int i = 1, j = 0; i <= 10; i++) {
while(acc < i/10.0 && j < waitingTimes.size()) {
acc += 1.0/waitingTimes.size();
j++;
}
System.out.println(waitingTimes.get(j-1) + "\t" + i/10.0);
}
}
}
The following test program produced the output below:
public class Main {
public static void main(String[] args) {
GG1Q q = new GG1Q(new ExponentialDistribution(100.0/3600), new UniformDistribution(1, 10));
QueueObserver observer = new QueueLengthQueueObserver();
q.setObserver(observer);
q.run(1000);
observer.printStats();
}
}

QueueLengthQueueObserver
Queue length PDF:
1 0.47058823529411764
2 0.23529411764705882
3 0.058823529411764705
4 0.058823529411764705
5 0.11764705882352941
6 0.058823529411764705