matrix10.rb

require 'matrix'
require 'cmath'


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, t)
	
	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 '#{t}'; 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 log2(x)
	
	return Math.log(x) if (x > 0)
	return 0
end

def exp(x, t, v, d, vi, x0)
	
	dt = d * t
	return  v * Matrix.diagonal(*(0...d.to_a.size).map { |i| CMath::exp(dt.row(i)[i]) }) * vi * x0
end

def f(x, vi, d, x0)
	a = vi * x
	b = vi * x0
	Matrix[(0...d.to_a.size).map { |i| (a.row(i)[0] / b.row(i)[0])**(1/d.row(i)[i]) * x0.row(i)[0] }].transpose
end

def load(s)
	
	raise if (s.nil?)

	d = "out/#{s}"
	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)
	return a
end


a = load(s = ARGV[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'] }
$stderr.puts(["##{s}", x1.select { |k, | k != 'l'}].inspect)

x = data(x1['l'][0]).merge({'wr' => 1.0})


x0 = Matrix[vx.map { |k| x[k] }].transpose - xz


am2 = (Matrix.identity(vx.size) - Matrix[*am]).to_a

v2, d2, vi2 = (Matrix[*am2]).eigen


m = d2.to_a.size

ei1 = (0...m).map { |i| d.row(i)[i] }
ei2 = (0...m).map { |i| d2.row(i)[i] }

n1 = (0...m).max_by { |x| ei1[x].abs }
n2 = (0...m).max_by { |x| ei2[x].abs }

$stderr.puts(['ei1', Hash[[vx, ei1].transpose]].inspect)
$stderr.puts(['ei2', Hash[[vx, ei2].transpose]].inspect)

$stderr.puts(['dom12', n1, n2].inspect)
raise if (n1 != n2)

x = (1..m).map { rand() - 0.5 }
t = rand(10)

y = exp(x, t, v2, d2, vi2, x0)

y1 = f(y, vi2, d2, x0)
y2 = Matrix[(0...m).map { |i| Math.exp(t) * x0.row(i)[0] }].transpose

$stderr.puts(['e12', (y1 - y2).column(0).norm].inspect)



l1 = []
l2 = []

ei2[n2] = -ei2[n2]

c = 110
(c + 1).times \
{
	|t|
	n = x1['l'][t]
	x = data(n).merge({'wr' => 1.0})

	y = Vector[*vx.map { |k| x[k] }]

	at = v * (t == 0 ? Matrix.identity(m) : (d ** t)) * vi
	y = at * x0

#	y += xz

	l1 << Hash[[vx, y.transpose.to_a[0]].transpose].merge({'t' => t})

	z = vi2 * y
	y2 = vi2 * x0
	
	r = (0...m).map { |i| z.row(i)[0] / y2.row(i)[0] }
	y1 = (0...m).map { |i| r[i] ** (1 / ei2[i]) * x0.row(i)[0] }
	
	l2 << Hash[[vx, y1].transpose].merge({'t' => t})
}

File.open('gnuplot.cmd', 'w').close

out('out1.txt', l1, "A^nx ##{s}")
out('out2.txt', l2, "Sum(A^nx) ##{s}")