GitHub

Source localization

This example illustrates source localization via multi-dimensional scaling using the Julia language.

This page comes from a single Julia file: source-local.jl.

You can access the source code for such Julia documentation using the 'Edit on GitHub' link in the top right. You can view the corresponding notebook in nbviewer here: source-local.ipynb, or open it in binder here: source-local.ipynb.

Setup

Add the Julia packages used in this demo. Change false to true in the following code block if you are using any of the following packages for the first time.

if false
    import Pkg
    Pkg.add([
        "InteractiveUtils"
        "LinearAlgebra"
        "MIRTjim"
        "Plots"
        "Random"
        "Statistics"
    ])
end

Tell Julia to use the following packages. Run Pkg.add() in the preceding code block first, if needed.

using InteractiveUtils: versioninfo
using LaTeXStrings
using LinearAlgebra: svd, norm, Diagonal
using MIRTjim: jim, prompt
using Plots: default, scatter, savefig
using Random: seed!
using Statistics: mean
default(); default(label="", markerstrokecolor=:auto, widen=true,
    guidefontsize=14, tickfontsize=12, legendfontsize=14)

The following line is helpful when running this jl-file as a script; this way it will prompt user to hit a key after each image is displayed.

isinteractive() && prompt(:prompt);
isinteractive() && jim(:prompt);

Generate data

Coordinates - can you identify the pattern? probably not...

C0 = [
3 2.5 2 2 2 2 2 1 0 0 1 1 1 1 0 0 1 2 2.5 3  4 4.5 5 6 7 7 6 6 6 6 7 7 6 5 5 5 5 5 4.5 4;
2 3.0 4 3 2 1 0 0 0 1 1 2 3 4 4 5 5 5 4.0 3  3 4.0 5 5 5 4 4 3 2 1 1 0 0 0 1 2 3 4 3.0 2
];

Interpolate the locations to make the final set more familiar

C2 = (C0[:,2:end] + C0[:,1:end-1])/2
Ce = (C0[:,1] + C0[:,end])/2
C3 = [C2 Ce]
C4 = [C0; C3]
C = reshape(C4, 2, :);

Compute J × J distance array and display it

J = size(C,2) # number of points
D = [norm(C[:,j] - C[:,i]) for i in 1:J, j in 1:J] # "comprehension" in julia!
pd = jim(D, L"D", color=:cividis, xlabel=L"j", ylabel=L"i")

# savefig(pd, "06_source_local1_d.pdf")
Example block output

MDS algorithm

Compute Gram matrix by de-meaning squared distance matrix

S = D.^2 # squared distances
G = S .- mean(S, dims=1) # trick: use "broadcasting" feature of julia
G = G .- mean(G, dims=2) # now we have de-meaned the columns and rows of S
G = -1/2 * G;

We still cannot determine visually the point locations:

pg = jim(G, L"G", color=:cividis, xlabel=L"j", ylabel=L"i")

# savefig(pg, "06_source_local1_g.pdf")
Example block output

Examine singular values

(_, σ, V) = svd(G) # svd returns singular values in descending order
ps = scatter(σ, label="singular values (noiseless case)",
    xlabel=L"k", ylabel=L"σ_k") # two nonzero (d=2)
Example block output
prompt()

# savefig(ps, "06_source_local1_eig.pdf")

Estimate the source locations using rank=2

Ch = Diagonal(sqrt.(σ[1:2])) * V[:,1:2]' # here is the key step
2×80 Matrix{Float64}:
 -0.5  -0.75  -1.0  -1.25  -1.5  -1.5  …   1.25   1.0  0.75  0.5  6.99074e-17
  0.6   0.1   -0.4  -0.9   -1.4  -0.9     -0.9   -0.4  0.1   0.6  0.6

Plot estimated source locations

pc = scatter(Ch[1,:], -Ch[2,:], xtick=-4:4, ytick=-3:3, aspect_ratio=1,
 title="Location estimates (noiseless case)")
Example block output
prompt()

Noisy case

seed!(0)
Dn = D + 0.3 * randn(size(D))
Sn = Dn.^2
Gn = Sn .- mean(Sn, dims=1)
Gn = Gn .- mean(Gn, dims=2) # de-meaned
Gn = -1/2 * Gn
pgn = jim(Gn, "G noisy", color=:cividis)
Example block output

Singular values

