| Andrew Cooke | Contents | Latest | RSS | Twitter | Previous | Next

C[omp]ute

Welcome to my blog, which was once a mailing list of the same name and is still generated by mail. Please reply via the "comment" links.

Always interested in offers/projects/new ideas. Eclectic experience in fields like: numerical computing; Python web; Java enterprise; functional languages; GPGPU; SQL databases; etc. Based in Santiago, Chile; telecommute worldwide. CV; email.

Personal Projects

Lepl parser for Python.

Colorless Green.

Photography around Santiago.

SVG experiment.

Professional Portfolio

Calibration of seismometers.

Data access via web services.

Cache rewrite.

Extending OpenSSH.

Last 100 entries

Programming, Business, Social Contracts; Distributed Systems for Fun and Profit; XML and Scheme; Internet Radio Stations (Curated List); Solid Data About Placebos; Half of Americans Think Climate Change Is a Sign of the Apocalypse; Saturday Surf Sessions With Juvenile Delinquents; Ssh, tty, stdout and stderr; Feathers falling in a vacuum; Santiago 30m Bike Route; Mapa de Ciclovias en Santiago; How Unreliable is UDP?; SE Santiago 20m Bike Route; Cameron's Rap; Configuring libxml with Eclipse; Reducing Combinatorial Complexity With Occam - AI; Sentidos Comunes (Chilean Online Magazine); Hilary Mantel: The Assassination of Margaret Thatcher - August 6th 1983; NSA Interceptng Gmail During Delivery; General IIR Filters; What's happening with Scala?; Interesting (But Largely Illegible) Typeface; Retiring Essentialism; Poorest in UK, Poorest in N Europe; I Want To Be A Redneck!; Reverse Racism; The Lost Art Of Nomography; IBM Data Center (Photo); Interesting Account Of Gamma Hack; The Most Interesting Audiophile In The World; How did the first world war actually end?; Ky - Restaurant Santiago; The Black Dork Lives!; The UN Requires Unaninmous Decisions; LPIR - Steganography in Practice; How I Am 6; Clear Explanation of Verizon / Level 3 / Netflix; Teenage Girls; Formalising NSA Attacks; Switching Brakes (Tektro Hydraulic); Naim NAP 100 (Power Amp); AKG 550 First Impressions; Facebook manipulates emotions (no really); Map Reduce "No Longer Used" At Google; Removing RAID metadata; New Bike (Good Bike Shop, Santiago Chile); Removing APE Tags in Linux; Compiling Python 3.0 With GCC 4.8; Maven is Amazing; Generating Docs from a GitHub Wiki; Modular Shelves; Bash Best Practices; Good Emergency Gasfiter (Santiago, Chile); Readings in Recent Architecture; Roger Casement; Integrated Information Theory (Or Not); Possibly undefined macro AC_ENABLE_SHARED; Update on Charges; Sunburst Visualisation; Spectral Embeddings (Distances -> Coordinates); Introduction to Causality; Filtering To Help Colour-Blindness; ASUS 1015E-DS02 Too; Ready Player One; Writing Clear, Fast Julia Code; List of LatAm Novels; Running (for women); Building a Jenkins Plugin and a Jar (for Command Line use); Headphone Test Recordings; Causal Consistency; The Quest for Randomness; Chat Wars; Real-life Financial Co Without ACID Database...; Flexible Muscle-Based Locomotion for Bipedal Creatures; SQL Performance Explained; The Little Manual of API Design; Multiple Word Sizes; CRC - Next Steps; FizzBuzz; Update on CRCs; Decent Links / Discussion Community; Automated Reasoning About LLVM Optimizations and Undefined Behavior; A Painless Guide To CRC Error Detection Algorithms; Tests in Julia; Dave Eggers: what's so funny about peace, love and Starship?; Cello - High Level C Programming; autoreconf needs tar; Will Self Goes To Heathrow; Top 5 BioInformatics Papers; Vasovagal Response; Good Food in Vina; Chilean Drug Criminals Use Subsitution Cipher; Adrenaline; Stiglitz on the Impact of Technology; Why Not; How I Am 5; Lenovo X240 OpenSuse 13.1; NSA and GCHQ - Psychological Trolls; Finite Fields in Julia (Defining Your Own Number Type); Julian Assange; Starting Qemu on OpenSuse

© 2006-2013 Andrew Cooke (site) / post authors (content).

Differentiating in Python

From: "andrew cooke" <andrew@...>

Date: Sun, 31 Dec 2006 20:21:52 -0300 (CLST)

I read this post -
http://wmfarr.blogspot.com/2006/10/automatic-differentiation-in-ocaml.html
- and realised I had written the same code in Python.  I hunted around and
found it on my (password protected) diary -
http://www.acooke.org/andrew/diary/2004/mar/4.html - so here it is in public.

It parses, differentiates, and then simplifies simple numerical terms.

Andrew


 output:

the original expression is a+3*b*a
it was parsed to (a+(3*(b*a)))
the differential wrt a is (1+(3*b))

