Standard Forth versus Python: a case study

[John Doty]

I realized that I have a little job on the table that is a fine test of
the Python versus Standard Forth code availability and reusability issue.

Note that I have little experience with either Python or Standard Forth
(but I have much experience with a very nonstandard Forth). I've noodled
around a bit with both gforth and Python, but I've never done a serious
application in either. In my heart, I'm more of a Forth fan: Python is a
bit too much of a black box for my taste. But in the end, I have work to
get done.

The problem:

I have a bunch of image files in FITS format. For each raster row in
each file, I need to determine the median pixel value and subtract it
from all of the pixels in that row, and then write out the results as
new FITS files.

This is a real problem I need to solve, not a made-up toy problem. I was
originally thinking of solving it in C (I know where to get the pieces
in that language), but it seemed like a good test problem for the Python
versus Forth issue.

I looked to import FITS reading/writing, array manipulation, and median
determination. From there, the solution should be pretty easy.

So, first look for a median function in Python. A little googling finds:

http://www.astro.cornell.edu/staff/loredo/statpy/

Wow! This is good stuff! An embarrassment of riches here! There are even
several FITS modules, and I wasn't even looking for those yet. And just
for further gratification, the page's author is an old student of mine
(but I'll try not to let this influence my attitude). So, I followed the
link to:

http://www.nmr.mgh.harvard.edu/Neural_Systems_Group/gary/python.html

From there, I downloaded stats.py, and the two other modules the page
says it requires, and installed them in my site-packages directory. Then
"from stats import median" got me a function to approximately determine
the median of a list. It just worked. The approximation is good enough
for my purposes here.

Pyfits required a little more resourcefulness, in part because STSCI's
ftp server was down yesterday, but I got it installed too. It helps that
when something is missing, the error message gives you a module name. It
needs the numarray module, so I got array manipulation as a side effect.

I haven't finished the program, but I've tried out the pieces and all
looks well here.

OK, now for Forth. Googling for "forth dup swap median" easily found:

http://www.taygeta.com/fsl/library/find.seq

At first blush, this looked really good for Forth. The search zeroed in
on just what I wanted, no extras. The algorithm is better than the one
in the Python stats module: it gives exact results, so there's no need
to check that an approximation is good enough. But then, the
disappointment came.

What do you do with this file? It documents the words it depends on, but
not where to get them. I'm looking at a significant effort to assemble
the pieces here, an effort I didn't suffer through with Python. So, my
first question was: "Is it worth it?".

The answer came from searching for FITS support in Forth. If it exists
in public, it must be really well hidden. That's a "show stopper", so
there was no point in pursuing the Forth approach further.

In the end, it was like comparing a muzzle-loading sharpshooter's rifle
with a machine gun: Forth got off one really good shot, but Python just
mowed the problems down.

The advocates of the idea that Standard Forth has been a boon to code
reusability seem mostly to be people with large private libraries of
Forth legacy code. No doubt to them it really has been a boon. But I
think this little experiment shows that for the rest of us, Python has a
published base of reusable code that puts Forth to shame.

 

[Paul Rubin]

I have a bunch of image files in FITS format. For each raster row in
each file, I need to determine the median pixel value and subtract it
from all of the pixels in that row, and then write out the results as
new FITS files.

I dunno what FITS is, but if you have a list of pixel values, that
calculation sounds like two lines:

median = sorted(pixels)[len(pixels)//2]
new_pixels = [p-median for p in pixels]

 

[John Doty]

I have a bunch of image files in FITS format. For each raster row in
each file, I need to determine the median pixel value and subtract it
from all of the pixels in that row, and then write out the results as
new FITS files.

I dunno what FITS is, but if you have a list of pixel values, that
calculation sounds like two lines:

median = sorted(pixels)[len(pixels)//2]
new_pixels = [p-median for p in pixels]

Yes. The efficient exact algorithms for this problem use *partial*
sorts. The Forth one from the FSL is of this class (although I know of
two better ones for big arrays). But it's tough to beat the efficiency
of the approximate histogram-based method the Python stats module
implements.

 

[bearophileHUGS]

John Doty:

Yes. The efficient exact algorithms for this problem use *partial*
sorts. The Forth one from the FSL is of this class (although I know of
two better ones for big arrays). But it's tough to beat the efficiency
of the approximate histogram-based method the Python stats module
implements.

The usual way to compute a true median with Python may be:

def median(inlist):
newlist = sorted(inlist)
index = len(newlist) // 2
if len(newlist) % 2:
return newlist[index]
else:
return (newlist[index] + newlist[index-1]) / 2.0

If you can use Psyco and your FITS lines are really long (well, maybe
too much, the treshold if about >~3000 in my PC) you can use something
like this instead the builtin timsort:
http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/466330
(To compute the median on a image line the median that uses sort is
probably better in most cases, expecially if you use the built in sort
of numerical libraries.)

Bye,
bearophile

 

[jacko]

John Doty:

Yes. The efficient exact algorithms for this problem use *partial*
sorts. The Forth one from the FSL is of this class (although I know of
two better ones for big arrays). But it's tough to beat the efficiency
of the approximate histogram-based method the Python stats module
implements.

The usual way to compute a true median with Python may be:

def median(inlist):
newlist = sorted(inlist)
index = len(newlist) // 2
if len(newlist) % 2:
return newlist[index]
else:
return (newlist[index] + newlist[index-1]) / 2.0

If you can use Psyco and your FITS lines are really long (well, maybe
too much, the treshold if about >~3000 in my PC) you can use something
like this instead the builtin timsort:
http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/466330
(To compute the median on a image line the median that uses sort is
probably better in most cases, expecially if you use the built in sort
of numerical libraries.)

Bye,
bearophile

partial bucket sort with quicksort of individual bucket needed for
index list.

APL would be fast, try a solution in J

it calculates need on demand i think, and so the calculation dependance
tree only does the one quicksort on the bucket needed, or both if on a
bucket boundry, but this can be avoided with clever bucket selection.

cheers

 

[idknow]

John Doty:

Yes. The efficient exact algorithms for this problem use *partial*
sorts. The Forth one from the FSL is of this class (although I know of
two better ones for big arrays). But it's tough to beat the efficiency
of the approximate histogram-based method the Python stats module
implements.

The usual way to compute a true median with Python may be:

def median(inlist):
newlist = sorted(inlist)
index = len(newlist) // 2
if len(newlist) % 2:
return newlist[index]
else:
return (newlist[index] + newlist[index-1]) / 2.0

[snip]

no sort() is needed to calculate the median of a list.

you just need one temp var.

 

[Terry Reedy]

(e-mail address removed) wrote:
no sort() is needed to calculate the median of a list.
you just need one temp var.

Of course. But for a short enough list, the builtin sort() method may be
faster than an O(n) algorithm coded in Python. And .sort() is already
present.

tjr

 

[Paul McGuire]

[snip]

no sort() is needed to calculate the median of a list.

you just need one temp var.

Ok, I'll bite. How do you compute the median of a list using just a single
temp var?

-- Paul

 

[werty]

Apples/oranges ? programmers are making very little $$ today .
Thats software ! No one is makin money on obsolete Forth ,
so why a comparisom ?

Ultimately the best OpSys will be free and millions of lines of code
obsoleted . Because no one can protect intellectual property , its
simply too costly for a Government to do this .

Notice the crypto fiasco , the USA govt forbad export of SSL and
in short order Australia gave the world a free SSLeasy !

This is software . No one will , for long , own software .
Microsoft and Linux will die in 24 months . No hardware will
produce money/profits using them .

Because we have a new hardware competitor , that obsoletes the BIOS
chip , North bridge/South bridge USB helpers , NICs .
No mo PNP (Plug_and_Pray) .
In 2 chips ( CPU /SoC and LCD controller) it will be faster than
a PC . 100's have allready written HDD drivers and haven't yet
free'd them . But when others give free , what good do it to
hold yours ! You look stupid selling what others improve and
free . Trying to sell what others find easy to create !
Intel made hardware too hard to program , ARM is easy .
Python wont last on ARM , for WCE and Linux will be tossed .
A simpler structured method will be used to customise
Browsers . There will be 2 columns , one on left will be main/orig
but on the Right will have hyperlink result . This way ya dont have
to go back to compare ! side by side .
Text editors will also work this way . You will read orig in left
columns
and 'edit'/annotate your stuff to the Right columns .

Sun microsystems et all ( the ones we seldom think about ) will
all be obsoleted , bankrupted , their hardware surplused .
No more h/w servers .
There will be no high speed 64 bit boxes in the future .
The Pocket PC will do work you could not imagine !
All will have 100 GB HDD , VGA LCD , USBH/USBD , WIFI N
and no WERTY keyboard ! full GUI ...
ethernet and firewire and rs232 will die

Why not "see" the future ?
No Linux ,no WXP , no C+ , no PC .....

 

[Ben Finney]

Browsers . There will be 2 columns , one on left will be main/orig
but on the Right will have hyperlink result . This way ya dont have
to go back to compare ! side by side .
Text editors will also work this way . You will read orig in left
columns
and 'edit'/annotate your stuff to the Right columns .

Hmm. Your ideas are intriguing to me and I wish to subscribe to your
newsletter.

 

[Paul Rubin]

Ok, I'll bite. How do you compute the median of a list using just a single
temp var?

Well there's an obvious quadratic-time method...

 

[bearophileHUGS]

no sort() is needed to calculate the median of a list.
you just need one temp var.

Can you show some actual code?

(There is the median of 5 algorithm too).

Bye,
bearophile

 

[Fredrik Lundh]

Well there's an obvious quadratic-time method...

that does it without modifying the list?

if you can modify the list, there are plenty of algorithms that does it
in expected O(n) or better, but I cannot think of a one that doesn't use
at least a few variables (e.g. two list indexes and a pivot).

but I haven't had enough coffee yet, so I'm probably missing something
simple here.

</F>

 

[Paul Rubin]

that does it without modifying the list?

if you can modify the list, there are plenty of algorithms that does
it in expected O(n) or better, but I cannot think of a one that
doesn't use at least a few variables (e.g. two list indexes and a pivot).

Hmm, whoops, I didn't count the list index for the quadratic time
version (but that version shouldn't need to modify the list).

If you can modify the list, let's see, you can swap two elements with
no temp vars:

a ^= a[i+1]
a[i+1] ^= a
a ^= a[i+1]

This assumes an indexed addressing mode so finding a[i+1] doesn't
require using a temp var to hold (i+1). Let's say the list length is
n, which is not a variable, and constant expressions like n-1 are also
not variables. I'm still envisioning some reasonable type of assembly
code. So now we can straightforwardly sort the list with one index
var and one flag bit:

flag = True
while flag:
flag = False
for i in 0..(n-2):
if a > a[i+1]:
swap a, a[i+1] as above
flag = True

and then pick the median out of the middle.

but I haven't had enough coffee yet, so I'm probably missing something
simple here.

Yeah, it's night here, maybe after I get some sleep I'll look for a
way to get rid of the flag bit above.

 

[Ian McConnell]

John Doty:

Yes. The efficient exact algorithms for this problem use *partial*
sorts. The Forth one from the FSL is of this class (although I know of
two better ones for big arrays). But it's tough to beat the efficiency
of the approximate histogram-based method the Python stats module
implements.

The usual way to compute a true median with Python may be:

def median(inlist):
newlist = sorted(inlist)
index = len(newlist) // 2
if len(newlist) % 2:
return newlist[index]
else:
return (newlist[index] + newlist[index-1]) / 2.0

If you can use Psyco and your FITS lines are really long (well, maybe
too much, the treshold if about >~3000 in my PC) you can use something
like this instead the builtin timsort:
http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/466330
(To compute the median on a image line the median that uses sort is
probably better in most cases, expecially if you use the built in sort
of numerical libraries.)

sort() sorts all of the data, but you're only after one or two numbers, so
the MODFIND method may be faster for the median:

http://www.geocities.com/SiliconValley/Garage/3323/algor.html

Ian

 

[idknow]

[snip]

no sort() is needed to calculate the median of a list.

you just need one temp var.

Ok, I'll bite. How do you compute the median of a list using just a single
temp var?

-- Paul

hi Paul; well when this was a stats-class assignment (back when pascal
was popular i just stepped through the vector and compared it

(pseudo-code)

ptr p = [with values].

fun median {
var x = 0.
while( *p++) {
if( (*p) > x) x = *p.
}
return x.
}

of course, pascal is more verbose but that's median()

 

[Paul Rubin]

fun median {
var x = 0.
while( *p++) {
if( (*p) > x) x = *p.
}
return x.
}

I count two variables, p and x.

 

[Paul Rubin]

I count two variables, p and x.

Also, that finds the maximum, not the median. I had stopped examining
it after seeing it used more than one variable.

 

[Andrew Haley]

I count two variables, p and x.

Isn't this the maximum?

Andrew.

 

[Fredrik Lundh]

hi Paul; well when this was a stats-class assignment (back when pascal
was popular i just stepped through the vector and compared it

(pseudo-code)

ptr p = [with values].

fun median {
var x = 0.
while( *p++) {
if( (*p) > x) x = *p.
}
return x.
}

of course, pascal is more verbose but that's median()

that's a rather unusual definition of median, though.

</F>