(_, sn, Vn) = svd(Gn)
psn = scatter(sn, label="singular values (noisy case)") # σ₂ ≫ σ₃
Example block output
prompt()

Plot estimated source locations from noisy distance measurements

Cn = Diagonal(sqrt.(sn[1:2])) * Vn[:,1:2]'
pcn = scatter(Cn[1,:], -Cn[2,:], xtick=-4:4, ytick=-3:3, aspect_ratio=1,
 title="Location estimates (noisy case)")
Example block output
prompt()

Constant bias case

seed!(0)
Db = D .+ 0.3 * maximum(D) # fairly large bias
Sb = Db.^2
Gb = Sb .- mean(Sb, dims=1)
Gb = Gb .- mean(Gb, dims=2) # de-meaned
Gb = -1/2 * Gb
pgb = jim(Gb, "G biased", color=:cividis)
Example block output

Singular values

(_, sb, Vb) = svd(Gb)
psb = scatter(sb, label="singular values (biased case)") # σ₂ ≫ σ₃
Example block output
prompt()

Plot estimated source locations from biased distance measurements

Cb = Diagonal(sqrt.(sb[1:2])) * Vb[:,1:2]'
pcb = scatter(Cb[1,:], -Cb[2,:], xtick=-4:4, ytick=-3:3, aspect_ratio=1,
 title="Location estimates (biased case)")
Example block output
prompt()

Equilateral triangle example

G = [-2 1 1; 1 -2 1; 1 1 -2] / (-6.)
(~, σ, V) = svd(G)
Ch = Diagonal(sqrt.(σ[1:2])) * V[:,1:2]'
scatter(Ch[1,:], Ch[2,:], aspect_ratio=1, title="Location estimates")
Example block output
prompt()

Reproducibility

This page was generated with the following version of Julia:

using InteractiveUtils: versioninfo
io = IOBuffer(); versioninfo(io); split(String(take!(io)), '\n')
12-element Vector{SubString{String}}:
 "Julia Version 1.10.1"
 "Commit 7790d6f0641 (2024-02-13 20:41 UTC)"
 "Build Info:"
 "  Official https://julialang.org/ release"
 "Platform Info:"
 "  OS: Linux (x86_64-linux-gnu)"
 "  CPU: 4 × AMD EPYC 7763 64-Core Processor"
 "  WORD_SIZE: 64"
 "  LIBM: libopenlibm"
 "  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)"
 "Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)"
 ""

And with the following package versions

import Pkg; Pkg.status()
Status `~/work/book-la-demo/book-la-demo/docs/Project.toml`
  [6e4b80f9] BenchmarkTools v1.5.0
  [aaaa29a8] Clustering v0.15.7
  [35d6a980] ColorSchemes v3.24.0
  [3da002f7] ColorTypes v0.11.4
⌅ [c3611d14] ColorVectorSpace v0.9.10
  [717857b8] DSP v0.7.9
  [72c85766] Demos v0.1.0 `~/work/book-la-demo/book-la-demo`
  [e30172f5] Documenter v1.2.1
  [4f61f5a4] FFTViews v0.3.2
  [7a1cc6ca] FFTW v1.8.0
  [587475ba] Flux v0.14.12
⌅ [a09fc81d] ImageCore v0.9.4
  [71a99df6] ImagePhantoms v0.7.2
  [b964fa9f] LaTeXStrings v1.3.1
  [7031d0ef] LazyGrids v0.5.0
  [599c1a8e] LinearMapsAA v0.11.0
  [98b081ad] Literate v2.16.1
  [7035ae7a] MIRT v0.17.0
  [170b2178] MIRTjim v0.23.0
  [eb30cadb] MLDatasets v0.7.14
  [efe261a4] NFFT v0.13.3
  [6ef6ca0d] NMF v1.0.2
  [15e1cf62] NPZ v0.4.3
  [0b1bfda6] OneHotArrays v0.2.5
  [429524aa] Optim v1.9.2
  [91a5bcdd] Plots v1.40.1
  [f27b6e38] Polynomials v4.0.6
  [2913bbd2] StatsBase v0.34.2
  [d6d074c3] VideoIO v1.0.9
  [b77e0a4c] InteractiveUtils
  [37e2e46d] LinearAlgebra
  [44cfe95a] Pkg v1.10.0
  [9a3f8284] Random
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

This page was generated using Literate.jl.