the original expression is a/b
it was parsed to (a/b)
the differential wrt b is (0-(a*(1/(b*b))))

the original expression is a+b+c+d
it was parsed to (a+(b+(c+d)))
the differential wrt b is 1

the original expression is a+b*c+d
it was parsed to (a+((b*c)+d))
the differential wrt b is c

code:

# calculate first derivatives of an arithmetic expression
# basic code using just +,-,/,*,(), integers and variables (lowercase), #
but functions aren't any harder conceptually

# here's the interesting bit
# walk the ast to calculate the derivative wrt some variable

# you'd add functions in the normal way - for example sin(...) would # map
to diff(...) * sin(...) + ... * cos(....)

def diffwrt(n,var):
    def ifn(n): return 0
    def vfn(v):
        if v is var: return 1
        else: return 0
    def nfn(op, n1, n2, dn1, dn2):
        if op is '+' or op is '-':
            return (op, dn1, dn2)
        elif op is '*':
            return ('+',
                    ('*', dn1, n2),
                    ('*', n1, dn2))
        elif op is '/':
            return ('-',
                    ('/', dn1, n2),
                    ('*', n1,
                     ('/', dn2, ('*', n2, n2))))
    return folddown(ifn, vfn, nfn, n)

def tidy(n):
    def ifn(n): return n
    def vfn(v): return v
    def nfn(op, n1, n2, dn1, dn2):
        if isinstance(dn1, int):
            if dn1 is 0:
                if op is '+': return dn2
                elif op is '*': return 0
                elif op is '/': return 0
            elif dn1 is 1 and op is '*': return dn2
        if isinstance(dn2, int):
            if dn2 is 0:
                if op is '+': return dn1
                elif op is '*': return 0
                elif op is '/': raise "division by zero"
            elif dn2 is 1 and op is '*': return dn1
        if isinstance(dn1, int) and isinstance(dn2, int):
            if op is '+': return dn1 + dn2
            elif op is '-': return dn1 - dn2
            elif op is '*': return dn1 * dn2
            elif op is '/': return dn1 / dn2
        return (op, dn1, dn2)
    return folddown(ifn, vfn, nfn, n)

def folddown(ifn, vfn, nfn, n):
    if isinstance(n, int): return ifn(n)
    elif isinstance(n, str): return vfn(n)
    else:
        (op,n1,n2) = n
        (dn1, dn2) = (folddown(ifn,vfn,nfn,n1), folddown(ifn,vfn,nfn,n2))
return nfn(op, n1, n2, dn1, dn2)

# a "simple" recursive descent parser

# grammar:
# expr: term ((+|-) term)*
# term: fact ((*|/) fact)*
# fact: '(' expr ')' | var | num

# the ast is just (operator, node, node) tuples

# utilities

def dropspace(s):
    if s and s[0] is ' ': return dropspace(s[1:])
    else: return s

def empty(s): return dropspace(s) is ""

# so these are a bunch of 'recognisers' (eg cousineau + mauny)
# (you can think of them as tokenizers - they either return a match plus #
the remaning text or None)

def add(s): return mkonechar('+')(s)
def subtract(s): return mkonechar('-')(s)
def multiply(s): return mkonechar('*')(s)
def divide(s): return mkonechar('/')(s)
def openbracket(s): return mkonechar('(')(s)
def closebracket(s): return mkonechar(')')(s)

def variable(s): return mkmanychar(lambda c : c >= 'a' and c <= 'z')(s)
def number(s): return mkmanychar(lambda c : c >= '0' and c <= '9')(s)

def mkonechar(c):
    def localonechar(s):
        ss = dropspace(s)
        if ss and ss[0] is c: return (c,ss[1:])
        else: return None
    return lambda s: localonechar(s)

def mkmanychar(p):
    def accum(s,id=""):
        if s and p(s[0]): return accum(s[1:],id+s[0])
        elif id != "": return (id,s)
        else: return None
    return lambda s: accum(dropspace(s))

# this handles the "nxt ((p1|p2) nxt)*" structure in the grammar

def mkextend(p1,p2,nxt):
    def localextend(n1,s):
        if p1(s):
            (x,s) = p1(s)
            (n2,s) = nxt(s)
            return ((x,n1,n2),s)
        elif p2(s):
            (x,s) = p2(s)
            (n2,s) = nxt(s)
            return ((x,n1,n2),s)
        else: return (n1,s)
    return lambda n,s: localextend(n,s)

# and this is the parser itself

def expr(s):
    if term(s):
        (n,s) = term(s)
        return extendexpr(n,s)
    else: raise ("cannot parse " + s)

def extendexpr(n,s): return mkextend(add,subtract,expr)(n,s)

def term(s):
    if fact(s):
        (n,s) = fact(s)
        return extendterm(n,s)
    else: raise ("cannot parse " +s)

def extendterm(n,s): return mkextend(multiply,divide,term)(n,s)

