oAxiom - Biotech,

Reimplemntation of specrend.c in Python

2011-04-08T07:01:00.000-07:00

Reimplementation of Specrend.c in Python.

This, (broadly) converts Temperature into and r, g, b value. Or, it gets you close enough that the rest is trivial and is left as an exercise to the reader.

The original is taken from:

specrend.c

And the original website is:

www.fourmilab.ch

Here's some Python:

"""
                Colour Rendering of Spectra

                       by John Walker
                  http://www.fourmilab.ch/
          
         Last updated: March 9, 2003

    Converted to Python by Andrew Hutchins, sometime in early
    2011.

           This program is in the public domain.
           The modifications are also public domain. (AH)

    For complete information about the techniques employed in
    this program, see the World-Wide Web document:

             http://www.fourmilab.ch/documents/specrend/
         
    The xyz_to_rgb() function, which was wrong in the original
    version of this program, was corrected by:
    
        Andrew J. S. Hamilton 21 May 1999
        Andrew.Hamilton@Colorado.EDU
        http://casa.colorado.edu/~ajsh/

    who also added the gamma correction facilities and
    modified constrain_rgb() to work by desaturating the
    colour by adding white.
    
    A program which uses these functions to plot CIE
    "tongue" diagrams called "ppmcie" is included in
    the Netpbm graphics toolkit:
        http://netpbm.sourceforge.net/
    (The program was called cietoppm in earlier
    versions of Netpbm.)

"""
import math

"""
/* A colour system is defined by the CIE x and y coordinates of
   its three primary illuminants and the x and y coordinates of
   the white point. */
"""
GAMMA_REC709 = 0

NTSCsystem  =  {"name": "NTSC", 
    "xRed": 0.67, "yRed": 0.33, 
    "xGreen": 0.21, "yGreen": 0.71, 
    "xBlue": 0.14, "yBlue": 0.08, 
    "xWhite": 0.3101, "yWhite": 0.3163, "gamma": GAMMA_REC709}
EBUsystem  =  {"name": "SUBU (PAL/SECAM)", 
    "xRed": 0.64, "yRed": 0.33, 
    "xGreen": 0.29, "yGreen": 0.60, 
    "xBlue": 0.15, "yBlue": 0.06, 
    "xWhite": 0.3127, "yWhite": 0.3291, "gamma": GAMMA_REC709 }
SMPTEsystem  =  {"name": "SMPTE", 
    "xRed": 0.63, "yRed": 0.34, 
    "xGreen": 0.31, "yGreen": 0.595, 
    "xBlue": 0.155, "yBlue": 0.07, 
    "xWhite": 0.3127, "yWhite": 0.3291, "gamma": GAMMA_REC709 }
HDTVsystem  =  {"name": "HDTV", 
    "xRed": 0.67, "yRed": 0.33, 
    "xGreen": 0.21, "yGreen": 0.71, 
    "xBlue": 0.15, "yBlue": 0.06, 
    "xWhite": 0.3127, "yWhite": 0.3291, "gamma": GAMMA_REC709 }
CIEsystem  =  {"name": "CIE", 
    "xRed": 0.7355, "yRed": 0.2645, 
    "xGreen": 0.2658, "yGreen": 0.7243, 
    "xBlue": 0.1669, "yBlue": 0.0085, 
    "xWhite": 0.3333333333, "yWhite": 0.3333333333, "gamma": GAMMA_REC709 }
Rec709system  =  {"name": "CIE REC709", 
    "xRed": 0.64, "yRed": 0.33, 
    "xGreen": 0.30, "yGreen": 0.60, 
    "xBlue": 0.15, "yBlue": 0.06, 
    "xWhite": 0.3127, "yWhite": 0.3291, "gamma": GAMMA_REC709 }

def upvp_to_xy(up, vp):
    xc = (9 * up) / ((6 * up) - (16 * vp) + 12)
    yc = (4 * vp) / ((6 * up) - (16 * vp) + 12)
    return(xc, yc)

def xy_toupvp(xc, yc):
    up = (4 * xc) / ((-2 * xc) + (12 * yc) + 3);
    vp = (9 * yc) / ((-2 * xc) + (12 * yc) + 3);  
    return(up, vp)

