Skip to content

Instantly share code, notes, and snippets.

@vznvzn vznvzn/matrix9b.rb
Created Dec 4, 2017

Embed
What would you like to do?
require 'matrix'
def f2(n)
return n.odd? ? (n * 3 + 1) / 2 : n / 2
end
def dense(w, d)
w2 = w - 1
a = (0...w2).to_a
s = '0' * w2
(1..(d * w - 1)).map { a.delete_at(rand(a.size)) }.each { |x| s[x, 1] = '1' }
return ('1' + s).to_i(2)
end
def stat(l)
l = [0] if (l.empty?)
t = t2 = 0
l.each \
{
|x|
t += x
t2 += x ** 2
}
c = l.size
a = t.to_f / c
z = t2.to_f / c - a ** 2
sd = Math.sqrt(z < 0 ? 0 : z)
return a, sd, l.max.to_f, l.min.to_f
end
def stat2(l, t, n)
return Hash[[["a#{n}", "sd#{n}", "mx#{n}"], stat(l)[0..2].map { |x| x / t }].transpose]
end
def d(s)
c = s.split('').select { |x| x == '1' }.size
d = c.to_f / s.length
return d
end
def count(n)
c = 0
n1 = n
l = [n]
c2 = nil
begin
n = f2(n)
l << n
c += 1
c2 = c if (c2.nil? && n <= n1)
end while n != 1
j = (0...l.size).max_by { |x| l[x] }
return {'c' => c, 'c2' => c2, 'n' => n1 , 'l' => l, 'j' => j}
end
def data(n)
ns = n.to_s(2)
nl = ns.length
m = nl / 2
nsh = ns[0..m]
nsl = ns[m..-1]
asdm1 = stat2(ns.split(/0+/).map { |x| x.length }, nl, 1)
l1 = ns.split(/1+/)
l1.shift
asdm0 = stat2(l1.map { |x| x.length }, nl, 0)
return {'n' => n, 'ns' => ns, 'nl' => nl, 'd' => d(ns), 'dh' => d(nsh), 'dl' => d(nsl)}.merge(asdm1).merge(asdm0)
end
def data2(n)
return data(n).merge(count(n))
end
def dist(c)
return (0...c).map { |i| count(dense(100, i.to_f / (c - 1))) }
end
def out(fn, l)
k = l[0].keys
f = File.open(fn, 'w')
f.puts(k.join("\t"))
l.each { |x| f.puts(x.values.join("\t")) }
f.close
f = File.open('gnuplot.cmd', 'a')
f.puts("set title '#{fn}'; plot \\")
l = []
(k - ['t']).each { |x| l << "'#{fn}' using (column('t')):(column('#{x}')) with line title '#{x}',\\" }
f.puts(l)
f.puts
f.puts("pause -1")
f.close
return l
end
def sum(l)
t = 0.0
l.each { |x| t += x }
return t
end
def int(x0, am2, vx, c, h)
yt = x0
t = 0
l = []
while (t < c)
l << Hash[[vx, yt.transpose.to_a[0]].transpose].merge({'t' => t})
yt = (h * Matrix[*am2] + Matrix.identity(vx.size)) * yt
t += h
end
return l
end
raise if (ARGV[0].nil?)
d = "out/#{ARGV[0]}"
l = Dir.glob("#{d}/*").select { |x| x =~ /out\d+.txt$/}
n = l.size
i = n - 1
a = Kernel.eval(File.open("#{d}/out#{i}.txt").readlines[0])
l = a['lb'].map { |x| x.size }
raise if (sum(l) != 0)
vx = a['v'] + ['wr']
$vxn = Hash[[vx + ['c'], (0..vx.size).to_a].transpose]
am = (1..vx.size).map { [0] * vx.size }
vx.each { |x| (vx + ['c']).each { |y| am[$vxn[x]][$vxn[y]] = a[x][y] if (!a[x][y].nil?) } }
am = am.transpose
b = am.pop
am = am.transpose
v, d, vi = Matrix[*am].eigen
xz = (Matrix.identity(vx.size) - Matrix[*am]).inverse * Matrix[b].transpose
x1 = dist(500).max_by { |x| x['c'] }
p(x1['c'])
x = data(x1['l'][0]).merge({'wr' => 1.0})
x0 = Matrix[vx.map { |k| x[k] }].transpose - xz
am2 = (Matrix[*am] - Matrix.identity(vx.size)).to_a
v2, d2, vi2 = (Matrix[*am2]).eigen
ei = (0...d2.to_a.size).map { |i| d2.row(i)[i] }
c0 = (vi2 * x0).column(0).to_a
l1 = []
l3 = []
c = 20
c.times \
{
|t|
n = x1['l'][t]
x = data(n).merge({'wr' => 1.0})
y = Vector[*vx.map { |k| x[k] }]
at = v * (d ** t) * vi
y = at * x0
# y += xz
l1 << Hash[[vx, y.transpose.to_a[0]].transpose].merge({'t' => t})
e = Matrix[(0...ei.size).map { |i| c0[i] * Math.exp(ei[i] * t.to_f) }].transpose
l3 << Hash[[vx, (v2 * e).to_a].transpose].merge({'t' => t})
}
l2 = int(x0, am2, vx, c, 1.0)
l4 = int(x0, am2, vx, c, 0.05)
File.open('gnuplot.cmd', 'w').close
l1p = out('out1.txt', l1)
l2p = out('out2.txt', l2)
l3p = out('out3.txt', l3)
l4p = out('out4.txt', l4)
f = File.open('gnuplot.cmd', 'a')
[l1p + l2p, l3p + l4p].each_with_index \
{
|l, i|
f.puts("set title '#{['l1/l2', 'l3/l4'][i]}'; plot \\")
f.puts(l)
f.puts
f.puts("pause -1")
f.puts
}
f.close
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.