Location along actual (beveled) curve instead of points

Question

The following code generates a tree root structure, by making each sub-root start at a point along the previous sub-root. They are all NURBS paths. The hitch, as you will see, is that the points are not inside the visible beveled curve, and therefore the sub-root's origin point gets off-set. What code should I write to make sure that the sub-roots always start inside the visible part of their parent root, instead of the points?

I've made a screencast to show what the problem is. I need a way to resolve it in the code itself, to make the whole model generative without assistance.

Update: here's the blender file

import bpy
import random
import mathutils

def pathPointLoc(cpath, points):        #uses the locations in the lsit (points) to shape the curve

    cpath.points.add(len(points)-1)

    for (index, point) in enumerate(points):
        cpath.points[index].co = point

    cpath.use_endpoint_u = True

    return

cu = bpy.data.curves.new("root", "CURVE")       #lines 16-26: generates the NURBS path curve for the tap (main) root
ob = bpy.data.objects.new("taproot", cu)        #these are not in a for loop because there only needs to be one tap root
polyline = cu.splines.new('NURBS')

scn = bpy.context.scene
scn.objects.link(ob)
scn.objects.active = ob

cu.dimensions = '3D'
cu.bevel_object = bpy.data.objects["circle1"]   #bezier circle, for bevel and consistent taper; circle1 is largest, circle3 is smallest
cu.taper_object = bpy.data.objects["circle1"]

x = random.randint(1,5)
y = random.randint(1,5)
z = 1                                           #keeping z axis constant, only need it to operate on x and y for now
w = 1                                           #honestly no idea what this element does but 4 values are necessary for the syntax

print("X and Y are "+str(x)+" and "+str(y))     #just to keep track of the values, not essential

if y > x:                                       #the curves behave strangely when y is less than x, hence this statement
    x = random.randint(y,y+2)
    print ("Revised values of above X and Y: "+str(x)+" and "+str(y))

tap_ur = (x,y,z,w)
tap_pts = [tap_ur]                              #ur for origin, because this is the base of each curve

for tapinc in range(1,9):                       #loop apends location tuples for each point on the curve, after the first one is defined

    y_noise = round(random.uniform(-1,1))*3     #since this curve has to move further on the x axis than the y, the y has a small - and + motion
    while y_noise == 0:
        y_noise = round(random.uniform(-1,1))*3 #0 avoided to keep the curve meandering, not constant

    x_noise = round(random.uniform(1,2))*3
    while x_noise == 0:
        x_noise = round(random.uniform(1,2))*3

    y += y_noise                                #meanders the curve
    x += x_noise                                #ensure that the next point on the curve is further along

    tap_nxt = list(tap_ur)                      #converting tuple to list in order to change x and y
    tap_nxt[0] = x
    tap_nxt[1] = y

    tap_pts.append(tap_nxt)

tap_pts.reverse()                               #because the curve takes locations from right left in the list; could have switched the path's direction instead but this is simpler

pathPointLoc(polyline, tap_pts)

for sub1 in range (1,len(tap_pts)):             #subroot, tier 1: looped version of the code from 16-61; note that the number of subroots depends on the length of the parent root

    cu = bpy.data.curves.new("root", "CURVE")
    ob = bpy.data.objects.new("subroot_1_"+str(sub1), cu)
    polyline = cu.splines.new('NURBS')

    scn = bpy.context.scene
    scn.objects.link(ob)
    scn.objects.active = ob

    cu.dimensions = '3D'
    cu.bevel_object = bpy.data.objects["circle2"]
    cu.taper_object = bpy.data.objects["circle2"]

    sub1_ur = random.choice(tap_pts[2:len(tap_pts)-1])
    sub1_pts = [sub1_ur]

    x = sub1_ur[0]                              #inheriting/resetting x and y values
    y = sub1_ur[1]

    for sub1inc in range(1,5):

        x_noise = round(random.uniform(-1,1))*2 #since this root tier is perpendicular to the parent one, x meanders it and y pushes it along
        while x_noise == 0:
            x_noise = round(random.uniform(-1,1))*2

        y_noise = round(random.uniform(1,2))*2
        while y_noise == 0:
            y_noise = round(random.uniform(1,2))*2

        x += x_noise

        if sub1%2 != 0:
            y += y_noise
        else:
            y -= y_noise

        sub1_nxt = list(sub1_ur)
        sub1_nxt[0] = x
        sub1_nxt[1] = y

        sub1_pts.append(sub1_nxt)

    x = sub1_ur[0]                              #resetting x and y to prevent carry-over
    y = sub1_ur[1]

    sub1_pts.reverse()

    pathPointLoc(polyline, sub1_pts)

    for sub2 in range (1,len(sub1_pts)):

        cu = bpy.data.curves.new("root", "CURVE")
        ob = bpy.data.objects.new("subroot_2_"+str((sub1-1)+sub2), cu)
        polyline = cu.splines.new('NURBS')

        scn = bpy.context.scene
        scn.objects.link(ob)
        scn.objects.active = ob

        cu.dimensions = '3D'
        cu.bevel_object = bpy.data.objects["circle3"]
        cu.taper_object = bpy.data.objects["circle3"]

        sub2_ur = random.choice(sub1_pts[2:len(sub1_pts)-1])
        sub2_pts = [sub2_ur]

        x = sub2_ur[0]
        y = sub2_ur[1]

        for sub2inc in range(1,5):

            y_noise = round(random.uniform(-1,1))
            while y_noise == 0:
                y_noise = round(random.uniform(-1,1))

            x_noise = round(random.uniform(1,2))
            while x_noise == 0:
                x_noise = round(random.uniform(1,2))        
            y += y_noise

            if sub2%2 != 0:
                x += x_noise
            else:
                x -= x_noise

            sub2_nxt = list(sub2_ur)
            sub2_nxt[0] = x
            sub2_nxt[1] = y

            sub2_pts.append(sub2_nxt)

        x = sub2_ur[0]
        y = sub2_ur[1]

        sub2_pts.reverse()

        pathPointLoc(polyline, sub2_pts)