def xyz_to_rgb(cs, xc, yc, zc):
    """
    Given an additive tricolour system CS, defined by the CIE x
    and y chromaticities of its three primaries (z is derived
    trivially as 1-(x+y)), and a desired chromaticity (XC, YC,
    ZC) in CIE space, determine the contribution of each
    primary in a linear combination which sums to the desired
    chromaticity.  If the  requested chromaticity falls outside
    the Maxwell  triangle (colour gamut) formed by the three
    primaries, one of the r, g, or b weights will be negative. 

    Caller can use constrain_rgb() to desaturate an
    outside-gamut colour to the closest representation within
    the available gamut and/or norm_rgb to normalise the RGB
    components so the largest nonzero component has value 1.
    """

    xr = cs["xRed"]
    yr = cs["yRed"]
    zr = 1 - (xr + yr)
    xg = cs["xGreen"]
    yg = cs["yGreen"]
    zg = 1 - (xg + yg)
    xb = cs["xBlue"]
    yb = cs["yBlue"]
    zb = 1 - (xb + yb)
    xw = cs["xWhite"]
    yw = cs["yWhite"]
    zw = 1 - (xw + yw)
    
    rx = (yg * zb) - (yb * zg)  
    ry = (xb * zg) - (xg * zb)  
    rz = (xg * yb) - (xb * yg)
    gx = (yb * zr) - (yr * zb)  
    gy = (xr * zb) - (xb * zr)  
    gz = (xb * yr) - (xr * yb)
    bx = (yr * zg) - (yg * zr)  
    by = (xg * zr) - (xr * zg)  
    bz = (xr * yg) - (xg * yr)
       
    rw = ((rx * xw) + (ry * yw) + (rz * zw)) / yw
    gw = ((gx * xw) + (gy * yw) + (gz * zw)) / yw
    bw = ((bx * xw) + (by * yw) + (bz * zw)) / yw
    
    rx = rx / rw;  ry = ry / rw;  rz = rz / rw
    gx = gx / gw;  gy = gy / gw;  gz = gz / gw
    bx = bx / bw;  by = by / bw;  bz = bz / bw
    
    r = (rx * xc) + (ry * yc) + (rz * zc)
    g = (gx * xc) + (gy * yc) + (gz * zc)
    b = (bx * xc) + (by * yc) + (bz * zc)
    return(r,g,b)


def inside_gamut(r, g, b):
    """ 
     Test whether a requested colour is within the gamut
     achievable with the primaries of the current colour
     system.  This amounts simply to testing whether all the
     primary weights are non-negative. */
    """    
    return (r >= 0) and (g >= 0) and (b >= 0)


def constrain_rgb(r, g, b):
    """
    If the requested RGB shade contains a negative weight for
    one of the primaries, it lies outside the colour gamut 
    accessible from the given triple of primaries.  Desaturate
    it by adding white, equal quantities of R, G, and B, enough
    to make RGB all positive.  The function returns 1 if the
    components were modified, zero otherwise.
    """
    # Amount of white needed is w = - min(0, *r, *g, *b)
    w = -min([0, r, g, b]) # I think?

    # Add just enough white to make r, g, b all positive.
    if w > 0:
        r += w
        g += w
        b += w
    return(r,g,b)

def gamma_correct(cs, c):
    """
    Transform linear RGB values to nonlinear RGB values. Rec.
    709 is ITU-R Recommendation BT. 709 (1990) ``Basic
    Parameter Values for the HDTV Standard for the Studio and
    for International Programme Exchange'', formerly CCIR Rec.
    709. For details see
    
       http://www.poynton.com/ColorFAQ.html
       http://www.poynton.com/GammaFAQ.html
    """
    gamma = cs.gamma
    
    if gamma == GAMMA_REC709:
        cc = 0.018
        if c < cc:
            c = ((1.099 * math.pow(cc, 0.45)) - 0.099) / cc
        else:
            c = (1.099 * math.pow(c, 0.45)) - 0.099
    else:
        c = math.pow(c, 1.0 / gamma)
    return(c)

def gamma_correct_rgb(cs, r, g, b):
    r = gamma_correct(cs, r)
    g = gamma_correct(cs, g)
    b = gamma_correct(cs, b)
    return(r,g,b)

def norm_rgb(r, g, b):
    """
    Normalise RGB components so the most intense (unless all
    are zero) has a value of 1.
    """
    greatest = max([r, g, b])
    
    if greatest > 0:
        r /= greatest
        g /= greatest
        b /= greatest
    return(r, g, b)
    
