Source

is-solver / partial_derivative / numerical_diff_fps_runner.py

Full commit
#!/usr/bin/env python

## This is a simple runner to take a Five-point stencil numerical derivativen
##
## FPS: http://en.wikipedia.org/wiki/Five-point_stencil

##################################################
## CONSTANTS 
H = 0.0001
##################################################

import os
import optparse
import uuid
from OrderedDict import * 
##import subprocess
import subprocess32 as subprocess
import numpy as np

##################################################
##################################################

### this is a custom tailored function
def getMomentForAllPoints(uuid, moment):
    PTS_LIST = ['point-A', 'point-B', 'point-C', 'point-D']
    IFEXT = '___results.csv'
    
    outdict = OrderedDict()

    for pt in PTS_LIST:
        pttmp = uuid + "_" + pt + IFEXT
        fo = open(pttmp, 'r')
        for line in fo:
            lineAry = line.split(",")
            if lineAry[0] == moment:
                outdict[pt] = float(lineAry[1])
        fo.close()

    return outdict


### make the derivative approximation
def makeDerivativeApproximationFPS(pointDict, H):
    pA = pointDict['point-A']
    pB = pointDict['point-B']
    pC = pointDict['point-C']
    pD = pointDict['point-D']
    fps = ((-1 * pA) + (8 * pB) - (8 * pD) + pC)/ (12*H)
    return(fps)


def mprint(fhlist, mstr):
    for fh in fhlist:
        fh.write(mstr)
    
strlist = lambda x: str(x)

def makeCSV(mat, rownames, colnames, fhlist):
    ## rownames and colnames are a list
    colstr = ",".join(colnames)
    mprint(fhlist, "moment," + colstr + "\n")
    Nrows = mat.shape[0]
    Ncols = mat.shape[1]
    for i in range(Nrows):
        ltmp = list(mat[i,:])
        ltmp = map(strlist, ltmp)
        rstr = ",".join(ltmp)
        mprint(fhlist, rownames[i] + "," + rstr + "\n")


##################################################
## parse some options
parser = optparse.OptionParser()
parser.add_option('-n', help='specify a comma-separated list of parameters', dest='paramNameList')
parser.add_option('-x', help='specify a comma-separated list of points corresponding to parameter list (same order)', dest='xvec')
parser.add_option('-m', help='specify a comma-separated list of output moments', dest='momentNameList')
parser.add_option('-e', help='specify an executable name', dest='EXECNAME')
parser.add_option('-o', help='specify an output file (CSV) [REQUIRED]', dest='ofile')
(opts, args) = parser.parse_args()
##################################################

paramNameList = opts.paramNameList
xvec = opts.xvec
momentList = opts.momentNameList
EXECNAME = opts.EXECNAME
ofile = opts.ofile

# paramNameList = "a,b"
# xvec = "1.5,3.0"
# momentList = "mean,var"
# EXECNAME = "testpoly"
CSV_POSTFIX = "results.csv"
UUID_PREFIX = "RESULTS"

##################################################


paramAry = paramNameList.split(",")
xAry = xvec.split(",")
momentAry = momentList.split(",")


##################################################
## generate the exec base string

#exec_base_str = "./" + EXECNAME + " " + "--moments='" + momentList + "' "
#exec_base_str = EXECNAME + " " + "--momentlist='" + momentList + "' "
exec_base_str = EXECNAME + " "
guid = str(uuid.uuid4())
execStrDict = OrderedDict()
ugkey_list = []

OUTFOLDER = os.path.join(UUID_PREFIX, guid)
if not os.path.exists(OUTFOLDER):
    os.makedirs(OUTFOLDER)
    print "CREATED FOLDER: " + OUTFOLDER

### perturb each parameter in the following way and see how it affects
### each moment
for pp in range(len(paramAry)):
    paramNameTmp = paramAry[pp]
    param_xTmp = float(xAry[pp])
    point_A = param_xTmp + 2*H
    point_B = param_xTmp + H
    point_C = param_xTmp - 2*H
    point_D = param_xTmp - H
    pointsDict = OrderedDict()
    pointsDict['A'] = point_A
    pointsDict['B'] = point_B
    pointsDict['C'] = point_C
    pointsDict['D'] = point_D
    uuid_i = str(uuid.uuid4())
    ugkey =  UUID_PREFIX + "/" + guid + "/" + uuid_i
    ugkey_list.append(ugkey)
    for point in pointsDict.keys():
        xAry_tmp = xAry[:]
        xAry_tmp[pp] = pointsDict[point]
        
        exec_detailed = ""
        for ppi in range(len(paramAry)):
            exec_detailed += " --" + str(paramAry[ppi]) + "='" + str(xAry_tmp[ppi]) + "'"

        ugkey_pt = ugkey + "_point-" + point
        exec_detailed += " --uuid_i='" + ugkey_pt + "'"
    
        exec_full = exec_base_str + exec_detailed
        execStrDict[ugkey_pt] = exec_full



## at the end of this we have the following state
## -- execStrDict.values() are all of the exec strs we need to evaluate
## to get a result
## -- ugkey_list are unique keys, one per input parameter
## -- execStrDict.keys() are going to be where files are located vor
## the various points, currently points A,B,C, and D for each 

## now, for every line in execStrDict.values(), subprocess that
## then, get the results back from the csv file
count = 0        
for run_this_string in execStrDict.values():
    print(str(count) + ":  About to run: " + run_this_string)
    subprocess.check_output(run_this_string, shell=True)

exit


##
## CHECKPOINT
## check to make sure the results are there, then calculate the
## numerical derivatives. Record to a log file and then quit
##

for esd in execStrDict.keys():
    rfile = esd + "___results.csv"
    if os.path.exists(rfile):
        print "Found: " + rfile
    else:
        print "MISSING: " + rfile


## make the derivative matrix
DERR = np.empty((len(momentAry), len(paramAry)))

##################################################
##
## USE POINT ESTIMATES TO FILL IN THE MATRIX
##
## pp is the column count
## mm is the row count
for pp in range(len(paramAry)):
    paramTmp = paramAry[pp]
    print("Calculating DERR for paramter: " + paramTmp)
    ugkey_tmp = ugkey_list[pp]

    ## for ugkey_tmp, open A,B,C,D
    ## udict is an OrderedDict such that you can do
    ## udict[moment][point_A]


    udict = OrderedDict()

    ## populate udict such that I can do udict[moment][point-{A,B,C,D}
    for mm in range(len(momentAry)):
        moment_tmp = momentAry[mm]
        udict[moment_tmp] = getMomentForAllPoints(ugkey_tmp, moment_tmp)

    ## calculate the moments
    for mm in range(len(momentAry)):
        moment_tmp = momentAry[mm]
        fpsApprox = makeDerivativeApproximationFPS(udict[moment_tmp], H)
        print("Partial derivative of " + moment_tmp + " with respect to " + paramTmp)
        DERR[ mm, pp ] = fpsApprox

print(DERR)
otherOutputFile = os.path.join(UUID_PREFIX, guid, 'derivative.csv')

of = open(ofile, 'w')
oof = open(otherOutputFile, 'w')
fhlist = [of, oof]

## def makeCSV(mat, rownames, colnames, fhlist):                      
makeCSV(DERR, momentAry, paramAry, fhlist)

closeall =  lambda fh: fh.close()
map(closeall, fhlist)