问:
How to plot geodesic curves on a surface embedded in 3D?
I have in mind this video, or this simulation, and I would like to reproduce the geodesic lines on some sort of surface in 3D, given by a function f(x,y), from some starting point.

The midpoint method seems computationally and code intense, and so I'd like to ask if there is a way to generate an approximate geodesic curve based on the normal vector to the surface at different points. Each point has a tangent vector space associated with it, and therefore, it seems like knowing the normal vector does not determine a specific direction to move forward the curve.

I have tried working with Geogebra, but I realize that it may be necessary to shift to other software platforms, such as Python (or Poser?), Matlab, or others.

Is this idea possible, and can I get some ideas as to how to implement it?

In case it provides some ideas as to how to answer the question, there previously was an answer (now unfortunatley erased) suggesting the midpoint method for a terrain with the functional form z = F(x,y), starting with the straight line between the endpoints, splitting in short segments [I presume the straight line on the XY plane (?)], and lifting [I presume the nodes between segments on the XY plane (?)] on the surface. Next it suggested finding "a midpoint" [I guess a midpoint of the segments joining each consecutive pairs of projected points on the surface (?)], and projecting "it" [I guess each one of these midpoints close, but not quite on the surface(?)] orthogonally on the surface (in the direction of the normal), using the equation Z + t = F(X + t Fx, Y + t Fy) [I guess this is a dot product meant to be zero...

(?)], where (X,Y,Z) are the coordinates of the midpoint, Fx, Fy the partial derivatives of F, and t the unknown [that is my main issue understanding this... What am I supposed to do with this t once I find it? Add it to each coordinate of (X,Y,Z) as in (X+t, Y+t, Z+t)? And then?]. This is a non-linear equation in t, solved via Newton's iterations.

As an update / bookmark, Alvise Vianello has kindly posted a Python computer simulation of geodesic lines inspired on this page on GitHub. Thank you very much!

答:
Using the midpoint search method:

applied to the function f(x,y) = x^3 + y^2, I am projecting the points of the line segment on the XY plane y = x from x = -1 to x = 1.

To get an idea, with one iteration and only 4 points on the line on the XY plane, the black spheres are these 4 original points of the line projected on the surface, while the red dots are the midpoints in a single iteration, and the yellow dots the result of the projection of the red dots along the normal to the surface:

Using Matlab fmincon() and after 5 iterations we can get a geodesic from point A to point B:

Here is the code:

```<code>% Creating the surface
x = linspace(-1,1);
y = linspace(-1,1);
[x,y] = meshgrid(x,y);
z = x.^3 + y.^2;
S = [x;y;z];
h = surf(x,y,z)
set(h,'edgecolor','none')
colormap summer

% Number of points
n = 1000;

% Line to project on the surface with n values to get a feel for it...
t = linspace(-1,1,n);
height = t.^3 + t.^2;
P = [t;t;height];

% Plotting the projection of the line on the surface:
hold on
%plot3(P(1,:),P(2,:),P(3,:),'o')

for j=1:5
% First midpoint iteration updates P...
P = [P(:,1), (P(:,1:end-1) + P(:,2:end))/2, P(:,end)];
%plot3(P(1,:), P(2,:), P(3,:), '.', 'MarkerSize', 20)

A = zeros(3,size(P,2));
for i = 1:size(P,2)
% Starting point will be the vertical projection of the mid-points:
    A(:,i) = [P(1,i), P(2,i), P(1,i)^3 + P(2,i)^2];
end

% Linear constraints:
nonlincon = @nlcon;

% Placing fmincon in a loop for all the points

for i = 1:(size(A,2))
    % Objective function:
    objective = @(x)(P(1,i) - x(1))^2 + (P(2,i) - x(2))^2 + (P(3,i)-x(3))^2;
    A(:,i) = fmincon(objective, A(:,i), [], [], [], [], [], [], nonlincon);
end

P = A;
end

plot3(P(1,:), P(2,:), P(3,:), '.', 'MarkerSize', 5,'Color','y')
</code>```

In a separate file with name nlcon.m:

```<code>function[c,ceq] = nlcon(x)
   c   = [];
   ceq = x(3) - x(1)^3 - x(2)^2;
</code>```

Same for a geodesic on a really cool surface with a straight, non-diagonal line on XY:

```<code>% Creating the surface
x = linspace(-1,1);
y = linspace(-1,1);
[x,y] = meshgrid(x,y);
z = sin(3*(x.^2+y.^2))/10;
S = [x;y;z];
h = surf(x,y,z)
set(h,'edgecolor','none')
colormap summer

% Number of points
n = 1000;

% Line to project on the surface with n values to get a feel for it...
t = linspace(-1,1,n);
height = sin(3*((.5*ones(1,n)).^2+ t.^2))/10;
P = [(.5*ones(1,n));t;height];

% Plotting the line on the surface:
hold on
%plot3(P(1,:),P(2,:),P(3,:),'o')

for j=1:2
% First midpoint iteration updates P...
P = [P(:,1), (P(:,1:end-1) + P(:,2:end))/2, P(:,end)];
%plot3(P(1,:), P(2,:), P(3,:), '.', 'MarkerSize', 20)

A = zeros(3,size(P,2));
for i = 1:size(P,2)
% Starting point will be the vertical projection of the first mid-point:
    A(:,i) = [P(1,i), P(2,i), sin(3*(P(1,i)^2+ P(2,i)^2))/10];
end

% Linear constraints:
nonlincon = @nonlincon;

% Placing fmincon in a loop for all the points

for i = 1:(size(A,2))
    % Objective function:
    objective = @(x)(P(1,i) - x(1))^2 + (P(2,i) - x(2))^2 + (P(3,i)-x(3))^2;
    A(:,i) = fmincon(objective, A(:,i), [], [], [], [], [], [], nonlincon);
end

P = A;
end

plot3(P(1,:), P(2,:), P(3,:), '.', 'MarkerSize',5,'Color','r')
</code>```

with the nonlinear constraint in nonlincon.m:

```<code>function[c,ceq] = nlcon(x)
   c   = [];
   ceq = x(3) - sin(3*(x(1)^2+ x(2)^2))/10;
</code>```

One nagging concern is the possibility of overfitting to the curve with this method, and this latter plot is an example of it. So I adjusted the code to just select one beginning and one ending point, and allowing the iterative process to find the rest of the curve, which for 100 iterations seemed to be heading in the right direction:

The above examples seem to follow a linear projection on the XY plane, but fortunately this is not a fixed pattern, which would cast further doubt on the method. See for instance the hyperbolic paraboloid x^2 - y^2:

Notice that there are algorithms to advance or push geodesic lines along a surface f(x,y) with small increments determined by the starting points and the normal vector to the surface, as in here. Thanks to the work of Alvise Vianello looking into the JS in that simulation and his sharing in GitHub, I was able to turn that algorithm into Matlab code, generating this plot for the first example, f(x,y) = x^3 + y^2:

Here is the Matlab code:

```<code>x = linspace(-1,1);
y = linspace(-1,1);
[x,y] = meshgrid(x,y);
z = x.^3 + y.^2;
S = [x;y;z];
h = surf(x,y,z)
set(h,'edgecolor','none')
colormap('gray');
hold on

f = @(x,y) x.^3 + y.^2; % The actual surface

dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y

N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.

C = {'k','b','r','g','y','m','c',[.8 .2 .6],[.2,.8,.1],[0.3010 0.7450 0.9330],[0.9290 0.6940 0.1250],[0.8500 0.3250 0.0980]}; % Color scheme

for s = 1:11     % No. of lines to be plotted.
start = -5:5;    % Distributing the starting points of the lines.
y0 = start(s)/5; % Fitting the starting pts between -1 and 1 along y axis.
x0 = 1;          % Along x axis always starts at 1.
dx0 = 0;         % Initial differential increment along x
dy0 = 0.05;      % Initial differential increment along y
step_size = 0.000008; % Will determine the progression rate from pt to pt.
eta =  step_size / sqrt(dx0^2 + dy0^2); % Normalization.
eps = 0.0001;          % Epsilon
max_num_iter = 100000; % Number of dots in each line.

x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values

for i = 2:(max_num_iter - 1)  % Creating the geodesic:
            xt = x(i);        % Values at point t of x, y and the function:
            yt = y(i);
            ft = f(xt,yt);

            xtm1 = x(i - 1);  % Values at t minus 1 (prior point) for x,y,f
            ytm1 = y(i - 1);
            ftm1 = f(xtm1,ytm1);

            xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
            ysymp = yt + (yt - ytm1);
            fsymp = ft + (ft - ftm1);

            df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
            n = N(xt,yt);                % Normal vector at point t
            gamma = df * n(3);           % Scalar x change f x z value of N

            xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
            ytp1 = ysymp - gamma * n(2);

            x(i + 1) = xtp1;
            y(i + 1) = ytp1;
end

P = [x; y; f(x,y)]; % Compiling results into a matrix.

indices = find(abs(P(1,:)) < 1); % Avoiding lines overshooting surface.
P = P(:,indices);
indices = find(abs(P(2,:)) < 1);
P = P(:,indices);

    units = 15; % Deternines speed (smaller, faster)
    packet = floor(size(P,2)/units);
    P = P(:,1: packet * units);

  for k = 1:packet:(packet * units)
        hold on
        plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))),...
            '.', 'MarkerSize', 3.5,'color',C{s})
        drawnow
  end

end
</code>```

And here is an earlier example from above, but now calculated differently, and with lines starting side by side, following geodesics (no point-to-point trajectory):

```<code>    x = linspace(-1,1);
    y = linspace(-1,1);
    [x,y] = meshgrid(x,y);
    z = sin(3*(x.^2+y.^2))/10;
    S = [x;y;z];
    h = surf(x,y,z)
    set(h,'edgecolor','none')
    colormap('gray');
    hold on

    f = @(x,y) sin(3*(x.^2+y.^2))/10; % The actual surface

    dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
    dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y

    N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.

    C = {'k','r','g','y','m','c',[.8 .2 .6],[.2,.8,.1],[0.3010 0.7450 0.9330],[0.7890 0.5040 0.1250],[0.9290 0.6940 0.1250],[0.8500 0.3250 0.0980]}; % Color scheme

    for s = 1:11     % No. of lines to be plotted.
    start = -5:5;    % Distributing the starting points of the lines.
    x0 = -start(s)/5; % Fitting the starting pts between -1 and 1 along y axis.
    y0 = -1;          % Along x axis always starts at 1.
    dx0 = 0;         % Initial differential increment along x
    dy0 = 0.05;      % Initial differential increment along y
    step_size = 0.00005; % Will determine the progression rate from pt to pt.
    eta =  step_size / sqrt(dx0^2 + dy0^2); % Normalization.
    eps = 0.0001;          % Epsilon
    max_num_iter = 100000; % Number of dots in each line.

    x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
    y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values

    for i = 2:(max_num_iter - 1)  % Creating the geodesic:
                xt = x(i);        % Values at point t of x, y and the function:
                yt = y(i);
                ft = f(xt,yt);

                xtm1 = x(i - 1);  % Values at t minus 1 (prior point) for x,y,f
                ytm1 = y(i - 1);
                ftm1 = f(xtm1,ytm1);

                xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
                ysymp = yt + (yt - ytm1);
                fsymp = ft + (ft - ftm1);

                df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
                n = N(xt,yt);                % Normal vector at point t
                gamma = df * n(3);           % Scalar x change f x z value of N

                xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
                ytp1 = ysymp - gamma * n(2);

                x(i + 1) = xtp1;
                y(i + 1) = ytp1;
    end

    P = [x; y; f(x,y)]; % Compiling results into a matrix.

    indices = find(abs(P(1,:)) < 1); % Avoiding lines overshooting surface.
    P = P(:,indices);
    indices = find(abs(P(2,:)) < 1);
    P = P(:,indices);
    units = 35; % Deternines speed (smaller, faster)
    packet = floor(size(P,2)/units);
    P = P(:,1: packet * units);

  for k = 1:packet:(packet * units)
        hold on

        plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))), '.', 'MarkerSize', 5,'color',C{s})
        drawnow
  end

    end
</code>```