def spectrum_to_xyz(spec_intens, temp): #spec_intens is a function
    """
    Calculate the CIE X, Y, and Z coordinates corresponding to
    a light source with spectral distribution given by  the
    function SPEC_INTENS, which is called with a series of
    wavelengths between 380 and 780 nm (the argument is 
    expressed in meters), which returns emittance at  that
    wavelength in arbitrary units.  The chromaticity
    coordinates of the spectrum are returned in the x, y, and z
    arguments which respect the identity:

            x + y + z = 1.
    
    CIE colour matching functions xBar, yBar, and zBar for
       wavelengths from 380 through 780 nanometers, every 5
       nanometers.  For a wavelength lambda in this range:

            cie_colour_match[(lambda - 380) / 5][0] = xBar
            cie_colour_match[(lambda - 380) / 5][1] = yBar
            cie_colour_match[(lambda - 380) / 5][2] = zBar

    AH Note 2011: This next bit is kind of irrelevant on modern 
    hardware. Unless you are desperate for speed.
    In which case don't use the Python version!
    
    To save memory, this table can be declared as floats
    rather than doubles; (IEEE) float has enough 
    significant bits to represent the values. It's declared
    as a double here to avoid warnings about "conversion
    between floating-point types" from certain persnickety
    compilers. */
    """

    cie_colour_match = [
        [0.0014,0.0000,0.0065], [0.0022,0.0001,0.0105], [0.0042,0.0001,0.0201],
        [0.0076,0.0002,0.0362], [0.0143,0.0004,0.0679], [0.0232,0.0006,0.1102],
        [0.0435,0.0012,0.2074], [0.0776,0.0022,0.3713], [0.1344,0.0040,0.6456],
        [0.2148,0.0073,1.0391], [0.2839,0.0116,1.3856], [0.3285,0.0168,1.6230],
        [0.3483,0.0230,1.7471], [0.3481,0.0298,1.7826], [0.3362,0.0380,1.7721],
        [0.3187,0.0480,1.7441], [0.2908,0.0600,1.6692], [0.2511,0.0739,1.5281],
        [0.1954,0.0910,1.2876], [0.1421,0.1126,1.0419], [0.0956,0.1390,0.8130],
        [0.0580,0.1693,0.6162], [0.0320,0.2080,0.4652], [0.0147,0.2586,0.3533],
        [0.0049,0.3230,0.2720], [0.0024,0.4073,0.2123], [0.0093,0.5030,0.1582],
        [0.0291,0.6082,0.1117], [0.0633,0.7100,0.0782], [0.1096,0.7932,0.0573],
        [0.1655,0.8620,0.0422], [0.2257,0.9149,0.0298], [0.2904,0.9540,0.0203],
        [0.3597,0.9803,0.0134], [0.4334,0.9950,0.0087], [0.5121,1.0000,0.0057],
        [0.5945,0.9950,0.0039], [0.6784,0.9786,0.0027], [0.7621,0.9520,0.0021],
        [0.8425,0.9154,0.0018], [0.9163,0.8700,0.0017], [0.9786,0.8163,0.0014],
        [1.0263,0.7570,0.0011], [1.0567,0.6949,0.0010], [1.0622,0.6310,0.0008],
        [1.0456,0.5668,0.0006], [1.0026,0.5030,0.0003], [0.9384,0.4412,0.0002],
        [0.8544,0.3810,0.0002], [0.7514,0.3210,0.0001], [0.6424,0.2650,0.0000],
        [0.5419,0.2170,0.0000], [0.4479,0.1750,0.0000], [0.3608,0.1382,0.0000],
        [0.2835,0.1070,0.0000], [0.2187,0.0816,0.0000], [0.1649,0.0610,0.0000],
        [0.1212,0.0446,0.0000], [0.0874,0.0320,0.0000], [0.0636,0.0232,0.0000],
        [0.0468,0.0170,0.0000], [0.0329,0.0119,0.0000], [0.0227,0.0082,0.0000],
        [0.0158,0.0057,0.0000], [0.0114,0.0041,0.0000], [0.0081,0.0029,0.0000],
        [0.0058,0.0021,0.0000], [0.0041,0.0015,0.0000], [0.0029,0.0010,0.0000],
        [0.0020,0.0007,0.0000], [0.0014,0.0005,0.0000], [0.0010,0.0004,0.0000],
        [0.0007,0.0002,0.0000], [0.0005,0.0002,0.0000], [0.0003,0.0001,0.0000],
        [0.0002,0.0001,0.0000], [0.0002,0.0001,0.0000], [0.0001,0.0000,0.0000],
        [0.0001,0.0000,0.0000], [0.0001,0.0000,0.0000], [0.0000,0.0000,0.0000]]

    X = 0
    Y = 0
    Z = 0
    for i, lamb in enumerate(range(380, 780, 5)): #lambda = 380; lambda < 780.1; i++, lambda += 5) {
        Me = spec_intens(lamb, temp);
        X += Me * cie_colour_match[i][0]
        Y += Me * cie_colour_match[i][1]
        Z += Me * cie_colour_match[i][2]
    XYZ = (X + Y + Z)
    x = X / XYZ;
    y = Y / XYZ;
    z = Z / XYZ;
    return(x, y, z)

