Offset points in space with python
Not perfect (but imperfections are like the one we have bevelling or insetting, etc. too far), and that can also be the case with a curve and bevel object.

The principle:


Start with the first two vertices of the curve, a segment S0 = (V0, V1)
Get the location of the starting position of the // curve P0
Calculate the displacement D0 = (V0, P0) 
Calculate the normal N0 = D0 cross (V0, V1)


Now iterate for 1 to n:


Get the new normal by D(i-1) cross S(i)
Get the new D(i) = S(i) cross N(i) 
The location is Pi = Vi + Di


Result with imperfection (the active curve with 4 empties as starting points):



Main part of the code is the following (complete code is in the blend file below, so this snippet won't work alone).

To use it from the blend file: the curve is the active object, starting points are selected object.

def create_parallel(curve_object, from_object):
    #Get the poly spline
    poly = curve_object.data.splines[0]

    #Create the target poly spline
    target_object = create_curve(curve_object, bpy.context.scene.collection)
    target_poly = target_object.data.splines[0]

    #Get the starting location in curve coordinates
    starting_location = curve_object.matrix_world.inverted() @ from_object.location

    #Delta from the first point
    delta = starting_location - poly.points[0].co.xyz
    #Its length
    distance = delta.length
    #Keep it normalized
    delta.normalize()
    #Orthogonal direction
    ortho = delta.cross( poly.points[1].co.xyz - poly.points[0].co.xyz ).normalized()

    #First vertex
    target_poly.points[0].co = (*starting_location, poly.points[0].co.w)

    for i in range(len(poly.points)-1):
        prev = poly.points[i].co
        co = poly.points[i+1].co

        axis = (co.xyz - prev.xyz)

        #Orthogonal direction is recalculated to be orthogonal to this delta and the next curve segment
        ortho = delta.cross(axis).normalized()

        #New delta is calculated as being orthogonal to the curve orientation and the current orthogonal direction
        delta = (axis).cross( ortho ).normalized()

        #Calculate the location of the current vertex        
        location = co.xyz + delta * distance
        #Assign it to the curve
        target_poly.points[i+1].co = (*location, co.w)




Improving it, should need to relax (or remove) vertices comparing the curvatures of the curves... another story...