[Math] Automatic calculation of the intersection of discrete curves

interpolationmathematicaMATLABpolynomials

first of all, let me apologize for a poor math-english translation, I'll try my very best.

I have the following situation:
I have over 16.000 data files which I generated from a biometric verification system.
Those files contain Error Rates for False-negative and False-positive at a specific threshold.
So every data file becomes one Plot to detect an equal error rate (intersection of both rates). Special is in this case, that it may also contain more than one intersection, which is not typical but possible in this case.

I used to use gnuplot before to generate plots, but now the amount is a new level.

I can't watch every plot by myself which would take days. And it's for sure that lots of the files may not contain any intersection in the result.

Therefore I wonder if there is any possibility to calculate for each file the intersection(s) with MATLAB, Mathematica or anything else?

And after having calculated maybe generate a plot (preferred svg format but not important in this step) for that file so that I end up with only "useful" data files and plots with a fixed threshold.

I hope I could explain what my problem is and seeing forward for your ideas.

(Note: Both examples do not belong together, just choosed quickly and randomly from the files)

Data file example:

> #k    NFA NIA NFR NEA FAR FRR 
> #Alpha Vector 0 0 -1 0 0 0 -1 0 -1 0 0 -1 
> 0.4   280.0   280.0   0.0 40.0    100 0
> 0.41  280.0   280.0   0.0 40.0    100 0
> 0.42  280.0   280.0   0.0 40.0    100 0
> 0.43  280.0   280.0   0.0 40.0    100 0
> 0.44  280.0   280.0   0.0 40.0    100 0
> 0.45  280.0   280.0   0.0 40.0    100 0
> 0.46  280.0   280.0   0.0 40.0    100 0
> 
> # k is the threshold, FAR and FRR are 2 curves. k is on x-axis, value of FAR/FRR on y.

Example Plot:
EER with gnuplot

addon information.

I have no experience with matlab and just worked with wolframalpha mathematica during calculus at university. I study computer sciences in Germany, so ideas for scripting and co are also welcome and probably welcomed 🙂

Best Answer

Let $k_1,\ldots, k_N$ be your points on the $x$-axis. Denote the two data sets you want to compare by $f_i := f(k_i)$ and $g_i := g(k_i)$.

First, it can happen that $f_i = g_i$ for some $i$. In this case you already found an intersection.

Second, you can have that $f_i > g_i$ and $f_{i+1} < g_{i+1}$ (or vice versa). This means that there is an intersection in the interval $(k_i, k_{i+1})$. Does this information suffice?

If not: Where exactly the point of intersection is located, depends on how you interpolate the data. The simplest type of interpolation is probably a linear one (as in your plot), in which case you only have to find the intersection of two lines. Can you do that?

Related Solutions

[Math] Finding points along a plot in Matlab

Maybe this code can help you:

function mse


x = linspace(0,10,100);
y = sin(x)+x-(x/5).^2;
hold off

plot(x,y)

yMax = max(y);
xMax = x(find(y==yMax));

hold on
plot(xMax,yMax,'*g')



y90 = yMax*0.9;
plot(x,y90,'-k')


for i=2:length(x)
    if y(i-1)<=y90 && y(i)>=y90 || y(i-1)>=y90 && y(i)<=y90
        plot(x(i),y(i),'og')
        plot(x(i-1),y(i-1),'og')

    end
end


y50 = yMax*0.5;
plot(x,y50,'-k')
for i=2:length(x)
    if y(i-1)<=y50 && y(i)>=y50 || y(i-1)>=y50 && y(i)<=y50
        plot(x(i),y(i),'og')
        plot(x(i-1),y(i-1),'og')

    end
end
end

Here $x,y$ are my data. Finding the maximum is easy, but for the 90% or 50% values I iterated through the data. As we do not have continous data we check, whether the 90% or 50% values lies between to points. The result is shown below. If you only want to highlight a single point instead of the two you could perform some (linear) interpolation. enter image description here

EDIT: This would highlight all points that match the X% criteria. If you only want to highlight the closest points to the maximum, you could change the loop to something like

 for i=find(y==ymax):length(x)

Which finds the point on the right and for the left side

 for i=find(y==ymax):-1:length(x)

To mark only one point you could use a break command.

EDIT2: The following code performs the same thing as the for loops, but it is shorter.

y90 = yMax*0.9;
xId= find(diff(sign(y-y90))~=0);
plot(x(xId),y(xId),'*k')
plot(x(xId+1),y(xId+1),'*k')

EDIT3: For your date you need to interpolate between the two values. Replace the loop with something like

for i=2:length(x)
    if y(i-1)<=y90 && y(i)>=y90 || y(i-1)>=y90 && y(i)<=y90
        plot(x(i),y(i),'og')
        plot(x(i-1),y(i-1),'og')

        xl = x(i-1);
        yl = y(i-1);
        xr = x(i);
        yr = y(i);

        xMid  = (y90-yl)*(xr-xl)/(yr-yl)+xl;
        plot(xMid,y90,'or');

    end
end

[Math] Tricky discrete math problems

A solution for $N=7$ with the smallest value for $\omega_N$ is this:

$$ W = \begin{bmatrix} 25 & 26 & 28 & 29 & 34 & 36 & 37 \end{bmatrix} $$

$$M = \begin{bmatrix} 0 & 24 & 22 & 21 & 16 & 14 & 13 \\ 27 & 0 & 24 & 23 & 18 & 16 & 15 \\ 31 & 30 & 0 & 27 & 22 & 20 & 19 \\ 33 & 32 & 30 & 0 & 24 & 22 & 21 \\ 43 & 42 & 40 & 39 & 0 & 32 & 31 \\ 47 & 46 & 44 & 43 & 38 & 0 & 35 \\ 49 & 48 & 46 & 45 & 40 & 38 & 0 \end{bmatrix}$$

This solution is generated by this z3py script:

""" Find 0 < w1 < w2 < ... < wN < 1.5 w_1 such that all values of
    2wi - wj for i != j are distinct from all wk.
"""
from __future__ import print_function
from z3 import *

N = 7
omega = [Int('omega_%s' % (i+1)) for i in range(N)]

s = Solver()

s.add(omega[0] > 0)
for i in range(1,N):
    s.add(omega[i] > omega[i-1])
s.add(2*omega[-1] < 3*omega[0])

# Impose upper bound on omega[-1] to find "smallest" solution.
s.add(omega[-1] < 38)

for i in range(N):
    for j in range(i+1,N):
        for w in omega:
            s.add(2*omega[i] - omega[j] != w)
            s.add(2*omega[j] - omega[i] != w)

result = s.check()
if result == sat:
    mdl = s.model()
    print('W =', ' '.join([str(mdl.evaluate(w).as_long()) for w in omega]))

    mat = [[0] * N for _ in range(N)]
    for i in range(N):
        for j in range(i):
            orow = mdl.evaluate(omega[i]).as_long()
            ocol = mdl.evaluate(omega[j]).as_long()
            mat[i][j] = 2 * orow - ocol
            mat[j][i] = 2 * ocol - orow

    print('\n'.join([' '.join(['{:3}'.format(item) for item in row])
                    for row in mat]))
else:
    print(result)

Data file example:

Best Answer

Related Solutions

[Math] Finding points along a plot in Matlab

[Math] Tricky discrete math problems

Related Question