def bb_spectrum(wavelength, bbTemp=5000):
    """
    Calculate, by Planck's radiation law, the emittance of a black body
    of temperature bbTemp at the given wavelength (in metres).  */
    """
    wlm = wavelength * 1e-9 # Convert to metres 
    return (3.74183e-16 * math.pow(wlm, -5.0)) / (math.exp(1.4388e-2 / (wlm * bbTemp)) - 1.0)

    """  Built-in test program which displays the x, y, and Z and RGB
    values for black body spectra from 1000 to 10000 degrees kelvin.
    When run, this program should produce the following output:

    Temperature       x      y      z       R     G     B
    -----------    ------ ------ ------   ----- ----- -----
       1000 K      0.6528 0.3444 0.0028   1.000 0.007 0.000 (Approximation)
       1500 K      0.5857 0.3931 0.0212   1.000 0.126 0.000 (Approximation)
       2000 K      0.5267 0.4133 0.0600   1.000 0.234 0.010
       2500 K      0.4770 0.4137 0.1093   1.000 0.349 0.067
       3000 K      0.4369 0.4041 0.1590   1.000 0.454 0.151
       3500 K      0.4053 0.3907 0.2040   1.000 0.549 0.254
       4000 K      0.3805 0.3768 0.2428   1.000 0.635 0.370
       4500 K      0.3608 0.3636 0.2756   1.000 0.710 0.493
       5000 K      0.3451 0.3516 0.3032   1.000 0.778 0.620
       5500 K      0.3325 0.3411 0.3265   1.000 0.837 0.746
       6000 K      0.3221 0.3318 0.3461   1.000 0.890 0.869
       6500 K      0.3135 0.3237 0.3628   1.000 0.937 0.988
       7000 K      0.3064 0.3166 0.3770   0.907 0.888 1.000
       7500 K      0.3004 0.3103 0.3893   0.827 0.839 1.000
       8000 K      0.2952 0.3048 0.4000   0.762 0.800 1.000
       8500 K      0.2908 0.3000 0.4093   0.711 0.766 1.000
       9000 K      0.2869 0.2956 0.4174   0.668 0.738 1.000
       9500 K      0.2836 0.2918 0.4246   0.632 0.714 1.000
      10000 K      0.2807 0.2884 0.4310   0.602 0.693 1.000
"""

if __name__ == "__main__":
    print "Temperature       x      y      z       R     G     B\n"
    print "-----------    ------ ------ ------   ----- ----- -----\n"
    
    for t in range(1000, 10000, 500):  # (t = 1000; t <= 10000; t+= 500) {
        x, y, z = spectrum_to_xyz(bb_spectrum, t)
        
        r, g, b = xyz_to_rgb(SMPTEsystem, x, y, z)
        
        print "  %5.0f K      %.4f %.4f %.4f   " % (t, x, y, z),
        
        r, g, b = constrain_rgb(r, g, b) # I omit the approximation bit here.
        r, g, b = norm_rgb(r, g, b) 
        print "%.3f %.3f %.3f" % (r, g, b)

Drip drip drip

2009-12-16T04:25:00.001-08:00

Drip drip drip, originally uploaded by afflictedmonkey.

Saline bag in the hospital. Not mine luckily.

By The Shipyard

2009-12-16T04:19:00.001-08:00

By The Shipyard, originally uploaded by afflictedmonkey.

Perhaps the last pictures of the docks didn't do full justice. This is the docks/ship fit out works on West Coast park in Singapore. It's a neat place, and I was lucky to be there when the sun went down.

Singapore's natural resources.

2009-09-27T00:06:00.000-07:00

