In [1]:
using Laplacians
In [2]:
n = 100000
x1 = []
x2 = []
for i in 1:10
a = wtedChimera(n)
t1 = randishPrim(a)
t2 = randishKruskal(a)
s1 = sum(compStretches(t1,a))/nnz(a)
s2 = sum(compStretches(t2,a))/nnz(a)
push!(x1,s1)
push!(x2,s2)
end
In [5]:
[x1 x2]
Out[5]:
In [4]:
include("../src/KMPSolver.jl")
include("../devel/kmpTest.jl")
Out[4]:
In [5]:
T = 6.0
p = T^(-2)
pList = [KMPparams(p,T,1/4,600,:rand)]
T = 4.0
p = T^(-2)
push!(pList,KMPparams(p,T,1/4,600,:rand))
push!(pList,KMPparams(p,T,1/4,600,:akpw))
out = manyRuns(10000,5,pList)
Out[5]:
In [38]:
a = chimera(10000,4)
la = lap(a)
fsub = lapWrapSolver(augTreeSolver,lap(a),tol=0.01,maxits=10000)
In [45]:
la = lap(a)
c = components(la[1:9999,1:9999])
maximum(c)
Out[45]:
In [17]:
n = 100000
a = wtedChimera(n)
b = randn(n)
b = b - mean(b)
la = lap(a)
b2 = la*b
f = KMPLapSolver(a,verbose=true)
@time x = f(b);
@time x2 = f(b2);
In [12]:
out = []
push!(out,[1 2 3])
push!(out,[4 5 6])
out
Out[12]:
In [18]:
f = lapWrapSolver(augTreeSolver,la,tol=0.1,maxits=1000)
#x = f(b)
#norm(la*x-b)/norm(b)
Out[18]:
In [ ]:
n = 1000000
a = chimera(n,4)
la = lap(a)
b = randn(n);
b = b - mean(b);
#t = akpw(a)
#heavy = makeHeavy(a,t)
#lheavy = lap(heavy);
A test of the solver code, with some param settings thrown in
In [ ]:
KMP_LOGGING = true
KMP_SAVEMATS = true
t = akpw(a)
T = 6.0
p = (16/9)*T^(-2)
params=KMPparams(p,T,1,500)
@time f = KMPLapSolver(a,params=params,tol=0.1,maxits=1000)
@time x = f(b)
@show norm(la*x-b)/norm(b);
In [ ]:
Profile.clear()
@profile f(b);
ProfileView.view()
In [ ]:
(353546 - 3*(99999))/2
In [ ]:
KMP_MATS
In [ ]:
KMP_FS
In [ ]:
la2 = KMP_MATS[1];
f2 = KMP_FS[2];
f1 = KMP_FS[2];
n2 = size(la2,1)
In [ ]:
b = randn(n2);
b = b - mean(b);
x = f1(b);
norm(la2*x-b)/norm(b)
A comparison with our old solver
In [ ]:
@time fa = lapWrapSolver(augTreeSolver,la,tol=0.1,maxits=1000)
@time xa = fa(b)
@show norm(la*xa-b)/norm(b);
In [ ]:
Profile.clear()
@profile fa(b);
ProfileView.view()
In [ ]:
To understand what's going on, let's look at how quickly it solves the precon system.
In [ ]:
KMP_LOGGING = true
# t = akpw(a)
T = 6.0
p = 2*T^(-2)
ord = Laplacians.dfsOrder(t)
aord = a[ord,ord]
tord = t[ord,ord]
params=KMPparams(p,T,2,500)
heavy = makeHeavy(aord, t=tord, params=params)
lh = lap(heavy)
@time fsub = KMPLapPrecon(aord, tord, params)
@time xh = fsub(b)
@show norm(lh*xh-b)/norm(b);
In [ ]:
xh = pcg(lh, b, fsub, tol=1e-2, maxits=60)
norm(lh*xh-b)/norm(b)
In [ ]:
nnz(triu(a,1))
In [ ]:
Profile.clear()
In [ ]:
@profile fsub(b);
In [ ]:
using ProfileView
ProfileView.view()
In [ ]:
#Profile.clear()
n = 1000000
a = chimera(n,1)
b = randn(n)
@profile a*b;
ProfileView.view()
In [ ]:
function ms1(x)
mn = mean(x)
for i in 1:length(x)
x[i] = x[i] - mn
end
end
In [ ]:
n = 10000900
x = randn(n)
@time x = x - mean(x);
y = randn(n)
@time ms1(y)
mean(y)
In [ ]:
function manyRuns(n,numruns,pList)
tot = zeros(length(pList))
for it in 1:numruns
a = chimera(n,n+it)
b = randn(n)
b = b - mean(b)
for i in 1:length(pList)
fsub = KMPLapPrecon(a, akpw(a), pList[i])
tic()
xh = fsub(b)
tot[i] += toq()
end
end
for i in 1:length(pList)
println(string(pList[i]), " : ", tot[i])
end
return tot
end
function manyRunsW(n,numruns,pList)
tot = zeros(length(pList))
for it in 1:numruns
a = wtedChimera(n,n+it)
b = randn(n)
b = b - mean(b)
for i in 1:length(pList)
fsub = KMPLapPrecon(a, akpw(a), pList[i])
tic()
xh = fsub(b)
tot[i] += toq()
end
end
for i in 1:length(pList)
println(string(pList[i]), " : ", tot[i])
end
return tot
end
The following gave a best n0=400, at least on a few small experiments. Another puts it at 800.
In [ ]:
T = 6.0
p = (16/9)*T^(-2)
pList = [KMPparams(p,T,2,100)]
push!(pList,KMPparams(p,T,2,200))
push!(pList,KMPparams(p,T,2,400))
push!(pList,KMPparams(p,T,2,800))
push!(pList,KMPparams(p,T,2,1600))
zeros(length(pList))
In [ ]:
KMP_LOGGING = true
totx = manyRuns(1000,2,pList);
In [ ]:
tot2 = manyRunsW(40000,50,pList)
In [ ]:
tot
In [ ]:
tot2
The following experiment suggests T = 6.0
In [ ]:
T = 4.0
pList = [KMPparams((16/9)*T^(-2),T,2,600)]
T = 6.0
push!(pList,KMPparams((16/9)*T^(-2),T,2,600))
T = 8.0
push!(pList,KMPparams((16/9)*T^(-2),T,2,600))
tot3 = manyRuns(30000,70,pList)
In [ ]:
tot3
In [ ]:
print(string(pList[1]))
In [ ]:
In [ ]:
marked = 0 + (rand(n) .< 0.6)
b = randn(n)
b = b - mean(b)
lt = lap(t)
lt[1,1] = lt[1,1] + 1
x = recTest12x(t, b, marked)
norm(lt*x-b);
In [ ]:
@time F = factorize(lt)
@time y = F \ b
norm(y - x)
In [ ]:
elims1, elims2 = elimDeg12(t, marked)
In [ ]:
(rand(5) .< 0.6) + 0
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
elims1
In [ ]:
elims2
In [ ]:
ind = find(marked.<2)
In [ ]:
ts = t[ind,ind];
lt = lap(t)
lt[1,1] = lt[1,1] + 1
lts = lap(ts)
lts[1,1] = lts[1,1] + 1
b = randn(n);
In [ ]:
cumb = copy(b)
for i in 1:length(elims1)
cumb[elims1[i].parent] += cumb[elims1[i].nodeid]
end
cumbs = cumb[ind]
sum(cumbs)
In [ ]:
xs = lts \ cumbs
norm(lts*xs-cumbs)
In [ ]:
x = zeros(n)
x[ind] = xs
for i in length(elims1):-1:1
node = elims1[i].nodeid
x[node] = x[elims1[i].parent] + cumb[node]/elims1[i].wtDeg
end
norm(lt*x-b)
In [ ]:
y = la \ b;
norm(la*y-b)
[y[ind] x[ind] xs]
In [ ]:
sum(las)
In [ ]:
In [ ]:
In [ ]:
function recTest(t,b, marked)
n = length(b)
elims1 = elimDeg1(t, marked)
cumb = copy(b)
for i in 1:length(elims1)
cumb[elims1[i].parent] += cumb[elims1[i].nodeid]
end
ind = find(marked.<2)
ts = t[ind,ind];
lts = lap(ts)
cumbs = cumb[ind]
xs = pinv(full(lts))*cumbs
x = zeros(n)
x[ind] = xs
for i in length(elims1):-1:1
node = elims1[i].nodeid
x[node] = x[elims1[i].parent] + cumb[node]/elims1[i].wtDeg
end
return x
end
In [ ]:
mkd = rand(0:1,n)
#elims1 = elimDeg1(t, mkd)
#ind = find(mkd .< 2)
In [ ]:
cumb = copy(b)
for i in 1:length(elims1)
cumb[elims1[i].parent] += cumb[elims1[i].nodeid]
end
cumbs = cumb[ind]
sum(cumbs)
In [ ]:
sum(cumbs)
In [ ]:
b = randn(n)
b = b - mean(b)
lt = lap(t)
x = pinv(full(lt)) * b
norm(lt*x-b)
In [ ]:
x = recTest(t, b, mkd)
norm(lt*x-b)
In [ ]:
sum(mkd.==2)
In [ ]:
norm(b)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
function makeElimList(t)
tr = matToTree(t)
n = size(tr.children,1)
elims = Array{elimTreeNode}(0)
for vi in n:-1:2
v = tr.children[vi]
push!(elims,elimTreeNode(v,tr.parent[v],tr.weights[vi]))
end
return elims
end
In [ ]:
function solveTree(elims,b)
cumb = copy(b)
n = size(b,1)
for i in 1:(n-1)
cumb[elims[i].parent] += cumb[elims[i].nodeid]
end
x = zeros(Float64,n)
for i in (n-1):-1:1
node = elims[i].nodeid
x[node] = x[elims[i].parent] + cumb[node]/elims[i].wtDeg
end
return x
end
In [ ]:
@time el = makeElimList(t);
In [ ]:
lt = lap(t)
lt[end,end] += 1
@time F = factorize(lt)
In [ ]:
b = randn(size(t,1));
b = b - mean(b);
@time x = solveTree(el,b);
x = x - mean(x);
norm(lap(t)*x-b)
In [ ]:
@time y = F \ b;
norm(lt*y-b)
In [ ]:
In [ ]:
In [ ]:
elims = Array{elimTreeNode}(0)
In [ ]:
for vi in n:-1:2
v = tr.children[vi]
push!(elims,elimTreeNode(v,tr.parent[v],tr.weights[vi]))
end
In [ ]:
elims
In [ ]:
b = randn(n)
b = b - mean(b)
In [ ]:
cumb = copy(b)
for i in 1:(n-1)
cumb[elims[i].parent] += cumb[elims[i].nodeid]
end
cumb
In [ ]:
x = zeros(Float64,n)
for i in (n-1):-1:1
node = elims[i].nodeid
x[node] = x[elims[i].parent] + cumb[node]/elims[i].wtDeg
end
x
In [ ]:
xs = pinv(full(lt))*b
In [ ]:
x-xs
In [ ]:
x = x - mean(x)
In [ ]:
lt*x - b
In [ ]:
b
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
tr.weights
In [ ]:
In [ ]:
In [ ]:
In [ ]:
st = compStretches(t,a);
aveStretch = sum(st)/nnz(a)
In [ ]:
rest = a-t;
fac = aveStretch*log(size(a,1))*2
heavy = rest+fac*t;
In [ ]:
fac
Sample with probability prop to stretch, but fix the number of edges we get to n/5
In [ ]:
n = size(a,1);
In [ ]:
strest = compStretches(t,rest)
probs = triu(strest);
probs = probs * (n/5)/ sum(probs);
In [ ]:
(pi,pj,pv) = findnz(probs)
In [ ]:
select = rand(size(pv)) .< pv;
sum(select)
In [ ]:
(ai,aj,av) = findnz(triu(rest))
In [ ]:
samp1 = sparse(ai,aj,av.*select./pv,n,n)
samp1 = samp1 + samp1';
In [ ]:
st1 = compStretches(t*fac, samp1);
In [ ]:
(s1i,s1j,s1v) = findnz(triu(samp1))
In [ ]:
maximum(x)
In [ ]:
a = chimera(100000,103+1111);
testGraph(a);
In [ ]:
include("../devel/kmpTest.jl")
In [ ]:
testAtSize(10000)
In [ ]:
x = testKMP(a,frac1 = 1/5, frac2 = 1/10)
In [ ]:
(pi,pj,pv) = findnz(probs)
select = rand(size(pv)) .< pv;
(ai,aj,av) = findnz(triu(rest))
samp1 = sparse(ai,aj,av.*select./pv,n,n)
# samp1 = samp1 + samp1';
# st1 = compStretches(t*fac, samp1);
# (s1i,s1j,s1v) = findnz(triu(samp1))
In [ ]:
ringGraph(1)
In [ ]:
a = wtedChimera(1000000)
@time t = akpw(a);
In [ ]:
n = size(t,1)
nnz(a)
In [ ]:
@time ord = Laplacians.dfsOrder(t);
In [ ]:
@time aord = symperm(a,ord);
@time aord = aord + aord';
In [ ]:
pe = randperm(n)
@time ape = symperm(a,pe)
@time ape = ape + ape';
In [ ]:
x = randn(n)
@time b1 = ape*x;
@time b2 = aord*x;
In [ ]:
@show nnz(ape)
@show nnz(aord)
In [ ]:
tord = t[ord,ord];
In [ ]:
@time st1 = compStretches(tord,aord);
In [ ]:
@time st2 = compStretchesDFS(tord,aord);
In [ ]:
sum(abs(st1-st2))
In [ ]:
@time tr, d1 = matToTreeDepth(tord);
In [ ]:
@time d2 = treeDepthDFS(tord);
In [ ]:
norm(d1-d2)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
function dfsOrder2(t::SparseMatrixCSC)
n = size(t,1)
seen = zeros(Bool,n)
ord = Array{Int64}(0)
stk = Array{Int64}(0)
push!(stk,1)
while ~isempty(stk)
u = pop!(stk)
push!(ord,u)
seen[u] = true
for vi in t.colptr[u]:(t.colptr[u+1]-1)
v = t.rowval[vi]
if !seen[v]
push!(stk,v)
end
end
end
return ord
end
In [ ]:
ord = dfsOrder(tr);
In [ ]:
stack = Array{Int64}(0)
In [ ]:
~isempty(stack)
In [ ]:
In [ ]:
In [ ]:
function dfsOrder2(tr)
ord = Array{Int64}(0)
stack = Array{Int64}(0)
push!(stack,tr.root)
while ~isempty(stack)
u = pop!(stack)
push!(ord,u)
for vi in tr.kidsPtr[u]:(tr.numKids[u] + tr.kidsPtr[u] - 1)
v = tr.children[vi]
push!(stack,v)
end
end
return ord
end
In [ ]:
ord2 = dfsOrder2(tr);
In [ ]:
norm(ord - ord2[size(tr.children,1):-1:1])
In [ ]:
# n = 100000
a = grid2(1000)
t = akpw(a)
aves = sum(compStretches(t,a))/nnz(a)
In [ ]:
avc = zeros(10)
for i in 1:10
ap = randperm(a)
tp = akpw(ap)
avc[i] = sum(compStretches(tp,ap))/nnz(a)
end
mean(avc)
In [ ]:
avc
In [ ]:
?akpw!
In [ ]:
a = grid2(700)
@time t = akpw(a,exponentialX=false)
@show sum(compStretches(t,a))/nnz(a)
@time t = akpw(a,exponentialX=true)
@show sum(compStretches(t,a))/nnz(a);
In [ ]:
nits = 100
df = Array{Float64}(0)
dt = Array{Float64}(0)
timef = Array{Float64}(0)
timet = Array{Float64}(0)
for i in 1:nits
a = wtedChimera(500000)
tic();
t = akpw(a,exponentialX=false)
push!(timef,toc())
push!(df,sum(compStretches(t,a))/nnz(a))
tic();
t = akpw(a,exponentialX=true)
push!(timet,toc())
push!(dt,sum(compStretches(t,a))/nnz(a))
a = unweight(a)
tic();
t = akpw(a,exponentialX=false)
push!(timef,toc())
push!(df,sum(compStretches(t,a))/nnz(a))
tic();
t = akpw(a,exponentialX=true)
push!(timet,toc())
push!(dt,sum(compStretches(t,a))/nnz(a))
end
In [ ]:
x = collect(1:length(dt));
plot(x,sort(timet),"rx",x,sort(timef),"bo")
In [ ]:
plot(x,sort(timet./timef),"rx",x,sort(timef./timet),"bo")
In [ ]:
plot(x,sort(dt./df),"rx",x,sort(df./dt),"bo")
In [ ]:
mean(df)
In [ ]:
mean(dt)
In [ ]:
plot(x,sort(df),"rx",x,sort(dt),"bo")
In [ ]:
In [ ]:
In [ ]:
In [ ]:
mean(timet)
In [ ]:
mean(df)
In [ ]:
mean(dt)
In [ ]:
function testSampler(a; t=akpw(a), frac1=1/5)
n = size(a,1);
rest = a-t;
st = compStretches(t,rest);
aveStretch = sum(st)/nnz(rest)
@show aveStretch
fac = aveStretch*log(size(a,1))*3
heavy = rest+fac*t;
(ai,aj,av) = findnz(triu(rest))
(si,sj,sv) = findnz(triu(st))
sv = sv ./ fac
restijvs = restIJVS(ai,aj,av,sv)
sampijvs = stretchSample2(restijvs,frac1)
samp1 = sparse(sampijvs.i,sampijvs.j,sampijvs.v,n,n)
samp1 = samp1 + samp1';
st1 = compStretches(t*fac, samp1);
(s1i,s1j,s1v) = findnz(triu(samp1))
add = speye(n)/10^6;
e = eigs(lap(heavy)+add,lap(samp1+t*fac)+add, nev=4);
return maximum(e[1])
end
In [ ]:
a = grid2(50)
t = akpw(a)
frac1 = 2/3
n = size(a,1);
rest = a-t;
st = compStretches(t,rest);
aveStretch = sum(st)/nnz(rest)
@show aveStretch
fac = aveStretch*log(size(a,1))
heavy = rest+fac*t;
(ai,aj,av) = findnz(triu(rest))
(si,sj,sv) = findnz(triu(st))
sv = sv ./ fac
restijvs = restIJVS(ai,aj,av,sv)
sampijvs = stretchSample(restijvs,frac1)
In [ ]:
include("../src/KMPSolver.jl")
include("../devel/kmpTest.jl")
In [ ]:
n = 1000000
a = chimera(n)
la = lap(a)
b = randn(n);
b = b - mean(b);
t = akpw(a)
heavy = makeHeavy(a,t)
lheavy = lap(heavy);
In [ ]:
In [ ]:
@profile f = KMPLapPrecon(a,t)
#@time xh = f(b)
#@show norm(lheavy*xh-b);
In [ ]:
Profile.clear()
In [ ]:
@time fa = lapWrapSolver(augTreeSolver,lheavy,tol=1e-2,maxits=1000)
@time xa = fa(b)
@show norm(lheavy*xa-b);
In [ ]:
function vecstats(s)
println("length : ", size(s,1), ", min : ", minimum(s), ", mean : ", mean(s), ", max : ", maximum(s), ", sum : ", sum(s))
end
In [ ]:
In [ ]:
vecstats(restijvs.s)
In [ ]:
a = grid2(200)
t = akpw(a)
frac1 = 1/5
n = size(a,1);
rest = a-t;
st = compStretches(t,rest);
aveStretch = sum(st)/nnz(rest)
@show aveStretch
targetStretch = 1/(2*log(n)/log(2))
In [ ]:
fac = aveStretch/targetStretch
heavy = rest+fac*t;
(ai,aj,av) = findnz(triu(rest))
(si,sj,sv) = findnz(triu(st))
sv = sv ./ fac
ijvs = IJVS(ai,aj,av,sv)
ijvs1 = stretchSample(ijvs,targetStretch,frac1)
ijvs2 = stretchSample(ijvs1,targetStretch,frac1)
samp1 = sparse(ijvs1.i,ijvs1.j,ijvs1.v,n,n)
samp1 = samp1 + samp1';
samp2 = sparse(ijvs2.i,ijvs2.j,ijvs2.v,n,n)
samp2 = samp2 + samp2';
In [ ]:
add = speye(n)/10^6;
@time e = eigs(lap(heavy)+add,lap(samp1+t*fac)+add, nev=1, maxiter = 10);
@show maximum(e[1])
@time e = eigs(lap(samp1)+add,lap(samp2+t*fac)+add, nev=1, maxiter = 10);
@show maximum(e[1])
In [ ]:
vecstats(ijvs.s)
vecstats(ijvs1.s)
vecstats(ijvs2.s)
In [ ]:
targetStretch
In [ ]:
In [ ]:
In [ ]:
vecstats(sampijvs.s)
In [ ]:
In [ ]:
mean(sampijvs.s)
In [ ]:
maximum(sampijvs.s)
In [ ]:
length(sampijvs.s)
In [ ]:
s2 = stretchSampleInner(sampijvs,1/2)
vecstats(s2.s)
In [ ]:
sum(s2.s .< 0.2)
In [ ]:
1/log(size(a,1))
In [ ]:
mean(restijvs.s)
In [ ]:
In [ ]:
In [ ]:
s3 = stretchSampleInner(s2,1/5)
[mean(s3.s), maximum(s3.s)]
In [ ]:
size(s3.i,1)
In [ ]:
In [ ]:
n = 10000000
p = randperm(n)
f = function()
x = b[p]
return x
end
f1 = function(b::Array{Float64,1})
n = size(b,1)
x = zeros(n)
for i in 1:n
x[p[i]] = b[i]
end
end
f3 = function(b::Array{Float64,1}, p::Array{Int64,1})
x = b[p]
return x
end
f4 = function(b::Array{Float64,1}, p::Array{Int64,1})
n = size(b,1)
x = zeros(n)
for i in 1:n
x[i] = b[p[i]]
end
end
In [ ]:
b = randn(n)
q = randperm(n)
@time x = f(b);
@time y = f1(b);
@time z = f3(b,q);
@time w = f4(b,q);
In [ ]:
@code_warntype f3(b)
In [ ]: