Fitting Data to Restriction of Known Function?

So I have a list of data points living in the $(x,y)$-plane. These data points come from a complicated curve in this plane whose exact form is not known! I am hoping to evaluate the Weirstrass-$\wp$ function $\wp(x+iy)$, on these data points and use some sort of fitting procedure to get the restriction of $\wp$, to this unknown curve, as a function of some length parameter along the curve. Numerics is not my strong suit, but there seems to be two ways potentially of doing this:

(I) First use a fit to find approximate closed form the the curve.

(This isn't ideal for me because in my case, the curve will have many turning points; certainly it won't be a function of $x$)

(II) To me, it seems that the better way would be to first evaluate the real and imaginary parts of $\wp$ on the known data points, and then apply a fitting procedure to this collection of data.

I've tried this second method using polynomial fitting, and needless to say, the results are terible given that Weirstrass-$\wp$ is a highly non-trivial beast. So my question is, given that I know exactly the function I'm trying to fit to, is there an effective way to carry out this fitting? I know exactly the function, I'm just restricting to a non-trivial subset of the domain given by known data points, and hoping to get an approximate function of length on the curve.

EDIT: Here is the code. I hesitated to include it because there's a step where the user has to manually "erase" data points using an interactive eraser. I want to erase all data points except the one "sawtooth" curve. Will be clear what I mean in the plot.

M = 1;
tau = (0.08) + (0.04)*I;
w1 = Pi/2;
w2 = Pi*(tau)/2;
inv = WeierstrassInvariants[{w1, w2}];
E2[t_] :=
1 - 24*Sum[(n*Exp[2*Pi*I*(t)*n])/(1 - Exp[2*Pi*I*(t)*n]), {n, 1,
300}];
z[u_] := (I*
M/2)*(WeierstrassZeta[u, inv] - ((1/3)*N[E2[tau], 50]*(u)));
WP[x_, y_] := WeierstrassP[w1*x + w2*y, inv];
L = -(1/3)*N[E2[tau], 50];
f[x_, y_] := Re[WP[x, y] - L];
g[x_, y_] := Im[WP[x, y] - L];
V1 = Quiet[
FindRoot[{f[x, y] == 0, g[x, y] == 0}, {x, 0.25}, {y, -1.5},
WorkingPrecision -> 50]];
V2 = Quiet[
FindRoot[{f[x, y] == 0, g[x, y] == 0}, {x, 1.75}, {y, 3.5},
WorkingPrecision -> 50]];
V3 = Quiet[
FindRoot[{f[x, y] == 0, g[x, y] == 0}, {x, 0.25}, {y, -3.5},
WorkingPrecision -> 50]];
V4 = Quiet[
FindRoot[{f[x, y] == 0, g[x, y] == 0}, {x, 1.75}, {y, 1.5},
WorkingPrecision -> 50]];
A1 = x /. V1;
B1 = y /. V1;
A2 = x /. V2;
B2 = y /. V2;
A3 = x /. V3;
B3 = y /. V3;
A4 = x /. V4;
B4 = y /. V4;
Z1 = Quiet[N[z[w1*A1 + w2*B1], 50]]
Z2 = Quiet[N[z[w1*A2 + w2*B2], 50]]
Z3 = Quiet[N[z[w1*A3 + w2*B3], 50]]
Z4 = Quiet[N[z[w1*A4 + w2*B4], 50]]
m = (Im[Z1] - Im[Z2])/((Re[Z1] - Re[Z2]));
Zed[x_, y_] := z[w1*x + w2*y];
ContourPlot[
Im[N[Zed[x, y]]] - (m*(Re[N[Zed[x, y] - (M/2)]]) + Im[M/2]) == 0, {x,
0, 2}, {y, -5, 7}, MaxRecursion -> 4,
RegionFunction ->
Function[{x, y, f}, N@Abs[Zed[x, y] - M/2] <= Abs[Z2 - M/2]]]
A = Union[%[[1, 1]]];
DynamicModule[{pt = {0, 0}, r = 3},
Column[{Slider[Dynamic[r], {0.1, 5}],
Button["Delete", A = Select[A, EuclideanDistance[#, pt] > r &]],
LocatorPane[Dynamic[pt],
Dynamic[ListPlot[A, Frame -> True, Axes -> False,
AspectRatio -> 1, Epilog -> Circle[pt, r], ImageSize -> 300,
PlotRange -> {{0, 2}, {-5, 5}}]], Appearance -> None]}]]
cc = ConstantArray[1, Length[A]];
Do[cc[[n]] = A[[n]][], {n, 1, Length[A]}];
xx = ConstantArray[1, Length[A]];
yy = ConstantArray[1, Length[A]];
Do[xx[[n]] = Re[N[Zed[A[[n]][], A[[n]][]]] - M/2], {n, 1,
Length[A]}];
Do[yy[[n]] = Im[N[Zed[A[[n]][], A[[n]][]]] - M/2], {n, 1,
Length[A]}];
dataxx = Transpose@{cc, xx};
datayy = Transpose@{cc, yy};
FitPolynomial[data_] := Fit[data, Table[x^n, {n, 0, 50}], x];
LLL = FitPolynomial[dataxx];
MMM = FitPolynomial[datayy];
FFF[x_] := Evaluate[LLL]
GGG[x_] := Evaluate[MMM]
zz1 = ConstantArray[1, Length[A]];
zz2 = ConstantArray[1, Length[A]];
Do[zz1[[n]] = Re[N[WP[A[[n]][], A[[n]][]]]], {n, 1,
Length[A]}];
Do[zz2[[n]] = Im[N[WP[A[[n]][], A[[n]][]]]], {n, 1,
Length[A]}];
datazz1 = Transpose@{cc, zz1};
datazz2 = Transpose@{cc, zz2};
ZZZ1 = FitPolynomial[datazz1];
ZZZZ1[x_] := Evaluate[ZZZ1];
ZZZ2 = FitPolynomial[datazz2];
ZZZZ2[x_] := Evaluate[ZZZ2];
Plot[ZZZZ1[x], {x, 0, 2}]
Plot[ZZZZ2[x], {x, 0, 2}]

It's datazz1 and datazz2which correspond to Weirstrass-$\wp$ (real and imaginary parts) evaluated on the data points. It is precisely these two collections that I'm trying to fit! You can see my terrible results in the final two plots.

• Can you provide sample data and the code where you transform it with the Weierstrass function? – C. E. Aug 6 '16 at 23:34
• @C.E. Done, sorry about that. I hesitated because it's kind of long, and there's a step where one needs to use an interactive "eraser" to get rid of transient data points. If you had thoughts, I'd appreciate it! But I understand it's kind of a bear. – Benighted Aug 6 '16 at 23:47
• It's usually best to avoid loops where possible, and run the code like: cc = Table[A[[n]][], {n, Length[A]}] Not just because it is quicker, but it's a lot clearer too, that way you also don't have to pre-create an array. Btw Something is giving me "Infinite expression $\frac{1}{0. +0.\mathbb{i}}$ encountered." – Feyre Aug 7 '16 at 8:36