Some more examples:

```<code>    x = linspace(-1,1);
    y = linspace(-1,1);
    [x,y] = meshgrid(x,y);
    z = x.^2 - y.^2;
    S = [x;y;z];
    h = surf(x,y,z)
    set(h,'edgecolor','none')
    colormap('gray');

    f = @(x,y) x.^2 - y.^2; % The actual surface

    dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
    dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y

    N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.

    C = {'b','w','r','g','y','m','c',[0.75, 0.75, 0],[0.9290, 0.6940, 0.1250],[0.3010 0.7450 0.9330],[0.1290 0.6940 0.1250],[0.8500 0.3250 0.0980]}; % Color scheme

    for s = 1:11     % No. of lines to be plotted.
    start = -5:5;    % Distributing the starting points of the lines.
    x0 = -start(s)/5; % Fitting the starting pts between -1 and 1 along y axis.
    y0 = -1;          % Along x axis always starts at 1.
    dx0 = 0;         % Initial differential increment along x
    dy0 = 0.05;      % Initial differential increment along y
    step_size = 0.00005; % Will determine the progression rate from pt to pt.
    eta =  step_size / sqrt(dx0^2 + dy0^2); % Normalization.
    eps = 0.0001;          % Epsilon
    max_num_iter = 100000; % Number of dots in each line.

    x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
    y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values

    for i = 2:(max_num_iter - 1)  % Creating the geodesic:
                xt = x(i);        % Values at point t of x, y and the function:
                yt = y(i);
                ft = f(xt,yt);

                xtm1 = x(i - 1);  % Values at t minus 1 (prior point) for x,y,f
                ytm1 = y(i - 1);
                ftm1 = f(xtm1,ytm1);

                xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
                ysymp = yt + (yt - ytm1);
                fsymp = ft + (ft - ftm1);

                df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
                n = N(xt,yt);                % Normal vector at point t
                gamma = df * n(3);           % Scalar x change f x z value of N

                xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
                ytp1 = ysymp - gamma * n(2);

                x(i + 1) = xtp1;
                y(i + 1) = ytp1;
    end

    P = [x; y; f(x,y)]; % Compiling results into a matrix.

    indices = find(abs(P(1,:)) < 1); % Avoiding lines overshooting surface.
    P = P(:,indices);
    indices = find(abs(P(2,:)) < 1);
    P = P(:,indices);
    units = 45; % Deternines speed (smaller, faster)
    packet = floor(size(P,2)/units);
    P = P(:,1: packet * units);

  for k = 1:packet:(packet * units)
        hold on
        plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))), '.', 'MarkerSize', 5,'color',C{s})
        drawnow
  end

  end
</code>```

Or this one:

```<code>    x = linspace(-1,1);
    y = linspace(-1,1);
    [x,y] = meshgrid(x,y);
    z = .07 * (.1 + x.^2 + y.^2).^(-1);
    S = [x;y;z];
    h = surf(x,y,z)
    zlim([0 8])
    set(h,'edgecolor','none')
    colormap('gray');
    axis off
    hold on

    f = @(x,y) .07 * (.1 + x.^2 + y.^2).^(-1);    % The actual surface

    dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
    dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y

    N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.

     C = {'w',[0.8500, 0.3250, 0.0980],[0.9290, 0.6940, 0.1250],'g','y','m','c',[0.75, 0.75, 0],'r',...
         [0.56,0,0.85],'m'}; % Color scheme

    for s = 1:10     % No. of lines to be plotted.
    start = -9:2:9;
    x0 = -start(s)/10;
    y0 = -1;          % Along x axis always starts at 1.
    dx0 = 0;         % Initial differential increment along x
    dy0 = 0.05;      % Initial differential increment along y
    step_size = 0.00005; % Will determine the progression rate from pt to pt.
    eta =  step_size / sqrt(dx0^2 + dy0^2); % Normalization.
    eps = 0.0001;          % EpsilonA
    max_num_iter = 500000; % Number of dots in each line.

    x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
    y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values

    for i = 2:(max_num_iter - 1)  % Creating the geodesic:
                xt = x(i);        % Values at point t of x, y and the function:
                yt = y(i);
                ft = f(xt,yt);

                xtm1 = x(i - 1);  % Values at t minus 1 (prior point) for x,y,f
                ytm1 = y(i - 1);
                ftm1 = f(xtm1,ytm1);

                xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
                ysymp = yt + (yt - ytm1);
                fsymp = ft + (ft - ftm1);

                df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
                n = N(xt,yt);                % Normal vector at point t
                gamma = df * n(3);           % Scalar x change f x z value of N

                xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
                ytp1 = ysymp - gamma * n(2);

                x(i + 1) = xtp1;
                y(i + 1) = ytp1;
    end

     P = [x; y; f(x,y)]; % Compiling results into a matrix.

    indices = find(abs(P(1,:)) < 1.5); % Avoiding lines overshooting surface.
    P = P(:,indices);
    indices = find(abs(P(2,:)) < 1);
    P = P(:,indices);

    units = 15; % Deternines speed (smaller, faster)
    packet = floor(size(P,2)/units);
    P = P(:,1: packet * units);

  for k = 1:packet:(packet * units)
        hold on
        plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))),...
            '.', 'MarkerSize', 3.5,'color',C{s})
        drawnow
  end

    end
</code>```

Or a sinc function:

```<code>    x = linspace(-10, 10);
    y = linspace(-10, 10);
    [x,y] = meshgrid(x,y);
    z = sin(1.3*sqrt (x.^ 2 + y.^ 2) + eps)./ (sqrt (x.^ 2 + y.^ 2) + eps);
    S = [x;y;z];
    h = surf(x,y,z)
    set(h,'edgecolor','none')
    colormap('gray');
    axis off
    hold on

    f = @(x,y) sin(1.3*sqrt (x.^ 2 + y.^ 2) + eps)./ (sqrt (x.^ 2 + y.^ 2) + eps);   % The actual surface

    dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
    dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y

    N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.

    C = {'w',[0.8500, 0.3250, 0.0980],[0.9290, 0.6940, 0.1250],'g','y','r','c','m','w',...
         [0.56,0,0.85],[0.8500, 0.7250, 0.0980],[0.2290, 0.1940, 0.6250],'w',...
         [0.890, 0.1940, 0.4250],'y',[0.2290, 0.9940, 0.3250],'w',[0.1500, 0.7250, 0.0980],...
         [0.8500, 0.3250, 0.0980],'m','w'}; % Color scheme

    for s = 1:12     % No. of lines to be plotted.

    x0 = 10;
    y0 = 10;          % Along x axis always starts at 1.
    dx0 = -0.001*(cos(pi /2 *s/11));         % Initial differential increment along x
    dy0 = -0.001*(sin(pi /2 *s/11));         % Initial differential increment along y
    step_size = 0.0005; % Will determine the progression rate from pt to pt.
    % Making it smaller increases the length of the curve.
    eta =  step_size / sqrt(dx0^2 + dy0^2); % Normalization.
    eps = 0.0001;          % EpsilonA
    max_num_iter = 500000; % Number of dots in each line.

    x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
    y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values

    for i = 2:(max_num_iter - 1)  % Creating the geodesic:
                xt = x(i);        % Values at point t of x, y and the function:
                yt = y(i);
                ft = f(xt,yt);

                xtm1 = x(i - 1);  % Values at t minus 1 (prior point) for x,y,f
                ytm1 = y(i - 1);
                ftm1 = f(xtm1,ytm1);

                xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
                ysymp = yt + (yt - ytm1);
                fsymp = ft + (ft - ftm1);

                df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
                n = N(xt,yt);                % Normal vector at point t
                gamma = df * n(3);           % Scalar x change f x z value of N

                xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
                ytp1 = ysymp - gamma * n(2);

                x(i + 1) = xtp1;
                y(i + 1) = ytp1;
    end

     P = [x; y; f(x,y)]; % Compiling results into a matrix.

    indices = find(abs(P(1,:)) < 10); % Avoiding lines overshooting surface.
    P = P(:,indices);
    indices = find(abs(P(2,:)) < 10);
    P = P(:,indices);

    units = 15; % Deternines speed (smaller, faster)
    packet = floor(size(P,2)/units);
    P = P(:,1: packet * units);

  for k = 1:packet:(packet * units)
        hold on
        plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))),...
            '.', 'MarkerSize', 3.5,'color',C{s})
        drawnow
  end

    end
</code>```

And one very last one:

```<code>    x = linspace(-1.5,1.5);
    y = linspace(-1,1);
    [x,y] = meshgrid(x,y);
    z = 0.5 *y.*sin(5 * x) - 0.5 * x.*cos(5 * y)+1.5;
    S = [x;y;z];
    h = surf(x,y,z)
    zlim([0 8])
    set(h,'edgecolor','none')
    colormap('gray');
    axis off
    hold on

    f = @(x,y) 0.5 *y.* sin(5 * x) - 0.5 * x.*cos(5 * y)+1.5;     % The actual surface

    dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
    dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y

    N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.

     C = {'w',[0.8500, 0.3250, 0.0980],[0.9290, 0.6940, 0.1250],'g','y','k','c',[0.75, 0.75, 0],'r',...
         [0.56,0,0.85],'m'}; % Color scheme

    for s = 1:11     % No. of lines to be plotted.
    start = [0, 0.7835,  -0.7835, 0.5877, -0.5877, 0.3918, -0.3918, 0.1959, -0.1959, 0.9794, -0.9794];
    x0 = start(s);
    y0 = -1;          % Along x axis always starts at 1.
    dx0 = 0;         % Initial differential increment along x
    dy0 = 0.05;      % Initial differential increment along y
    step_size = 0.00005; % Will determine the progression rate from pt to pt.
    % Making it smaller increases the length of the curve.
    eta =  step_size / sqrt(dx0^2 + dy0^2); % Normalization.
    eps = 0.0001;          % EpsilonA
    max_num_iter = 500000; % Number of dots in each line.

    x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
    y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values

    for i = 2:(max_num_iter - 1)  % Creating the geodesic:
                xt = x(i);        % Values at point t of x, y and the function:
                yt = y(i);
                ft = f(xt,yt);

                xtm1 = x(i - 1);  % Values at t minus 1 (prior point) for x,y,f
                ytm1 = y(i - 1);
                ftm1 = f(xtm1,ytm1);

                xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
                ysymp = yt + (yt - ytm1);
                fsymp = ft + (ft - ftm1);

                df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
                n = N(xt,yt);                % Normal vector at point t
                gamma = df * n(3);           % Scalar x change f x z value of N

                xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
                ytp1 = ysymp - gamma * n(2);

                x(i + 1) = xtp1;
                y(i + 1) = ytp1;
    end

     P = [x; y; f(x,y)]; % Compiling results into a matrix.

    indices = find(abs(P(1,:)) < 1.5); % Avoiding lines overshooting surface.
    P = P(:,indices);
    indices = find(abs(P(2,:)) < 1);
    P = P(:,indices);

    units = 15; % Deternines speed (smaller, faster)
    packet = floor(size(P,2)/units);
    P = P(:,1: packet * units);

  for k = 1:packet:(packet * units)
        hold on
        plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))),...
            '.', 'MarkerSize', 3.5,'color',C{s})
        drawnow
  end

    end
</code>```