It's a common refrain in Singapore that the country has no natural resources.

This is not true, and has never been true. Singapore has several natural resources, but the one that occupies prime position is the docks. Singapore has one of the finest deep-water docks in the world. Often taking the busiest port in the world title.

Google is fast

2009-09-18T06:17:00.001-07:00

45 minutes to index and report. Possibly less - that's just the first time I checked.

(See the previous entry for details)

Chinese Country Names

2009-09-18T05:31:00.000-07:00

Chinese, of course, doesn't have an alphabet, instead it uses ideograms each one unique and represents some sort of meaning. There are several tens of thousands of characters. Although for every day speech and to be considered fluent you need to be conversant in about 3000 characters (which seems some sort of magic number, about 3000 English words gets you to fluent, 5000 you are learned and educated whereas 10,000 and you are a smart arse).

But, this poses something of a problem. In English when you need a word you just stitch together a bunch of Latin or root words or, hell, you just make it up (gruzbumble - there, a new word for you, 0 hits in google, it means "a bumble bee who's really pissed off"). But what do you do in Chinese? Well, you can actually transliterate - converting the sound, from say, English, into a set of Chinese characters. This works quite okay, although the pronunciation can be quite a bit off it's usually recognisable.

However, that new transliterated word will be made up of characters that are separate from the new word - they have their own inherent meaning.

So (and this is a round about way of getting to my point), for example, when China makes up a name for a country, what do they do? Well they transliterate the sound, but they inadvertently (and some times, uh, vertantly) leave behind a meaning in the characters, for example:

英国, pronounced 'ying guo' means literally 'hero country' = United Kingdom
美国 = meiguo = beautiful country = USA.
法国 = Fa guo = Law Country = France
泰国 = Tai Guo = Lucky Country = Thailand

But some times ironic or silly meanings can be left behind, one example:

德国 = De Guo = Moral Country = Germany
日本 = Ri ben = Sunny book = Japan
黑山 = Hei Shan = Black Mountain = Montenegro

But then sometimes, by accident really silly stuff gets left behind in the characters:

匈牙利 = Xiong ya li = Hungary = Fearsome Teeth Are sharp!
柬埔寨 = Jian Pu Jia = Cambodia = Crappy Village (Okay, this is a pun, actually)
马来西亚 = Ma Lai Xi Ya = Malaysia = Horse comes to West Asia (?!?!)
葡萄牙 = Pu Tao Ya = Portugal = Grape Teeth (?!?!!)
马里 = ma li = Mali = Horse Inside (?)
索马里 = Suo Ma li = Somalia = Search for horse Inside
多哥 = Duo Ge = Togo = Many brothers

Heh, sometimes it also goes wrong with other things. My home town (Northampton) is transliterated as Bei am dun, the 'bei' means literally 'North'. There are two major problems with this name - Northampton is not in the North, it's in the middle, (the midlands, to be exact, the home of rugby, but little else). And crucially, the 'North' in Northampton doesn't mean 'North'. It means 'Norse' as in Norse hampton = Norse town = People from Scandinavia live here in a town.

It only stops when you get to Bangkok

2009-09-14T04:36:00.001-07:00

It only stops when you get to Bangkok, originally uploaded by afflictedmonkey.

Railway line near Buona Vista. It's technically part of Malaysia so I suppose I was trespassing.

There's a ruined building just to the right of the picture, and occasionally trains traveling at about 20 mph go by. Embarrassingly I nearly got run over by one once! (It's the fastest way to Colbar)

Cat on the Path

2009-09-09T06:27:00.001-07:00

Cat on the Path, originally uploaded by afflictedmonkey.

There's a load of cats at the condo here, they tend to hang out around the base of our block waiting for food. Each cat has its own territory. One day I should map it out and see which cat owns which part of the condo.

Our block belongs to the white cat relaxing near the path light.

Fingers and Claws

2009-08-22T06:17:00.003-07:00

Fingers and Claws
Originally uploaded by afflictedmonkey

Transit Hub

2009-08-22T06:17:00.001-07:00

Transit Hub
Originally uploaded by afflictedmonkey

Playing with my camera, a few shots from flickr.

To War!

2009-08-22T06:16:00.001-07:00

To War!
Originally uploaded by afflictedmonkey

At the Turtle museum in Chinese Garden, they all live free in the pool and go crazy when people feed them.

On the folly of investors advice.

2009-08-15T07:21:00.000-07:00