def fact(s):
    if openbracket(s):
        (x,s) = openbracket(s)
        if expr(s):
            (n,s) = expr(s)
            if closebracket(s):
                (x,s) = closebracket(s)
                return (n,s)
            else: raise ("missing ): " + s)
        else: raise ("cannot parse: " + s)
    elif variable(s):
        (x,s) = variable(s)
        return (x,s)
    elif number(s):
        (x,s) = number(s)
        return (int(x),s)
    else: raise ("cannot parse: " + s)

# finally, a pretty printer

def asttostring(n):
    if isinstance(n, int): return str(n)
    elif isinstance(n, str): return n
    else: return astoptostring(n)

def astoptostring((op,n1,n2)):
    return "("+asttostring(n1)+op+asttostring(n2)+")"

# and test

def demo(text, var):
    (ast,x) = expr(text)
    ast2 = tidy(ast)
    print
    print "the original expression is", text
    print "it was parsed to", asttostring(ast2)
    diff = tidy(diffwrt(ast2, var))
    print "the differential wrt", var, "is", asttostring(diff)
    print

demo("a+3*b*a", "a")
demo("a/b", "b")
demo("a+b+c+d", "b")
demo("a+b*c+d", "b")

This Differentiates Strings

From: "andrew cooke" <andrew@...>

Date: Sun, 31 Dec 2006 21:53:52 -0300 (CLST)

I just realised that the code I linked to differentiates actual functions
(although it's not quite as cool as that sounds - rather than using
reflection to handle expressions directly (which is probably impossible in
OCaml, but would work in Lisp, I supposed)) it *appears* (I don't follow
understand the code) to define its own datatypes for the different
functions.

Since I don't fully understand things I'm not sure, but I think that means
it's equivalent to my code less the parsing.  If that's right I'm
surprised at how much code is needed (most of my code is parsing, although
I don't handle as wide a range of functions).

Suspect I'm missing something.  Maybe the idea is that once can import the
defined module rather than some standard moduls, without changing any
code, and get differentitation "for free".  Which there was more
explanation with the code...

Andrew

Differentiating Functions is Important

From: Will M Farr <farr@...>

Date: Thu, 4 Jan 2007 10:42:26 -0500

Andrew,

I think you have it correct in your comment above, but I'd like to  
emphasize that differentiating functions (as I do) is very  
important.  Using reflection (or lisp-style macros) to obtain a  
textual representation of a function (as you do) doesn't work well  
when you compose functions together, or use non-mathematical  
operators in a function.  For example, what's the derivative of

fun x ->
	if x > 0 then
		x
	else
		0 - x

at x = 3?  (You could do this textually if your text processor knew  
enough to avoid processing the if statement, but in order to make  
this work in general, you would have to process the whole language.)   
How about this:

fun x ->
	let y = 0 - x in
	if y < 0 then
		x
	else
		0 - y

It gets even worse if you define a function like

fun x ->
	let y = other_fun x in
	let z = another_fun y in
	x *. y *. z

To process this textually, you need a database which stores the text  
of other_fun and another_fun so they can be re-differentiated w.r.t x  
and then y (respectively).  Don't forget to apply the chain rule to  
the derivative of another_fun (since the argument is y)!  It's really  
a mess.

By the way, your last paragraph is entirely correct---I just take  
some existing code, add "open Deriv" to the top, place (C ...) in  
front of all numerical constants, and then I get derivatives for free.

A disadvantage of my method, which yours doesn't share, is that it  
would compute the derivative of this function:

fun x ->
	x +. 3.0 -. x

as

fun x ->
	1.0 +. 0.0 -. 1.0

I don't do any simplification, because I don't have the text of the  
expressions available at all.

Thanks for the interesting post---it's fun to see other people doing  
this kind of stuff!

Will

Re: Differentiating Functions is Important

From: "andrew cooke" <andrew@...>

Date: Thu, 4 Jan 2007 13:10:12 -0300 (CLST)

Hi,

Thanks for replying!  I just re-read my comment and it sounded more
negative than I intended because when I wrote the first half I still
hadn't worked out what I wrote in the second half (if you see what I
mean!) - since that final guess is correct this really is pretty sweet.

The only thing we may disagree on is that I think Lisp could combine the
best of both worlds, in that a Lisp macro effectively gives you both the
text and the "real function" (but I don't have enough experience to know
how the library you refer to would be handled/avoided - perhaps it
cannot).  Maybe MetaML would allow something similar to Lisp in ML, but
then you probably wouldn't be able to integrate the new functions into
existing code so easily (I think that's an excellent example of how cool
ML's module system is).  And of course, with Lisp you don't have a decent
type system :o)

Cheers,
Andrew

Caseless (test)

From: "andrew cooke" <andrew@...>

Date: Thu, 4 Jan 2007 14:11:00 -0300 (CLST)

Just found a bug in the blog system - emails can be forced to lower case,
so I have changed things to make IDs caseless.  This is to test whether it
works...

Andrew

Comment on this post