On a torus this becomes more complicated, and so far I haven't been able to adapt this code to it. As a reference, the code provided by Paul Chesler in here does provide with an option. It involves saving a file tor.m as

```<code>function[c,ceq] = nlcon(x)
   c   = [];
   ceq = x(3) - x(1)^3 - x(2)^2;
</code>```0

and in another file torus.m plot the torus and call the function above to solve a system of second order differential equations. Here it is for different geodesics:

from below:

with torus.m:

```<code>function[c,ceq] = nlcon(x)
   c   = [];
   ceq = x(3) - x(1)^3 - x(2)^2;
</code>```1

or

```<code>function[c,ceq] = nlcon(x)
   c   = [];
   ceq = x(3) - x(1)^3 - x(2)^2;
</code>```2

or

from below:

tspan = linspace(pi,100*pi,1000000); defines the limits of integration along 1e^6 pts, and [3/4 * pi,.1,-pi/2,.2], the parametric values for a starting point (on the last example).

On a sphere:

code here.

Body:
Using the midpoint search method:

applied to the function f(x,y) = x^3 + y^2, I am projecting the points of the line segment on the XY plane y = x from x = -1 to x = 1.

To get an idea, with one iteration and only 4 points on the line on the XY plane, the black spheres are these 4 original points of the line projected on the surface, while the red dots are the midpoints in a single iteration, and the yellow dots the result of the projection of the red dots along the normal to the surface:

Using Matlab fmincon() and after 5 iterations we can get a geodesic from point A to point B:

Here is the code:

% Creating the surface
x = linspace(-1,1);
y = linspace(-1,1);
[x,y] = meshgrid(x,y);
z = x.^3 + y.^2;
S = [x;y;z];
h = surf(x,y,z)
set(h,'edgecolor','none')
colormap summer

% Number of points
n = 1000;

% Line to project on the surface with n values to get a feel for it...
t = linspace(-1,1,n);
height = t.^3 + t.^2;
P = [t;t;height];

% Plotting the projection of the line on the surface:
hold on
%plot3(P(1,:),P(2,:),P(3,:),'o')

for j=1:5
% First midpoint iteration updates P...
P = [P(:,1), (P(:,1:end-1) + P(:,2:end))/2, P(:,end)];
%plot3(P(1,:), P(2,:), P(3,:), '.', 'MarkerSize', 20)

A = zeros(3,size(P,2));
for i = 1:size(P,2)
% Starting point will be the vertical projection of the mid-points:
    A(:,i) = [P(1,i), P(2,i), P(1,i)^3 + P(2,i)^2];
end

% Linear constraints:
nonlincon = @nlcon;

% Placing fmincon in a loop for all the points

for i = 1:(size(A,2))
    % Objective function:
    objective = @(x)(P(1,i) - x(1))^2 + (P(2,i) - x(2))^2 + (P(3,i)-x(3))^2;
    A(:,i) = fmincon(objective, A(:,i), [], [], [], [], [], [], nonlincon);
end

P = A;
end

plot3(P(1,:), P(2,:), P(3,:), '.', 'MarkerSize', 5,'Color','y')

In a separate file with name nlcon.m:

function[c,ceq] = nlcon(x)
   c   = [];
   ceq = x(3) - x(1)^3 - x(2)^2;

Same for a geodesic on a really cool surface with a straight, non-diagonal line on XY:

% Creating the surface
x = linspace(-1,1);
y = linspace(-1,1);
[x,y] = meshgrid(x,y);
z = sin(3*(x.^2+y.^2))/10;
S = [x;y;z];
h = surf(x,y,z)
set(h,'edgecolor','none')
colormap summer

% Number of points
n = 1000;

% Line to project on the surface with n values to get a feel for it...
t = linspace(-1,1,n);
height = sin(3*((.5*ones(1,n)).^2+ t.^2))/10;
P = [(.5*ones(1,n));t;height];

% Plotting the line on the surface:
hold on
%plot3(P(1,:),P(2,:),P(3,:),'o')

for j=1:2
% First midpoint iteration updates P...
P = [P(:,1), (P(:,1:end-1) + P(:,2:end))/2, P(:,end)];
%plot3(P(1,:), P(2,:), P(3,:), '.', 'MarkerSize', 20)

A = zeros(3,size(P,2));
for i = 1:size(P,2)
% Starting point will be the vertical projection of the first mid-point:
    A(:,i) = [P(1,i), P(2,i), sin(3*(P(1,i)^2+ P(2,i)^2))/10];
end

% Linear constraints:
nonlincon = @nonlincon;

% Placing fmincon in a loop for all the points

for i = 1:(size(A,2))
    % Objective function:
    objective = @(x)(P(1,i) - x(1))^2 + (P(2,i) - x(2))^2 + (P(3,i)-x(3))^2;
    A(:,i) = fmincon(objective, A(:,i), [], [], [], [], [], [], nonlincon);
end

P = A;
end

plot3(P(1,:), P(2,:), P(3,:), '.', 'MarkerSize',5,'Color','r')

with the nonlinear constraint in nonlincon.m:

function[c,ceq] = nlcon(x)
   c   = [];
   ceq = x(3) - sin(3*(x(1)^2+ x(2)^2))/10;

One nagging concern is the possibility of overfitting to the curve with this method, and this latter plot is an example of it. So I adjusted the code to just select one beginning and one ending point, and allowing the iterative process to find the rest of the curve, which for 100 iterations seemed to be heading in the right direction:

The above examples seem to follow a linear projection on the XY plane, but fortunately this is not a fixed pattern, which would cast further doubt on the method. See for instance the hyperbolic paraboloid x^2 - y^2:

Notice that there are algorithms to advance or push geodesic lines along a surface f(x,y) with small increments determined by the starting points and the normal vector to the surface, as in here. Thanks to the work of Alvise Vianello looking into the JS in that simulation and his sharing in GitHub, I was able to turn that algorithm into Matlab code, generating this plot for the first example, f(x,y) = x^3 + y^2:

Here is the Matlab code:

x = linspace(-1,1);
y = linspace(-1,1);
[x,y] = meshgrid(x,y);
z = x.^3 + y.^2;
S = [x;y;z];
h = surf(x,y,z)
set(h,'edgecolor','none')
colormap('gray');
hold on

f = @(x,y) x.^3 + y.^2; % The actual surface

dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y

N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.

C = {'k','b','r','g','y','m','c',[.8 .2 .6],[.2,.8,.1],[0.3010 0.7450 0.9330],[0.9290 0.6940 0.1250],[0.8500 0.3250 0.0980]}; % Color scheme

for s = 1:11     % No. of lines to be plotted.
start = -5:5;    % Distributing the starting points of the lines.
y0 = start(s)/5; % Fitting the starting pts between -1 and 1 along y axis.
x0 = 1;          % Along x axis always starts at 1.
dx0 = 0;         % Initial differential increment along x
dy0 = 0.05;      % Initial differential increment along y
step_size = 0.000008; % Will determine the progression rate from pt to pt.
eta =  step_size / sqrt(dx0^2 + dy0^2); % Normalization.
eps = 0.0001;          % Epsilon
max_num_iter = 100000; % Number of dots in each line.

x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values

for i = 2:(max_num_iter - 1)  % Creating the geodesic:
            xt = x(i);        % Values at point t of x, y and the function:
            yt = y(i);
            ft = f(xt,yt);

            xtm1 = x(i - 1);  % Values at t minus 1 (prior point) for x,y,f
            ytm1 = y(i - 1);
            ftm1 = f(xtm1,ytm1);

            xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
            ysymp = yt + (yt - ytm1);
            fsymp = ft + (ft - ftm1);

            df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
            n = N(xt,yt);                % Normal vector at point t
            gamma = df * n(3);           % Scalar x change f x z value of N

            xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
            ytp1 = ysymp - gamma * n(2);

            x(i + 1) = xtp1;
            y(i + 1) = ytp1;
end

P = [x; y; f(x,y)]; % Compiling results into a matrix.

indices = find(abs(P(1,:)) < 1); % Avoiding lines overshooting surface.
P = P(:,indices);
indices = find(abs(P(2,:)) < 1);
P = P(:,indices);

    units = 15; % Deternines speed (smaller, faster)
    packet = floor(size(P,2)/units);
    P = P(:,1: packet * units);

  for k = 1:packet:(packet * units)
        hold on
        plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))),...
            '.', 'MarkerSize', 3.5,'color',C{s})
        drawnow
  end

end

And here is an earlier example from above, but now calculated differently, and with lines starting side by side, following geodesics (no point-to-point trajectory):

    x = linspace(-1,1);
    y = linspace(-1,1);
    [x,y] = meshgrid(x,y);
    z = sin(3*(x.^2+y.^2))/10;
    S = [x;y;z];
    h = surf(x,y,z)
    set(h,'edgecolor','none')
    colormap('gray');
    hold on

    f = @(x,y) sin(3*(x.^2+y.^2))/10; % The actual surface

    dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
    dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y

    N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.

    C = {'k','r','g','y','m','c',[.8 .2 .6],[.2,.8,.1],[0.3010 0.7450 0.9330],[0.7890 0.5040 0.1250],[0.9290 0.6940 0.1250],[0.8500 0.3250 0.0980]}; % Color scheme

    for s = 1:11     % No. of lines to be plotted.
    start = -5:5;    % Distributing the starting points of the lines.
    x0 = -start(s)/5; % Fitting the starting pts between -1 and 1 along y axis.
    y0 = -1;          % Along x axis always starts at 1.
    dx0 = 0;         % Initial differential increment along x
    dy0 = 0.05;      % Initial differential increment along y
    step_size = 0.00005; % Will determine the progression rate from pt to pt.
    eta =  step_size / sqrt(dx0^2 + dy0^2); % Normalization.
    eps = 0.0001;          % Epsilon
    max_num_iter = 100000; % Number of dots in each line.

    x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
    y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values

    for i = 2:(max_num_iter - 1)  % Creating the geodesic:
                xt = x(i);        % Values at point t of x, y and the function:
                yt = y(i);
                ft = f(xt,yt);

                xtm1 = x(i - 1);  % Values at t minus 1 (prior point) for x,y,f
                ytm1 = y(i - 1);
                ftm1 = f(xtm1,ytm1);

                xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
                ysymp = yt + (yt - ytm1);
                fsymp = ft + (ft - ftm1);

                df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
                n = N(xt,yt);                % Normal vector at point t
                gamma = df * n(3);           % Scalar x change f x z value of N

                xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
                ytp1 = ysymp - gamma * n(2);

                x(i + 1) = xtp1;
                y(i + 1) = ytp1;
    end

    P = [x; y; f(x,y)]; % Compiling results into a matrix.

    indices = find(abs(P(1,:)) < 1); % Avoiding lines overshooting surface.
    P = P(:,indices);
    indices = find(abs(P(2,:)) < 1);
    P = P(:,indices);
    units = 35; % Deternines speed (smaller, faster)
    packet = floor(size(P,2)/units);
    P = P(:,1: packet * units);

  for k = 1:packet:(packet * units)
        hold on

        plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))), '.', 'MarkerSize', 5,'color',C{s})
        drawnow
  end

    end