Teh Hooi Ling in this mornings Business Times reiterates many of the common investment advice for novice investors (like myself). I guess, this morning my mind finally snapped and I just couldn't take the mindless thoughtless advice constantly repeated as if it is mantra from these professionals. Actually, I'm being a little unfair in singling out Teh Hooi Ling, who's column I often enjoy, but to distill it to it's base, investment advice commonly boils down to this simple set of vague pointers:

1. Compound interest is miraculous, at 8-10% you will be a millionaire in your 50's!
2. Research and think in detail about the stocks you buy,
3. Sell when things are too risky!

The first one really pisses me off, but I'll come back to that one, so I'll deal with objection 1 first: I think a lot of this sort of advice is really lacking in specifics - you know really practical specifics, e.g. buy a stock when the PE goes below 10, or buy the market after 3 straight days of losses. I can understand this totally - there are no specifics to give. Every rigid investment strategy known to man has at some time found wanting (value investing, dividend yield, technical, fundamental etc.), so there are no specifics to give. Then, you base your research on a companies last quarter of results (assuming it's available, and for the Singapore Exchange it often isn't). But that is in the past, and the exchange prices the share *on future earnings*. So, what, at the end of the day is research - it all comes down to a gut feeling on whether a share is good or not.

Objection 3 is easy to dismiss - risk is obvious in hindsight. I know a lot of people who got out of the market in 2005 as it appeared too risky. It was, for sure, but in the meantime it continued to run up until 2008, only then did it crash. Risk and value are easy in hindsight. Examples from the past are meaningless. Currently as I write this I think that risk has possibly never been higher in the world, due to a variety of factors, essentially between China and the US, but perhaps in another entry.

Objection 1 is the one that always gets me though. Not because it is untrue, it is true - 100,000 units of currency invested at 8-10 % will become 1,000,000 units in about 30 years. That's maths. My objection is - where can I find and purchase these investment vehicles that pay me 8-10 % per year with little to no risk? Or, perhaps even with a lot of risk? The local banks pay a shocking 0.125% interest, whilst their 'special' 'super' 'amazing' accounts pay a whopping 0.8%! Oh. my. god. Hold the front page! The Straits Time Index (an index of approx 70% of the turnover on the singapore exchange) pays a dividend yield of around 3.7% (as of 2009, it's probably going to be lower in the near future due to reduced profitability).

I suppose you could buy preference shares, DBS will give you 6% and UOB 5%, or you could get into some of the REITs, which do start to get near the magical 8-10% (Suntec Reit pays ~7%, although it's share price is rising). And finally you could buy some of the really wierd esoteric dividend yield vehicles - like the shipping trusts, some of which pay ~12%. But you would be a fool to put your trust in such horrifyingly volatile instruments, they pay 12% or more for a reason (er... because sometimes, they don't pay anything at all).

Then there is the old chestnut that the stock exchange grows at 10% on average per year. But these measurements always seem to choose good windows for measurement. Let's look at the STI from 1987 to today. In 1987 it was ~800. Today it hovers around 2600 after a pretty robust rebound (despite Singapore still being in a recession). That's a gain of 69%. Or, 3.13% per year (actually it's less than that due to compounding, which I can't be bothered to work out). So, adding in a generous dividend yield of 3.7% (it has not always been so high) and the generous growth of 3.13%, that still only makes 6.9%. Still a couple of percents short of the magical 8-10%. Actually there is a reason it's 6.9%, and it's terribly simple - when all is said and done, the exchange can go up and down (and it sure can go up and down a hell of a lot), but it must return to its trend line - GDP growth.

Well, that's probably enough for now. I haven't talked about bonds (appalling, even risky companies can get away with ~4% and a lot of the juicy stuff is not available to the retail investor - you need ~100,000 SGD to even get started), inflation linked treasuries (miserly, inflation +2-3%) hedge funds (not always available to the retail investor, and anyway fund of hedge funds only get the same as the market over the long term).

Rant over.

2009-05-27T18:41:00.001-07:00

GLbase

Analysis of ChIP-seq data.

From the preamble:
This is GLbase, a set of tools to deal with heterogeneous genomic data. It's particularly useful for intersecting any data that can be expressed as a CSV file. Examples are microarray data, gene lists, Chip-seq location lists. GLbase excels at integrating these data and then outputting graphical representations of the resulting data:

You can get GLbase here:

Some screenshots of the output it can generate (via R):

Google Earth