Some more examples:

    x = linspace(-1,1);
    y = linspace(-1,1);
    [x,y] = meshgrid(x,y);
    z = x.^2 - y.^2;
    S = [x;y;z];
    h = surf(x,y,z)
    set(h,'edgecolor','none')
    colormap('gray');

    f = @(x,y) x.^2 - y.^2; % The actual surface

    dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
    dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y

    N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.

    C = {'b','w','r','g','y','m','c',[0.75, 0.75, 0],[0.9290, 0.6940, 0.1250],[0.3010 0.7450 0.9330],[0.1290 0.6940 0.1250],[0.8500 0.3250 0.0980]}; % Color scheme

    for s = 1:11     % No. of lines to be plotted.
    start = -5:5;    % Distributing the starting points of the lines.
    x0 = -start(s)/5; % Fitting the starting pts between -1 and 1 along y axis.
    y0 = -1;          % Along x axis always starts at 1.
    dx0 = 0;         % Initial differential increment along x
    dy0 = 0.05;      % Initial differential increment along y
    step_size = 0.00005; % Will determine the progression rate from pt to pt.
    eta =  step_size / sqrt(dx0^2 + dy0^2); % Normalization.
    eps = 0.0001;          % Epsilon
    max_num_iter = 100000; % Number of dots in each line.

    x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
    y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values

    for i = 2:(max_num_iter - 1)  % Creating the geodesic:
                xt = x(i);        % Values at point t of x, y and the function:
                yt = y(i);
                ft = f(xt,yt);

                xtm1 = x(i - 1);  % Values at t minus 1 (prior point) for x,y,f
                ytm1 = y(i - 1);
                ftm1 = f(xtm1,ytm1);

                xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
                ysymp = yt + (yt - ytm1);
                fsymp = ft + (ft - ftm1);

                df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
                n = N(xt,yt);                % Normal vector at point t
                gamma = df * n(3);           % Scalar x change f x z value of N

                xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
                ytp1 = ysymp - gamma * n(2);

                x(i + 1) = xtp1;
                y(i + 1) = ytp1;
    end

    P = [x; y; f(x,y)]; % Compiling results into a matrix.

    indices = find(abs(P(1,:)) < 1); % Avoiding lines overshooting surface.
    P = P(:,indices);
    indices = find(abs(P(2,:)) < 1);
    P = P(:,indices);
    units = 45; % Deternines speed (smaller, faster)
    packet = floor(size(P,2)/units);
    P = P(:,1: packet * units);

  for k = 1:packet:(packet * units)
        hold on
        plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))), '.', 'MarkerSize', 5,'color',C{s})
        drawnow
  end

  end

Or this one:

    x = linspace(-1,1);
    y = linspace(-1,1);
    [x,y] = meshgrid(x,y);
    z = .07 * (.1 + x.^2 + y.^2).^(-1);
    S = [x;y;z];
    h = surf(x,y,z)
    zlim([0 8])
    set(h,'edgecolor','none')
    colormap('gray');
    axis off
    hold on

    f = @(x,y) .07 * (.1 + x.^2 + y.^2).^(-1);    % The actual surface

    dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
    dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y

    N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.

     C = {'w',[0.8500, 0.3250, 0.0980],[0.9290, 0.6940, 0.1250],'g','y','m','c',[0.75, 0.75, 0],'r',...
         [0.56,0,0.85],'m'}; % Color scheme

    for s = 1:10     % No. of lines to be plotted.
    start = -9:2:9;
    x0 = -start(s)/10;
    y0 = -1;          % Along x axis always starts at 1.
    dx0 = 0;         % Initial differential increment along x
    dy0 = 0.05;      % Initial differential increment along y
    step_size = 0.00005; % Will determine the progression rate from pt to pt.
    eta =  step_size / sqrt(dx0^2 + dy0^2); % Normalization.
    eps = 0.0001;          % EpsilonA
    max_num_iter = 500000; % Number of dots in each line.

    x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
    y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values

    for i = 2:(max_num_iter - 1)  % Creating the geodesic:
                xt = x(i);        % Values at point t of x, y and the function:
                yt = y(i);
                ft = f(xt,yt);

                xtm1 = x(i - 1);  % Values at t minus 1 (prior point) for x,y,f
                ytm1 = y(i - 1);
                ftm1 = f(xtm1,ytm1);

                xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
                ysymp = yt + (yt - ytm1);
                fsymp = ft + (ft - ftm1);

                df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
                n = N(xt,yt);                % Normal vector at point t
                gamma = df * n(3);           % Scalar x change f x z value of N

                xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
                ytp1 = ysymp - gamma * n(2);

                x(i + 1) = xtp1;
                y(i + 1) = ytp1;
    end

     P = [x; y; f(x,y)]; % Compiling results into a matrix.

    indices = find(abs(P(1,:)) < 1.5); % Avoiding lines overshooting surface.
    P = P(:,indices);
    indices = find(abs(P(2,:)) < 1);
    P = P(:,indices);

    units = 15; % Deternines speed (smaller, faster)
    packet = floor(size(P,2)/units);
    P = P(:,1: packet * units);

  for k = 1:packet:(packet * units)
        hold on
        plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))),...
            '.', 'MarkerSize', 3.5,'color',C{s})
        drawnow
  end

    end

Or a sinc function:

    x = linspace(-10, 10);
    y = linspace(-10, 10);
    [x,y] = meshgrid(x,y);
    z = sin(1.3*sqrt (x.^ 2 + y.^ 2) + eps)./ (sqrt (x.^ 2 + y.^ 2) + eps);
    S = [x;y;z];
    h = surf(x,y,z)
    set(h,'edgecolor','none')
    colormap('gray');
    axis off
    hold on

    f = @(x,y) sin(1.3*sqrt (x.^ 2 + y.^ 2) + eps)./ (sqrt (x.^ 2 + y.^ 2) + eps);   % The actual surface

    dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
    dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y

    N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.

    C = {'w',[0.8500, 0.3250, 0.0980],[0.9290, 0.6940, 0.1250],'g','y','r','c','m','w',...
         [0.56,0,0.85],[0.8500, 0.7250, 0.0980],[0.2290, 0.1940, 0.6250],'w',...
         [0.890, 0.1940, 0.4250],'y',[0.2290, 0.9940, 0.3250],'w',[0.1500, 0.7250, 0.0980],...
         [0.8500, 0.3250, 0.0980],'m','w'}; % Color scheme

    for s = 1:12     % No. of lines to be plotted.

    x0 = 10;
    y0 = 10;          % Along x axis always starts at 1.
    dx0 = -0.001*(cos(pi /2 *s/11));         % Initial differential increment along x
    dy0 = -0.001*(sin(pi /2 *s/11));         % Initial differential increment along y
    step_size = 0.0005; % Will determine the progression rate from pt to pt.
    % Making it smaller increases the length of the curve.
    eta =  step_size / sqrt(dx0^2 + dy0^2); % Normalization.
    eps = 0.0001;          % EpsilonA
    max_num_iter = 500000; % Number of dots in each line.

    x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
    y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values

    for i = 2:(max_num_iter - 1)  % Creating the geodesic:
                xt = x(i);        % Values at point t of x, y and the function:
                yt = y(i);
                ft = f(xt,yt);

                xtm1 = x(i - 1);  % Values at t minus 1 (prior point) for x,y,f
                ytm1 = y(i - 1);
                ftm1 = f(xtm1,ytm1);

                xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
                ysymp = yt + (yt - ytm1);
                fsymp = ft + (ft - ftm1);

                df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
                n = N(xt,yt);                % Normal vector at point t
                gamma = df * n(3);           % Scalar x change f x z value of N

                xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
                ytp1 = ysymp - gamma * n(2);

                x(i + 1) = xtp1;
                y(i + 1) = ytp1;
    end

     P = [x; y; f(x,y)]; % Compiling results into a matrix.

    indices = find(abs(P(1,:)) < 10); % Avoiding lines overshooting surface.
    P = P(:,indices);
    indices = find(abs(P(2,:)) < 10);
    P = P(:,indices);

    units = 15; % Deternines speed (smaller, faster)
    packet = floor(size(P,2)/units);
    P = P(:,1: packet * units);

  for k = 1:packet:(packet * units)
        hold on
        plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))),...
            '.', 'MarkerSize', 3.5,'color',C{s})
        drawnow
  end

    end

And one very last one:

    x = linspace(-1.5,1.5);
    y = linspace(-1,1);
    [x,y] = meshgrid(x,y);
    z = 0.5 *y.*sin(5 * x) - 0.5 * x.*cos(5 * y)+1.5;
    S = [x;y;z];
    h = surf(x,y,z)
    zlim([0 8])
    set(h,'edgecolor','none')
    colormap('gray');
    axis off
    hold on

    f = @(x,y) 0.5 *y.* sin(5 * x) - 0.5 * x.*cos(5 * y)+1.5;     % The actual surface

    dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
    dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y

    N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.

     C = {'w',[0.8500, 0.3250, 0.0980],[0.9290, 0.6940, 0.1250],'g','y','k','c',[0.75, 0.75, 0],'r',...
         [0.56,0,0.85],'m'}; % Color scheme

    for s = 1:11     % No. of lines to be plotted.
    start = [0, 0.7835,  -0.7835, 0.5877, -0.5877, 0.3918, -0.3918, 0.1959, -0.1959, 0.9794, -0.9794];
    x0 = start(s);
    y0 = -1;          % Along x axis always starts at 1.
    dx0 = 0;         % Initial differential increment along x
    dy0 = 0.05;      % Initial differential increment along y
    step_size = 0.00005; % Will determine the progression rate from pt to pt.
    % Making it smaller increases the length of the curve.
    eta =  step_size / sqrt(dx0^2 + dy0^2); % Normalization.
    eps = 0.0001;          % EpsilonA
    max_num_iter = 500000; % Number of dots in each line.

    x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
    y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values

    for i = 2:(max_num_iter - 1)  % Creating the geodesic:
                xt = x(i);        % Values at point t of x, y and the function:
                yt = y(i);
                ft = f(xt,yt);

                xtm1 = x(i - 1);  % Values at t minus 1 (prior point) for x,y,f
                ytm1 = y(i - 1);
                ftm1 = f(xtm1,ytm1);

                xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
                ysymp = yt + (yt - ytm1);
                fsymp = ft + (ft - ftm1);

                df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
                n = N(xt,yt);                % Normal vector at point t
                gamma = df * n(3);           % Scalar x change f x z value of N

                xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
                ytp1 = ysymp - gamma * n(2);

                x(i + 1) = xtp1;
                y(i + 1) = ytp1;
    end

     P = [x; y; f(x,y)]; % Compiling results into a matrix.

    indices = find(abs(P(1,:)) < 1.5); % Avoiding lines overshooting surface.
    P = P(:,indices);
    indices = find(abs(P(2,:)) < 1);
    P = P(:,indices);

    units = 15; % Deternines speed (smaller, faster)
    packet = floor(size(P,2)/units);
    P = P(:,1: packet * units);

  for k = 1:packet:(packet * units)
        hold on
        plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))),...
            '.', 'MarkerSize', 3.5,'color',C{s})
        drawnow
  end

    end