2009-05-04T17:11:00.000-07:00

Google Earth is Fun

2004-2009 - Singapore

2000-2004 - Norwich, UK.

Renvy 0.1

2009-04-26T03:42:00.000-07:00

Renvy 0.1, in its most basic form is usable. It will load and save R scripts, and you can run them. Plenty of things to do, it only deals with single tabs, no cut copy paste, no undo redo, no syntax highlighting or other nice bells and whistles.

Renvy requires rpy, python, and R. Click main.py to run Renvy. Or type:

$ python main.py

to run Renvy. You can get Renvy here. Or clone the mercurial repository with:

$ hg clone https://oaxiom@bitbucket.org/oaxiom/renvy/

Renvy

2009-04-19T18:33:00.000-07:00

My general disappointment with the available R GUI's has finally led me to develop my own, Renvy, a simple GUI front-end to R without much of the bells and whistles. Think of Geany for coding, Renvy is to R. [Actually Geany can integrate R, it just does it poorly].

R is a statistical programming language that has gained a lot of traction amongst biologists due to the excellent Bioconductor package which allows researchers to do amazing amounts of analysis with the minimum of fuss.

You can follow the progress of Renvy here. Suffice to say it is not usable currently!

Economics - China and the reserve currency.

2009-04-08T17:52:00.000-07:00

Hmm... Probably not the best entry to start with, 200 lines of Python code. It's just I'd only seen this code as fragments on the web and had never seen it pieced together and integrated into a GUI.

I commented on a blog entry here http://www.riseupforthenewday.net/2009/03/24/china-moves-to-position-itself-as-a-power-player-in-time-of-crisis/comment-page-1/#comment-36 talking about China's push for a reserve currency.

I've edited it a little for clarity.

I’m not really sure the Chinese could really see this financial crisis coming. If they could then they certainly would not have parked massive reserves in the dollar nor their disastrous investment in Morgan Stanley and Blackstone.

As for the crisis not affecting China’s economy - I think it is to early to declare that. There is a large lag between import/export declines and direct effects on the economy. World-wide logistic chains take a while to unravel.

6.5% growth forecast in China next year sounds impressive, but is actually considered a catastrophe in China (re: protect 8, ie preserve 8% growth), this forecast is likely to be revised downwards as well, unlikely to go negative, but in a country of 1.2bn people with a dubious social contract with the people (don’t do anything political and we’ll provide the growth) it is a recipe for combustion. If China really does escape this crisis unscathed, then it is likely just saving up problems in their banking system for the future.

As for the new reserve currency - this is just China grandstanding to get a place at the table. The chances of it getting off the ground are zero. And commodity linked currencies would be a disaster, plus whom would guarantee the currency? The USD is the reserve currency because it is the largest economy with a massive and liquid trade in dollars. If China was serious about a new reserve currency then it should concentrate on opening up its own currency - as the third largest economy in the world it could form a bipolar reserve currency with the US. But it would not do that because the China government (rightly) fears massive currency movements overseas and is still chastened by 1997.

Integrating Cairo into wxWidgets...

2009-03-27T08:02:00.000-07:00

Some python code to integrate Cairo into wxWidgets:

This goes into the top of your code:


#----------------------------------------------------------------------
# The Cairo interface.
import wx.lib.wxcairo

import ctypes
import cairo
from ctypes.util import find_library

cairo_dll = ctypes.CDLL(find_library("cairo"))