On a torus this becomes more complicated, and so far I haven't been able to adapt this code to it. As a reference, the code provided by Paul Chesler in here does provide with an option. It involves saving a file tor.m as

function xp=tor(t,x)
xp=zeros(4,1);
xp(1)=x(2);
xp(2)=-(2+cos(x(1)))*sin(x(1))*x(4)^2;
xp(3)=x(4);
xp(4)=2*(sin(x(1))/(2+cos(x(1))))*x(2)*x(4);

and in another file torus.m plot the torus and call the function above to solve a system of second order differential equations. Here it is for different geodesics:

from below:

with torus.m:

[u,v]=meshgrid(linspace(0,2*pi,100),linspace(0,2*pi,100));
x=(2+cos(u)).*cos(v);
y=(2+cos(u)).*sin(v);
z=sin(u);

h=surf(x,y,z);
set(h,'edgecolor','none');
colormap('gray');
daspect([1 1 1])

tspan = linspace(pi,100*pi,1000000);
[t,X]=ode45('tor',tspan ,[pi,.1,-pi/2,.2]);
u=X(:,1);
v=X(:,3);
x=(2+cos(u)).*cos(v);
y=(2+cos(u)).*sin(v);
z=sin(u);

P=[x y z];
disp(size(P))

units = 10; % Deternines speed (smaller, faster)
packet = floor(size(P,1)/units);

  for k = 1:packet:(packet * units)
        hold on
        plot3(P(k:(k+packet-1),1), P((k:(k+packet-1)),2), P((k:(k+packet-1)),3),...
            '.-', 'MarkerSize', 3.5,'color','r', 'LineWidth', 3)
        drawnow
        pause(2)
  end

or

[u,v]=meshgrid(linspace(0,2*pi,100),linspace(0,2*pi,100));
x=(2+cos(u)).*cos(v);
y=(2+cos(u)).*sin(v);
z=sin(u);

h=surf(x,y,z);
set(h,'edgecolor','none');
colormap('gray');
daspect([1 1 1])

tspan = linspace(pi,100*pi,1000000);
[t,X]=ode45('tor',tspan ,[pi/6,.1,-pi/2,.2]);
u=X(:,1);
v=X(:,3);
x=(2+cos(u)).*cos(v);
y=(2+cos(u)).*sin(v);
z=sin(u);

P=[x y z];
disp(size(P))

units = 10; % Deternines speed (smaller, faster)
packet = floor(size(P,1)/units);

  for k = 1:packet:(packet * units)
        hold on
        plot3(P(k:(k+packet-1),1), P((k:(k+packet-1)),2), P((k:(k+packet-1)),3),...
            '.-', 'MarkerSize', 3.5,'color','m', 'LineWidth', 3)
        drawnow
        pause(2)
  end

or

from below:

tspan = linspace(pi,100*pi,1000000); defines the limits of integration along 1e^6 pts, and [3/4 * pi,.1,-pi/2,.2], the parametric values for a starting point (on the last example).

On a sphere:

code here.


来源: stackexchange

问:
Why does include statement return an error that basically states that the language (C++) has a bug?
Some info: I am a person who has gotten into C++ development recently, which means i really don't know how to troubleshoot things like this. You have been warned.

Question: See title

I believe I have what looks to me to be a compiler error, but have not been able to find anything like it anywhere on stack overflow.

Problematic code:

```<code>#include <map> // This is the line that triggered the error.

using namespace std;

class studySet
{

public:

    studySet(int t, double v, string n, string a,map<int, string[2]> s)
    {

        type = t;
        version = v;
        name = n;
        author = a;
        setTerms = s;

    }
    int getType()
    {

        return type;

    }
    double getVersion()
    {

        return version;

    }
    string getName()
    {

        return name;

    }
    string getAuthor()
    {

        return author;

    }

    void study()
    {

    }

private:

    int type;
    double version;
    string name;
    string author;
    map<int, string[2]> setTerms;

};
</code>```

Error shown:

```<code>No viable overloaded '=' //Shown in a file included from above code.
</code>```

OS: MacOS 10.14.3 Mojave

IDE: Xcode

Compiler: IDK, whatever Xcode uses.

Please comment if you need more code from me as i am fairly new to stack overflow as well.

答:
Thanks to the comments, I was able to determine that the line

```<code>setTerms = s;
</code>```

triggered the error. That is because you can't pass a vector to fill another vector.

ContentLicense: CC BY-SA 4.0

Body:
Some info: I am a person who has gotten into C++ development recently, which means i really don't know how to troubleshoot things like this. You have been warned.

Question: See title

I believe I have what looks to me to be a compiler error, but have not been able to find anything like it anywhere on stack overflow.

Problematic code:

#include <map> // This is the line that triggered the error.

using namespace std;

class studySet
{

public:

    studySet(int t, double v, string n, string a,map<int, string[2]> s)
    {

        type = t;
        version = v;
        name = n;
        author = a;
        setTerms = s;

    }
    int getType()
    {

        return type;

    }
    double getVersion()
    {

        return version;

    }
    string getName()
    {

        return name;

    }
    string getAuthor()
    {

        return author;

    }

    void study()
    {

    }

private:

    int type;
    double version;
    string name;
    string author;
    map<int, string[2]> setTerms;

};

Error shown:

No viable overloaded '=' //Shown in a file included from above code.

OS: MacOS 10.14.3 Mojave

IDE: Xcode

Compiler: IDK, whatever Xcode uses.

Please comment if you need more code from me as i am fairly new to stack overflow as well.

Body:
Thanks to the comments, I was able to determine that the line

setTerms = s;

triggered the error. That is because you can't pass a vector to fill another vector.


问:
Pandas - Fast way of accessing a column of objects' attribute
Let's say I have a custom class in python, that has the attribute val. If I have a pandas dataframe with a column of these objects, how can I access this attribute and make a new column with this value?

Example data:

```<code>df
Out[46]:
row   custom_object
1     foo1
2     foo2
3     foo3
4     foo4
Name: book, dtype: object
</code>```

Where the custom objects are of class Foo:

```<code>class Foo:
    def __init__(self, val):
        self.val = val
</code>```

The only way I know of to create a new column with the instance attributes is using an apply and lambda combo which is slow on large datasets:

df['custom_val'] = df['custom_object'].apply(lambda x: x.val)

Is there a more efficient way?

答:
Setup code:

import operator
import random
from dataclasses import dataclass

import numpy as np
import pandas as pd

@dataclass
class SomeObj:
    val: int

df = pd.DataFrame(data={f"col_1": [SomeObj(random.randint(0, 10000)) for _ in range(10000000)]})
```
<hr />
<h2>Solution 1</h2>
<pre class="lang-py prettyprint-override"><code>df['col_1'].map(lambda elem: elem.val)
</code>```

Time: ~3.2 seconds

Solution 2

df['col_1'].map(operator.attrgetter('val'))
```
<p>Time: <strong>~2.7 seconds</strong></p>
<h2>Solution 3</h2>
<pre class="lang-py prettyprint-override"><code>[elem.val for elem in df['col_1']]
</code>```

Time: ~1.4 seconds

Note: Keep in mind that this solution produces a different result type, which may be an issue in certain situations.

ContentLicense: CC BY-SA 3.0

Body:
Let's say I have a custom class in python, that has the attribute val. If I have a pandas dataframe with a column of these objects, how can I access this attribute and make a new column with this value?

Example data:

df
Out[46]:
row   custom_object
1     foo1
2     foo2
3     foo3
4     foo4
Name: book, dtype: object

Where the custom objects are of class Foo:

class Foo:
    def __init__(self, val):
        self.val = val

The only way I know of to create a new column with the instance attributes is using an apply and lambda combo which is slow on large datasets:

df['custom_val'] = df['custom_object'].apply(lambda x: x.val)

Is there a more efficient way?

Body:
Setup code:

import operator
import random
from dataclasses import dataclass

import numpy as np
import pandas as pd

@dataclass
class SomeObj:
    val: int

df = pd.DataFrame(data={f"col_1": [SomeObj(random.randint(0, 10000)) for _ in range(10000000)]})

Solution 1

df['col_1'].map(lambda elem: elem.val)

Time: ~3.2 seconds

Solution 2

df['col_1'].map(operator.attrgetter('val'))

Time: ~2.7 seconds

Solution 3

[elem.val for elem in df['col_1']]

Time: ~1.4 seconds

Note: Keep in mind that this solution produces a different result type, which may be an issue in certain situations.


问:
How can I use the FTS5 extension with the sqlite3 python module with Python 3.7?
How can I use the FTS5 extension with the sqlite3 python module with Python 3.7?

I tried to run the following code in Python with python testFTS5.py:

```<code>import sqlite3
conn = sqlite.connect('some_db.db')
sqlite.enable_load_extension(True)
sqlite.load_extension('fts5')
</code>```

which results in the error message:

```<code>Traceback (most recent call last):
  File "./src/test.py", line 3, in <module>
    sqlite.enable_load_extension(True)
AttributeError: module 'sqlite3' has no attribute 'enable_load_extension'
</code>```

I tried sqlite.load_extension('FTS5') and sqlite.load_extension('ENABLE_FTS5') but it unsurprisingly yields the same error message (with the corresponding filename being not found). I also tried running the code with LD_LIBRARY_PATH=/usr/local/bin python testFTS5.py but I get the same error message.

I checked the sqlite3 location by running the following code in the terminal:

```<code>derno@ompn:/mnt/ilcompn0d1/user/dernonco/fts-test$ which sqlite3
/usr/local/bin/sqlite3
</code>```

and I listed the installed sqlite3 extensions:

```<code>derno@ompn:/mnt/ilcompn0d1/user/dernonco/fts-test$ sqlite3
SQLite version 3.18.0 2017-03-28 18:48:43
Enter ".help" for usage hints.
Connected to a transient in-memory database.
Use ".open FILENAME" to reopen on a persistent database.
sqlite> pragma compile_options;
COMPILER=gcc-5.4.0 20160609
DEFAULT_SYNCHRONOUS=2
DEFAULT_WAL_SYNCHRONOUS=2
ENABLE_FTS5
ENABLE_RTREE
SYSTEM_MALLOC
THREADSAFE=1
</code>```

which seems to indicate that the FTS5 is available in my /usr/local/bin/sqlite3 version.

However when I run

```<code>import sqlite3

con = sqlite3.connect(':memory:')
cur = con.cursor()
cur.execute('pragma compile_options;')
available_pragmas = cur.fetchall()
con.close()

print(available_pragmas)
</code>```

it outputs:

```<code>[('COMPILER=gcc-5.4.0 20160609',), ('DEFAULT_SYNCHRONOUS=2',), ('DEFAULT_WAL_SYNCHRONOUS=2',),
('ENABLE_FTS3',), ('ENABLE_RTREE',), ('SYSTEM_MALLOC',), ('THREADSAFE=1',)]
</code>```

There is no ENABLE_FTS5 in that list.

I tried with Python 3.7.6 (default, Dec 19 2019, 23:49:42)  and  Python 3.6.7 (default, Oct 25 2018, 09:16:13).

答:
You would call .enable_load_extension(True) and .load_extension('fts5') on the connection object, and not on the module.

However, this shouldn't be necessary since -- as you've seen -- your installation supports full-text search.

Here's a way you can test that out just to be sure:

```<code>import sqlite3

conn = sqlite3.connect(':memory:')
conn.execute("""create virtual table fts5test using fts5 (data);""")
conn.execute("""insert into fts5test (data)
                values ('this is a test of full-text search');""")
conn.execute("""select * from fts5test where data match 'full';""").fetchall()
</code>```

Result:

```<code>In [67]: conn.execute("""select * from fts5test where data match 'full';""").fetchall()
Out[67]: [('this is a test of full-text search',)]
</code>```

Body:
How can I use the FTS5 extension with the sqlite3 python module with Python 3.7?

I tried to run the following code in Python with python testFTS5.py:

import sqlite3
conn = sqlite.connect('some_db.db')
sqlite.enable_load_extension(True)
sqlite.load_extension('fts5')

which results in the error message:

Traceback (most recent call last):
  File "./src/test.py", line 3, in <module>
    sqlite.enable_load_extension(True)
AttributeError: module 'sqlite3' has no attribute 'enable_load_extension'

I tried sqlite.load_extension('FTS5') and sqlite.load_extension('ENABLE_FTS5') but it unsurprisingly yields the same error message (with the corresponding filename being not found). I also tried running the code with LD_LIBRARY_PATH=/usr/local/bin python testFTS5.py but I get the same error message.

I checked the sqlite3 location by running the following code in the terminal:

derno@ompn:/mnt/ilcompn0d1/user/dernonco/fts-test$ which sqlite3
/usr/local/bin/sqlite3

and I listed the installed sqlite3 extensions:

derno@ompn:/mnt/ilcompn0d1/user/dernonco/fts-test$ sqlite3
SQLite version 3.18.0 2017-03-28 18:48:43
Enter ".help" for usage hints.
Connected to a transient in-memory database.
Use ".open FILENAME" to reopen on a persistent database.
sqlite> pragma compile_options;
COMPILER=gcc-5.4.0 20160609
DEFAULT_SYNCHRONOUS=2
DEFAULT_WAL_SYNCHRONOUS=2
ENABLE_FTS5
ENABLE_RTREE
SYSTEM_MALLOC
THREADSAFE=1

which seems to indicate that the FTS5 is available in my /usr/local/bin/sqlite3 version.

However when I run

import sqlite3

con = sqlite3.connect(':memory:')
cur = con.cursor()
cur.execute('pragma compile_options;')
available_pragmas = cur.fetchall()
con.close()

print(available_pragmas)

it outputs:

[('COMPILER=gcc-5.4.0 20160609',), ('DEFAULT_SYNCHRONOUS=2',), ('DEFAULT_WAL_SYNCHRONOUS=2',),
('ENABLE_FTS3',), ('ENABLE_RTREE',), ('SYSTEM_MALLOC',), ('THREADSAFE=1',)]

There is no ENABLE_FTS5 in that list.

I tried with Python 3.7.6 (default, Dec 19 2019, 23:49:42)  and  Python 3.6.7 (default, Oct 25 2018, 09:16:13).

Body:
You would call .enable_load_extension(True) and .load_extension('fts5') on the connection object, and not on the module.

However, this shouldn't be necessary since -- as you've seen -- your installation supports full-text search.

Here's a way you can test that out just to be sure:

import sqlite3

conn = sqlite3.connect(':memory:')
conn.execute("""create virtual table fts5test using fts5 (data);""")
conn.execute("""insert into fts5test (data)
                values ('this is a test of full-text search');""")
conn.execute("""select * from fts5test where data match 'full';""").fetchall()

Result:

In [67]: conn.execute("""select * from fts5test where data match 'full';""").fetchall()
Out[67]: [('this is a test of full-text search',)]


问:
Laravel add calculate the sum where the value is in another table
Basically i have two tables. answers and choices . In my choices table, I got a column choices.value which is has values 0-4. My current query like this

```<code> $answers = \DB::table('answers')
                    ->join('choices','answers.choice_id','=','choices.choice_id')
                    ->select('answers.user_id','choices.choice_id','choices.choice','choices.value')
                    ->where('answers.user_id',\Auth::user()->id)
                    //->groupBy('answers.user_id')
                    ->get();
</code>```

My current response is like this

```<code>    "user_id": 2,
        "choice_id": 2,
        "choice": "I feel discourated about the future",
        "value": 1
    },
    {
        "user_id": 2,
        "choice_id": 2,
        "choice": "I don't enjoy things the way I used to",
        "value": 1
    },
    {
        "user_id": 2,
        "choice_id": 2,
        "choice": "I feel guilty a good part of time",
        "value": 1
</code>```

how do i add the values so that my result will be like this

```<code>"user_id":2,
"total_score":3
</code>```

I tried doing DB::raw(SUM(choices.values)as score) but i get a large amount. I guess its adding all the choices values from the choices table and not in the answers.

My answers db which i only select the answers of the user = 2. I limit to 5

```<code>+---------+-------------+-----------+
| user_id | question_id | choice_id |
+---------+-------------+-----------+
|       2 |           1 |         2 |
|       2 |           2 |         2 |
|       2 |           3 |         2 |
|       2 |           4 |         2 |
|       2 |           5 |         2 |
+---------+-------------+-----------+
</code>```

My choices table, I only select questions 1 and 2 and their choices.

```<code>+-----------+-------------+--------------------------------------------------------------+-------+
| choice_id | question_id | choice                                                       | value |
+-----------+-------------+--------------------------------------------------------------+-------+
|         1 |           1 | I do not feel sad                                            |     0 |
|         2 |           1 | I feel sad                                                   |     1 |
|         3 |           1 | I am sad all the time and I can't snap out of it             |     2 |
|         4 |           1 | I am so sad and unhappy that I can't stand it                |     3 |
|         1 |           2 | I am not particularly discouraged about the future           |     0 |
|         2 |           2 | I feel discourated about the future                          |     1 |
|         3 |           2 | I feel I have nothing to look forward to                     |     2 |
|         4 |           2 | I feel the future is hopeless and that things cannot improve |     3 |
+-----------+-------------+--------------------------------------------------------------+-------+
</code>```

Also i want to make a new table named scores and the columns will be the result of what i want. I want to add the choices.values in the answer in every answers.user_id so that when i view the scores table, it will display the total score in every user or when the user finished answering all 21 questions because i only got 21 items it will automatically add in the scores table. Could I possibly do that?. Is it okay to add a value in the answers table based in the choice_id? That is what im thinking but I think its redundant since the choice_id is there already. Thanks in advance.

PS: Tried writing these queries but always get 441 which is the total value of all the choices in the choices table

```<code>SELECT  answers.user_id,choices.choice_id,choices.value,COALESCE(sum(choices.value),0) as score  FROM
`answers`  JOIN `choices` ON  choices.choice_id = answers.choice_id where
answers.user_id = 2

SELECT answers.user_id,choices.choice_id,choices.value,SUM(choices.value),0 as score FROM `answers`
join choices on choices.choice_id = answers.choice_id
where answers.user_id = 2

SELECT answers.user_id,choices.choice_id, sum(choices.value) from answers
JOIN `choices` ON  choices.choice_id = answers.choice_id
group by answers.user_id
</code>```

答:
use group by maybe it will help</p>

```<code> $answers = \DB::table('answers')
                    ->join('choices','answers.choice_id','=','choices.choice_id')
                    ->select('answers.user_id',DB::raw('count(*) as total'))
                    ->where('answers.user_id',\Auth::user()->id)
                    ->groupBy('answers.user_id')
                    ->get();
</code>```

Body:
Basically i have two tables. answers and choices . In my choices table, I got a column choices.value which is has values 0-4. My current query like this

 $answers = \DB::table('answers')
                    ->join('choices','answers.choice_id','=','choices.choice_id')
                    ->select('answers.user_id','choices.choice_id','choices.choice','choices.value')
                    ->where('answers.user_id',\Auth::user()->id)
                    //->groupBy('answers.user_id')
                    ->get();