# Pycairo's API representation (from pycairo.h)
class Pycairo_CAPI(ctypes.Structure):
  _fields_ = [
     ('Context_Type', ctypes.py_object),
     ('Context_FromContext', ctypes.PYFUNCTYPE(ctypes.py_object,
                                               ctypes.c_void_p,
                                               ctypes.py_object,
                                               ctypes.py_object)),
     ('FontFace_Type', ctypes.py_object),
     ('FontFace_FromFontFace', ctypes.PYFUNCTYPE(ctypes.py_object, ctypes.c_void_p)),
     ('FontOptions_Type', ctypes.py_object),
     ('FontOptions_FromFontOptions', ctypes.PYFUNCTYPE(ctypes.py_object, ctypes.c_void_p)),
     ('Matrix_Type', ctypes.py_object),
     ('Matrix_FromMatrix', ctypes.PYFUNCTYPE(ctypes.py_object, ctypes.c_void_p)),
     ('Path_Type', ctypes.py_object),
     ('Path_FromPath', ctypes.PYFUNCTYPE(ctypes.py_object, ctypes.c_void_p)),
     ('Pattern_Type', ctypes.py_object),
     ('SolidPattern_Type', ctypes.py_object),
     ('SurfacePattern_Type', ctypes.py_object),
     ('Gradient_Type', ctypes.py_object),
     ('LinearGradient_Type', ctypes.py_object),
     ('RadialGradient_Type', ctypes.py_object),
     ('Pattern_FromPattern', ctypes.c_void_p),
     ('ScaledFont_Type', ctypes.py_object),
     ('ScaledFont_FromScaledFont', ctypes.PYFUNCTYPE(ctypes.py_object, ctypes.c_void_p)),
     ('Surface_Type', ctypes.py_object),
     ('ImageSurface_Type', ctypes.py_object),
     ('PDFSurface_Type', ctypes.py_object),
     ('PSSurface_Type', ctypes.py_object),
     ('SVGSurface_Type', ctypes.py_object),
     ('Win32Surface_Type', ctypes.py_object),
     ('XlibSurface_Type', ctypes.py_object),
     ('Surface_FromSurface', ctypes.PYFUNCTYPE(ctypes.py_object, ctypes.c_void_p)),
     ('Check_Status', ctypes.PYFUNCTYPE(ctypes.c_int, ctypes.c_int))]

# look up the API
ctypes.pythonapi.PyCObject_Import.restype = ctypes.POINTER(Pycairo_CAPI)
pycairo_api = ctypes.pythonapi.PyCObject_Import("cairo", "CAPI").contents;

ContextType = pycairo_api.Context_Type

def Context_FromSWIGObject(swigObj):
   """
   (Internal)
   get the Cairo object for drawing.
   """
   ptr = ctypes.c_void_p(int(swigObj))
   #increment the native context's ref count, since the Pycairo_Context decrements it
   #when it is finalised.
   cairo_dll.cairo_reference(ptr)
   return pycairo_api.Context_FromContext(ptr, ContextType, None)

Then inside your drawing class you will need to add an OnPaint event:


    def OnPaint(self, event):
        """
        Call to get Cairo to paint, this is done automatically if bound
        to a Panel.
        """
        try:
            dc = wx.PaintDC(self.panel) # make each time OnPaint is called.
            #dc = wx.AutoBufferedPaintDC(self) # not clear why this doesn't work ...
            self.size = dc.GetSizeTuple()
            self.setViewPortSize(self.size[0], self.size[1])
            gc = wx.GraphicsContext.Create(dc)
            nc = gc.GetNativeContext()
            ctx = Context_FromSWIGObject(nc)
        except:
            raise ErrorCairoAcquireDevice

        self._paint(ctx) # actual paint calls go here.
        return(True)

I found it useful to bind all of this into a panel on my wxWidgets GUI. I have a panel in the GUI already, then I bind another panel on top of that panel.
My panel needs a few special options (see below), so I can't just bind it into that panel.


    def bindPanel(self, panel):
        """
        bind a wx.Panel into gDraw.
        """
        self.panel = drawPanel(panel, self.OnPaint)

        return(self.panel)

And this is my drawPanel class, which I used to bind the events into and exists somewhere were I need it in the GUI.


class drawPanel(wx.Panel):
    """
    An empty Panel class that sets things up suitably for Cairo.
    """
    def __init__(self, parent, paint_Procedure):
        # This has to be set to full repaint to work, 
        # otherwise it will try to repaint only small parts.
        wx.Panel.__init__(self, parent, -1, size=(-1,-1), style=wx.FULL_REPAINT_ON_RESIZE) 
        self.Bind(wx.EVT_PAINT, paint_Procedure, self)
        self.Fit()

Welcome to oAxiom

2009-03-27T07:54:00.000-07:00

An occasional spot for me to write stuff I find interesting, and my faltering efforts to find funding for a biotech venture.

I sit at the interface of biology and bioinformatics - I'm in the unique position to be both a fully trained biologist (cell culture, Western blots, all that stuff). And I am also an experienced programmer (Python, C, R).

My project sits at the interface between these two people. As such this blog will contain snippets of code, snippets of biology and a dose of entrepreneurship, as wells as other miscellaneous nonsense.

Welcome to oAxiom, its only employee (if I could call myself that): Andrew Hutchins. I'm in Singapore (currently), and you can contact me at oaxiom@oaxiom.com.