My current response is like this

    "user_id": 2,
        "choice_id": 2,
        "choice": "I feel discourated about the future",
        "value": 1
    },
    {
        "user_id": 2,
        "choice_id": 2,
        "choice": "I don't enjoy things the way I used to",
        "value": 1
    },
    {
        "user_id": 2,
        "choice_id": 2,
        "choice": "I feel guilty a good part of time",
        "value": 1

how do i add the values so that my result will be like this

"user_id":2,
"total_score":3

I tried doing DB::raw(SUM(choices.values)as score) but i get a large amount. I guess its adding all the choices values from the choices table and not in the answers.

My answers db which i only select the answers of the user = 2. I limit to 5

+---------+-------------+-----------+
| user_id | question_id | choice_id |
+---------+-------------+-----------+
|       2 |           1 |         2 |
|       2 |           2 |         2 |
|       2 |           3 |         2 |
|       2 |           4 |         2 |
|       2 |           5 |         2 |
+---------+-------------+-----------+

My choices table, I only select questions 1 and 2 and their choices.

+-----------+-------------+--------------------------------------------------------------+-------+
| choice_id | question_id | choice                                                       | value |
+-----------+-------------+--------------------------------------------------------------+-------+
|         1 |           1 | I do not feel sad                                            |     0 |
|         2 |           1 | I feel sad                                                   |     1 |
|         3 |           1 | I am sad all the time and I can't snap out of it             |     2 |
|         4 |           1 | I am so sad and unhappy that I can't stand it                |     3 |
|         1 |           2 | I am not particularly discouraged about the future           |     0 |
|         2 |           2 | I feel discourated about the future                          |     1 |
|         3 |           2 | I feel I have nothing to look forward to                     |     2 |
|         4 |           2 | I feel the future is hopeless and that things cannot improve |     3 |
+-----------+-------------+--------------------------------------------------------------+-------+

Also i want to make a new table named scores and the columns will be the result of what i want. I want to add the choices.values in the answer in every answers.user_id so that when i view the scores table, it will display the total score in every user or when the user finished answering all 21 questions because i only got 21 items it will automatically add in the scores table. Could I possibly do that?. Is it okay to add a value in the answers table based in the choice_id? That is what im thinking but I think its redundant since the choice_id is there already. Thanks in advance.

PS: Tried writing these queries but always get 441 which is the total value of all the choices in the choices table

SELECT  answers.user_id,choices.choice_id,choices.value,COALESCE(sum(choices.value),0) as score  FROM
`answers`  JOIN `choices` ON  choices.choice_id = answers.choice_id where
answers.user_id = 2

SELECT answers.user_id,choices.choice_id,choices.value,SUM(choices.value),0 as score FROM `answers`
join choices on choices.choice_id = answers.choice_id
where answers.user_id = 2

SELECT answers.user_id,choices.choice_id, sum(choices.value) from answers
JOIN `choices` ON  choices.choice_id = answers.choice_id
group by answers.user_id

Body:
use group by maybe it will help

 $answers = \DB::table('answers')
                    ->join('choices','answers.choice_id','=','choices.choice_id')
                    ->select('answers.user_id',DB::raw('count(*) as total'))
                    ->where('answers.user_id',\Auth::user()->id)
                    ->groupBy('answers.user_id')
                    ->get();


问:
Django Python save database blog posts to a zipfile?
I have a django blog and want to download a backup zipfile with all the entries. The blog post text content is stored in the database.

I have written this code with the goal of trying to get the zipfile to save a bunch of .txt files in the main zip directory, but all this code does is outputs a single corrupted zip file. It cannot be unzipped but for some reason it can be opened in Word and it shows all of the blog post text mashed up.

```<code>def download_backups(request):
    zip_filename = "test.zip"
    s = BytesIO()
    zf = zipfile.ZipFile(s, "w")
    blogposts = Blog.objects.all()
    for blogpost in blogposts:
        filename = blogpost.title + ".txt"
        zf.writestr(filename, blogpost.content)
    resp = HttpResponse(s.getvalue())
    resp['Content-Disposition'] = 'attachment; filename=%s' % zip_filename
    return resp

</code>```

Any help is appreciated.

答:
Based on this answer to another question, you may be having an issue with the read mode. You'll also need to call zf.close(), either explicitly or implicitly, before the file will actually be complete.

I think there's a simpler way of handling this using a temporary file, which should have the advantage of not needing to fit all of the file's contents in memory.

```<code>from tempfile import TemporaryFile
from zipfile import ZipFile

with TemporaryFile() as tf:
    with ZipFile(tf, mode="w") as zf:
        zf.writestr("file1.txt", "The first file")
        zf.writestr("file2.txt", "A second file")
    tf.seek(0)
    print(tf.read())
</code>```

The with blocks here will result in your temp file going out of scope and being deleted, and zf.close being called implicitly before you attempt to read the file.

If the goal here is just to back up the data rather than using this specific format, though, I'd suggest using the built-in dumpdata management command. You can call it from code if you want to serve the results through a view like this.

Body:
I have a django blog and want to download a backup zipfile with all the entries. The blog post text content is stored in the database.

I have written this code with the goal of trying to get the zipfile to save a bunch of .txt files in the main zip directory, but all this code does is outputs a single corrupted zip file. It cannot be unzipped but for some reason it can be opened in Word and it shows all of the blog post text mashed up.

def download_backups(request):
    zip_filename = "test.zip"
    s = BytesIO()
    zf = zipfile.ZipFile(s, "w")
    blogposts = Blog.objects.all()
    for blogpost in blogposts:
        filename = blogpost.title + ".txt"
        zf.writestr(filename, blogpost.content)
    resp = HttpResponse(s.getvalue())
    resp['Content-Disposition'] = 'attachment; filename=%s' % zip_filename
    return resp

Any help is appreciated.

Body:
Based on this answer to another question, you may be having an issue with the read mode. You'll also need to call zf.close(), either explicitly or implicitly, before the file will actually be complete.

I think there's a simpler way of handling this using a temporary file, which should have the advantage of not needing to fit all of the file's contents in memory.

from tempfile import TemporaryFile
from zipfile import ZipFile

with TemporaryFile() as tf:
    with ZipFile(tf, mode="w") as zf:
        zf.writestr("file1.txt", "The first file")
        zf.writestr("file2.txt", "A second file")
    tf.seek(0)
    print(tf.read())

The with blocks here will result in your temp file going out of scope and being deleted, and zf.close being called implicitly before you attempt to read the file.

If the goal here is just to back up the data rather than using this specific format, though, I'd suggest using the built-in dumpdata management command. You can call it from code if you want to serve the results through a view like this.


问:
How to create a horizontal list view or recycler view with fixed size by changing the width of its items
I'm beginner in Android Development and trying to create a display of items in such a way that the size of these varieties according to the quantity, as they have done in this app. Any suggestions?

Screenshot 1

Screenshot 2

答:
```<code>LinearLayoutManager layoutManager = new LinearLayoutManager(this, LinearLayoutManager.HORIZONTAL, false);

RecyclerView myList = (RecyclerView) findViewById(R.id.my_recycler_view);

myList.setLayoutManager(layoutManager);
</code>```


问:
AttributeError: Everything was right from the tutorial video I used to make this, but still there is an error
I used a youtube video to make this platformer game. All of the attributes and parameters were defined, but still there appears an error. Here is my code block:

```<code>def update(self):
        # Game Loop Update
        self.all_sprites.update()

        # Check if player hits a platform only if falling
        if self.player.vel.y > 0:
            hits = pg.sprite.spritecollide(self.player, self.platforms, False)
            if hits:
                self.player.pos.y = hits[0].rect.top + 1
                self.player.vel.y = 0

        # If player reaches top 1/4 of screen
        if self.player.rect.top <= HEIGHT / 4:
            self.player.pos.y += abs(self.player.vel.y)
            for plat in self.platforms:
                plat.rect.y += abs(self.player.vel.y)
                if plat.rect.top >= HEIGHT:
                    plat.kill()

</code>```

This is the Player Class for extra info:

```<code>class Player(pg.sprite.Sprite):
    def __init__(self, game):
        pg.sprite.Sprite.__init__(self)
        self.game = game
        self.image = pg.Surface((30,40))
        self.image.fill(YELLOW)
        self.rect = self.image.get_rect()
        self.rect.center = (WIDTH/2, HEIGHT/2)
        self.pos = vec(WIDTH/2, HEIGHT/2)
        self.vel = vec(0,0)
        self.acc = vec(0,0)
</code>```

In the 6th line given in the code, the console gives error:
AttributeError: 'bool' object has no attribute 'vel'
Can anybody help?
Tell me if I need to provide more info.

答:
According to the error message, self.player contains a boolean value, not an object with a vel attribute as you expect.

The issue is not in the code you have provided. I suggest you review the surrounding code to ensure that this attribute is being handled properly.

Body:
I used a youtube video to make this platformer game. All of the attributes and parameters were defined, but still there appears an error. Here is my code block:

def update(self):
        # Game Loop Update
        self.all_sprites.update()

        # Check if player hits a platform only if falling
        if self.player.vel.y > 0:
            hits = pg.sprite.spritecollide(self.player, self.platforms, False)
            if hits:
                self.player.pos.y = hits[0].rect.top + 1
                self.player.vel.y = 0

        # If player reaches top 1/4 of screen
        if self.player.rect.top <= HEIGHT / 4:
            self.player.pos.y += abs(self.player.vel.y)
            for plat in self.platforms:
                plat.rect.y += abs(self.player.vel.y)
                if plat.rect.top >= HEIGHT:
                    plat.kill()

This is the Player Class for extra info:

class Player(pg.sprite.Sprite):
    def __init__(self, game):
        pg.sprite.Sprite.__init__(self)
        self.game = game
        self.image = pg.Surface((30,40))
        self.image.fill(YELLOW)
        self.rect = self.image.get_rect()
        self.rect.center = (WIDTH/2, HEIGHT/2)
        self.pos = vec(WIDTH/2, HEIGHT/2)
        self.vel = vec(0,0)
        self.acc = vec(0,0)

In the 6th line given in the code, the console gives error:
AttributeError: 'bool' object has no attribute 'vel'
Can anybody help?
Tell me if I need to provide more info.


问:
"sna" or "igraph" : Why do I get different degree values for undirected graph?
I am doing some basic network analysis using networks from the R package "networkdata". To this end, I use the package "igraph" as well as "sna". However, I realised that the results of descriptive network statistics vary depending on the package I use. Most variation is not too grave but the average degree of my undirected graph halved as soon as I switched from "sna" to "igraph".

```<code>library(networkdata)
n_1 <- covert_28

library(igraph)
library(sna)

n_1_adjmat <- as_adjacency_matrix(n_1)
n_1_adjmat2 <- as.matrix(n_1_adjmat)

mean(sna::degree(n_1_adjmat2, cmode = "freeman")) # [1] 23.33333
mean(igraph::degree(n_1, mode = "all")) # [1] 11.66667

</code>```

This doesn't happen in case of my directed graph. Here, I get the same results regardless of using "sna" or "igraph".

Is there any explanation for this phenomenon? And if so, is there anything I can do in order to prevent this from happening?

Thank you in advance!

答:
This is explained in the documentation for sna::degree.

indegree of a vertex, v, corresponds to the cardinality
  of the vertex set N^+(v) = {i in V(G) : (i,v) in E(G)};
  outdegree corresponds to the cardinality of the vertex
  set N^-(v) = {i in V(G) : (v,i) in E(G)}; and total
  (or “Freeman”) degree corresponds to |N^+(v)| + |N^-(v)|.

(Note that, for simple graphs,

indegree=outdegree=total degree/2.)

A simpler example than yours makes it clear.

```<code>library(igraph)
library(sna)

g = make_ring(3)
plot(g)
</code>```

```<code>AM = as.matrix(as_adjacency_matrix(g))
sna::degree(AM)
[1] 4 4 4

igraph::degree(g)
[1] 2 2 2
</code>```

Vertex 1 has links to both vertices 2 and 3. These count in the
in-degree and also count in the out-degree, so

Freeman = in + out = 2 + 2 = 4

The "Note" in the documentation states this.

Body:
I am doing some basic network analysis using networks from the R package "networkdata". To this end, I use the package "igraph" as well as "sna". However, I realised that the results of descriptive network statistics vary depending on the package I use. Most variation is not too grave but the average degree of my undirected graph halved as soon as I switched from "sna" to "igraph".

library(networkdata)
n_1 <- covert_28

library(igraph)
library(sna)

n_1_adjmat <- as_adjacency_matrix(n_1)
n_1_adjmat2 <- as.matrix(n_1_adjmat)

mean(sna::degree(n_1_adjmat2, cmode = "freeman")) # [1] 23.33333
mean(igraph::degree(n_1, mode = "all")) # [1] 11.66667

This doesn't happen in case of my directed graph. Here, I get the same results regardless of using "sna" or "igraph".

Is there any explanation for this phenomenon? And if so, is there anything I can do in order to prevent this from happening?

Thank you in advance!

Body:
This is explained in the documentation for sna::degree.

indegree of a vertex, v, corresponds to the cardinality
  of the vertex set N^+(v) = {i in V(G) : (i,v) in E(G)};
  outdegree corresponds to the cardinality of the vertex
  set N^-(v) = {i in V(G) : (v,i) in E(G)}; and total
  (or “Freeman”) degree corresponds to |N^+(v)| + |N^-(v)|.

(Note that, for simple graphs,

indegree=outdegree=total degree/2.)

A simpler example than yours makes it clear.

library(igraph)
library(sna)

g = make_ring(3)
plot(g)

AM = as.matrix(as_adjacency_matrix(g))
sna::degree(AM)
[1] 4 4 4

igraph::degree(g)
[1] 2 2 2

Vertex 1 has links to both vertices 2 and 3. These count in the
in-degree and also count in the out-degree, so

Freeman = in + out = 2 + 2 = 4

The "Note" in the documentation states this.


问:
Highlight rows depending on month

I have a table and a input where a person is supposed to enter a month 1-12, and when submitting, the rows from that month should turn red. How do I get the month out of these dates to determine which month they are from? I am unsure how to continue.

```<code> <input type="text" placeholder="Kuukausi" id="searchMonth"/>
            <button type="submit" id="searchbutton">
            <span>Etsi</span></button>
            </div>
            <table style="width: 100%" id="table">
            <tr>
                <th>Lintu</th>
                <th>Päivä</th>
            </tr>
            <tr>
                <td>Varis</td>
                <td>20.4.2018</td>
            </tr>
            <tr>
                <td>Huuhkaja</td>
                <td>14.9.2018</td>
            </tr>
            <tr>
                <td>Pääskynen</td>
                <td>24.4.2018</td>
            </tr>
            <tr>
                <td>Peippo</td>
                <td>30.3.2018</td>
            </tr>
</code>```

答:
Firstly, You can get month value by using split like currentText.split(".")[1]

Secondly, The input type of searchMonth should be number to allow number value only.

$("#searchbutton").on("click", function(){
   var searchMonth = $("#searchMonth").val();

   $("table tr").each(function(){
      var currentText = $(this).find(":eq(1)").text();
      var currentMonth = currentText.split(".")[1];

      var background_color = searchMonth == parseInt(currentMonth) ? "#00ff78" : "";
        $(this).css("background-color", background_color);
   });
});```
<pre class="snippet-code-html lang-html prettyprint-override"><code><script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/3.3.1/jquery.min.js"></script>
<input type="number" placeholder="Kuukausi" id="searchMonth" value = "4"/>
<button type="submit" id="searchbutton">
            <span>Etsi</span></button>
            </div>
            <table style="width: 100%" id="table">
            <tr>
                <th>Lintu</th>
                <th>Päivä</th>
            </tr>
            <tr>
                <td>Varis</td>
                <td>20.4.2018</td>
            </tr>
            <tr>
                <td>Huuhkaja</td>
                <td>14.9.2018</td>
            </tr>
            <tr>
                <td>Pääskynen</td>
                <td>24.4.2018</td>
            </tr>
            <tr>
                <td>Peippo</td>
                <td>30.3.2018</td>
            </tr>
            <tr>
                <td>Peippo</td>
                <td>20.02.2018</td>
            </tr></code>```

Title: Highlight rows depending on month

Body:
I have a table and a input where a person is supposed to enter a month 1-12, and when submitting, the rows from that month should turn red. How do I get the month out of these dates to determine which month they are from? I am unsure how to continue.

 <input type="text" placeholder="Kuukausi" id="searchMonth"/>
            <button type="submit" id="searchbutton">
            Etsi</button>

                Lintu
                Päivä

                Varis
                20.4.2018

                Huuhkaja
                14.9.2018

                Pääskynen
                24.4.2018

                Peippo
                30.3.2018

Body:
Firstly, You can get month value by using split like currentText.split(".")[1]

Secondly, The input type of searchMonth should be number to allow number value only.

$("#searchbutton").on("click", function(){
   var searchMonth = $("#searchMonth").val();

   $("table tr").each(function(){
      var currentText = $(this).find(":eq(1)").text();
      var currentMonth = currentText.split(".")[1];

      var background_color = searchMonth == parseInt(currentMonth) ? "#00ff78" : "";
        $(this).css("background-color", background_color);
   });
});

<script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/3.3.1/jquery.min.js"></script>
<input type="number" placeholder="Kuukausi" id="searchMonth" value = "4"/>
<button type="submit" id="searchbutton">
            <span>Etsi</span></button>
            </div>
            <table style="width: 100%" id="table">
            <tr>
                <th>Lintu</th>
                <th>Päivä</th>
            </tr>
            <tr>
                <td>Varis</td>
                <td>20.4.2018</td>
            </tr>
            <tr>
                <td>Huuhkaja</td>
                <td>14.9.2018</td>
            </tr>
            <tr>
                <td>Pääskynen</td>
                <td>24.4.2018</td>
            </tr>
            <tr>
                <td>Peippo</td>
                <td>30.3.2018</td>
            </tr>
            <tr>
                <td>Peippo</td>
                <td>20.02.2018</td>
            </tr>


问:
Regex get subfolder name from URL after domain extension
i'm new to regex, so I want to write a regExp to identify a special subfolder path directly after the tld top level domain extension.

http://sub.domain.tld/foo/bar

How to do it?

答:
Assuming your tld can be any of the usual - not literally tld, try

```<code>\.\w*\/([^\/\s]*)  // group 1 will have your desired folder
</code>```

regex101 demo with test cases

Body:
Assuming your tld can be any of the usual - not literally tld, try

\.\w*\/([^\/\s]*)  // group 1 will have your desired folder

regex101 demo with test cases


问:
C#: Populating combobox from stored procedure only displays first item
I am trying to populate a combobox at initialization by executing a stored procedure. When I run the procedure on my DB, there are 3 rows that come out:

But once I execute the procedure in my code:

```<code>public partial class AddInstitutionStart : Window
    {
        public AddInstitutionStart()
        {
            InitializeComponent();
            ComboBoxInstitutionCategory.Items.Add(FillCombo.fillInstCategory());
        }
    }

    public static class FillCombo
    {
        public static string fillInstCategory()
        {
            string connstr = ConfigurationManager.ConnectionStrings["connstr"].ConnectionString;
            SqlConnection Connect = null;
            SqlCommand Command;
            using (Connect = new SqlConnection(connstr))
            {
                Connect.Open();
                using (Command = new SqlCommand("[dbo].[spMetadataTb_GetCategory]", Connect))
                {
                    Command.CommandType = CommandType.StoredProcedure;
                    SqlDataReader dr = Command.ExecuteReader();
                    while (dr.Read())
                    {
                        string categoryInst = dr.GetString(0);
                        return categoryInst;
                    }
                }

                Command.Dispose();
                Connect.Close();
                return "-2";
            }
        }
    }
</code>```

Only the first row is displayed:

I have used the same logic to fill another combobox (not at initialization) and it works perfectly and displays multiple item. Any idea on what I might be doing wrong here?

答:
Tons of possible solutions for this, but please forget the static!

1st:

return a List from the other answer, use it like this:

```<code>public partial class AddInstitutionStart : Window
{
    public AddInstitutionStart()
    {
        InitializeComponent();
        foreach(string cat in new FillCombo().fillInstCategory())
            ComboBoxInstitutionCategory.Items.Add(cat);
    }
}
</code>```

or, if AddRange exists:

```<code>public partial class AddInstitutionStart : Window
{
    public AddInstitutionStart()
    {
        InitializeComponent();
        ComboBoxInstitutionCategory.Items.AddRange(new FillCombo().fillInstCategory());
    }
}
</code>```

This solution 1st fills the list in fillInstCategory(), then fills your combobox via foreach or addrange.

2nd:

make it to an enumerable function (use it same as above):

Only 2 change in your code, 1st:

```<code>public IEnumerable<string> GetInstCategory()
{
    ...
}
</code>```

then extend the return word with yield:

```<code>        while (dr.Read())
            yield return dr.GetString(0);
</code>```

this fills the combobox right after each data is read (not once, when all data is read as with solution 1), because yield causes jumping back and forth between your Items.Add and While(dr.Read()) code parts.

3rd:

inject a method (or lambda expression):

Only 2 change in your code, 1st:

```<code>public void GetInstCategory(Action<string> Injection)
{
    ...
}
</code>```

then:

```<code>        while (dr.Read())
            Injection(dr.GetString(0));
</code>```

Use it like:

```<code>public partial class AddInstitutionStart : Window
{
    public AddInstitutionStart()
    {
        InitializeComponent();
        new FillCombo().fillInstCategory(ComboBoxInstitutionCategory.Items.Add)) //note there is no "(" and ")" after Add, so the addition doesnt happen, just the function was passed to "fillInstCategory" to use it whenever necessary.
    }
}
</code>```

This fills the same way your combobox as solution 2: Reads 1 data, then Adds 1 to combobox, then reads again...

I hope you learn something from this.

Body:
I am trying to populate a combobox at initialization by executing a stored procedure. When I run the procedure on my DB, there are 3 rows that come out:

But once I execute the procedure in my code:

public partial class AddInstitutionStart : Window
    {
        public AddInstitutionStart()
        {
            InitializeComponent();
            ComboBoxInstitutionCategory.Items.Add(FillCombo.fillInstCategory());
        }
    }

    public static class FillCombo
    {
        public static string fillInstCategory()
        {
            string connstr = ConfigurationManager.ConnectionStrings["connstr"].ConnectionString;
            SqlConnection Connect = null;
            SqlCommand Command;
            using (Connect = new SqlConnection(connstr))
            {
                Connect.Open();
                using (Command = new SqlCommand("[dbo].[spMetadataTb_GetCategory]", Connect))
                {
                    Command.CommandType = CommandType.StoredProcedure;
                    SqlDataReader dr = Command.ExecuteReader();
                    while (dr.Read())
                    {
                        string categoryInst = dr.GetString(0);
                        return categoryInst;
                    }
                }

                Command.Dispose();
                Connect.Close();
                return "-2";
            }
        }
    }

Only the first row is displayed:

I have used the same logic to fill another combobox (not at initialization) and it works perfectly and displays multiple item. Any idea on what I might be doing wrong here?

Body:
Tons of possible solutions for this, but please forget the static!

1st:

return a List from the other answer, use it like this:

public partial class AddInstitutionStart : Window
{
    public AddInstitutionStart()
    {
        InitializeComponent();
        foreach(string cat in new FillCombo().fillInstCategory())
            ComboBoxInstitutionCategory.Items.Add(cat);
    }
}

or, if AddRange exists:

public partial class AddInstitutionStart : Window
{
    public AddInstitutionStart()
    {
        InitializeComponent();
        ComboBoxInstitutionCategory.Items.AddRange(new FillCombo().fillInstCategory());
    }
}

This solution 1st fills the list in fillInstCategory(), then fills your combobox via foreach or addrange.

2nd:

make it to an enumerable function (use it same as above):

Only 2 change in your code, 1st:

public IEnumerable<string> GetInstCategory()
{
    ...
}

then extend the return word with yield:

        while (dr.Read())
            yield return dr.GetString(0);

this fills the combobox right after each data is read (not once, when all data is read as with solution 1), because yield causes jumping back and forth between your Items.Add and While(dr.Read()) code parts.

3rd:

inject a method (or lambda expression):

Only 2 change in your code, 1st:

public void GetInstCategory(Action<string> Injection)
{
    ...
}

then:

        while (dr.Read())
            Injection(dr.GetString(0));

Use it like:

public partial class AddInstitutionStart : Window
{
    public AddInstitutionStart()
    {
        InitializeComponent();
        new FillCombo().fillInstCategory(ComboBoxInstitutionCategory.Items.Add)) //note there is no "(" and ")" after Add, so the addition doesnt happen, just the function was passed to "fillInstCategory" to use it whenever necessary.
    }
}

This fills the same way your combobox as solution 2: Reads 1 data, then Adds 1 to combobox, then reads again...

I hope you learn something from this.