matrix3b.rb

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 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 pair2(n, m)
	n1 = n
	
	m.times { n = f2(n) }
	x = {0 => data(n1), 1 => data(n)}

	x[1]['wr'] = x[1]['nl'].to_f / x[0]['nl'].to_f
	return x
end

def pair(d, m)
	w = 100
	return(pair2(dense(w, d), m))
end


def dist(m, c)
		
	w = 100
	d = 0.0
	l = []
	c.times \
	{
		|i|
		
		l << pair(i.to_f / (c - 1), m)
	}
#	$stderr.puts("#{l.size} pts")
	return l
end

def dot(x, z)

	t = z['c']
	
	(z.keys - ['c']).each { |v| t += z[v] * x[v] }
	return t
		
end

def detect(x1, a, l1, c)
	
	x1 = x1.dup
	i = 0
	x1['nl'] = 1
	loop \
	{
		x2 = {}
		(a['v'] + ['wr']).each \
		{
			|vy|
			
			x2[vy] = dot(x1, a[vy].select { |k, x| (a['v'] + ['c']).member?(k) }) 
			
			return c if (x2[vy].abs > 1e10)
		}
		x2['wr2'] = (x2['wr'] - 1.0) * 50
		x2['nl'] = x1['nl'] * x2['wr']
		x1 = x2
		l1 << x1
		i += 1
		break if (x1['nl'] < 1 || i == c)
	}
	return i

end

def sum(l)

	t = 0.0
	l.each { |x| t += x }
	return t
end

def detect2(x, v, d, vi, xz, vx, l2, e)
	
	ei = (0...d.to_a.size).map { |i| d.row(i)[i] }
	n = (0...ei.size).max_by { |x| ei[x].abs }
	$stderr.puts(x)
	$stderr.puts(ei.inspect)
	$stderr.puts(['dom eig, n', ei[n], n].inspect)
	
	
#	raise if (!(ed < 0))
	
	x = x.dup
	x['wr'] = 1.0
	nl = 1.0
	File.open(fn1 = 'out1.txt', 'w').close
	File.open(fn2 = 'out2.txt', 'w').close
	File.open(fn3 = 'out3.txt', 'w').close
	
	y0 = Matrix[vx.map { |k| x[k] }].transpose - xz
	
	c0 = vi * y0
	p(c0)
	
	c1 = c0.column(0)[n]
	c1a = c1 / c1.abs
	$stderr.puts(['c1, c1a', c1, c1a].inspect)
	
	t = 1
	begin
	
		c1a *= -1 if (ei[n] < 0)
		at = v * (d ** t) * vi
		y = at * y0

		yn = y.column(0).norm
		yd = y.column(0) * c1a / yn - v.column(n)
		
		yh = Hash[[vx, yd.to_a].transpose].merge({'y_n' => yn, 'yd_n' => yd.norm})
		yh['wr'] *= 100.0
		out(fn1, yh)
		out(t.odd? ? fn2 : fn3, yh) 
		
		y += xz
		x1 = Hash[[vx, y.column(0).to_a].transpose]
		
		nl *= x1['wr']
		x1['nl'] = nl
		l2 << x1
		
		break if (yh['yd_n'] < e)
		t += 1
	end while (t < 100)
		
	f = File.open('gnuplot.cmd', 'w')
	f.puts("set colors classic;")
	k = ['wr'] + (vx - ['wr']) + ['y_n', 'yd_n']
	
	op = {}
	k.each { |x| op[x] = 'with line ' }
	op['y_n'] += 'lw 3 '
	op['yd_n'] += 'lw 3 '
	op['wr'] += 'lw 4 lc 0'
	
	[fn1, fn2, fn3].each \
	{
		|fn|
		f.puts("plot \\")
		
		k.each { |x| f.puts("'#{fn}' using (column('#{x}')) #{op[x]},\\") }
		f.puts
		f.puts("pause -1")
	}
	
	f.close
	return t
end

def out(fn, a)
	
	f = File.open(fn, 'a')
	f.puts(a.keys.join("\t")) if (f.size == 0)
	f.puts(a.values.join("\t"))
	f.close
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
#$stderr.puts(['vx', vx].inspect)
#$stderr.puts(['am', am].inspect)
#$stderr.puts(['b', b].inspect)

v, d, vi = Matrix[*am].eigen

xz = (Matrix.identity(vx.size) - Matrix[*am]).inverse * Matrix[b].transpose

l = dist(1, 100)
$c = 500

l1 = []
l2 = []

[l[95]].each \
{
	|x, i|
	
	c = detect2(x[0], v, d, vi, xz, vx, l1, 0.05)
	detect(x[0], a, l2, c)
}

File.open(fn = 'out.txt', 'w').close
[l1.size, l2.size].min.times \
{
	|i|
	x = Vector[*vx.map { |k| l1[i][k] - l2[i][k] }]
	out(fn, l2[i].merge({ 'e' => x.norm}))
}

f = File.open('gnuplot2.cmd', 'w')
f.puts('set ytics nomirror; set y2tics; plot \\')

(vx - ['wr'] + ['wr2']).each { |x| f.puts("'#{fn}' using (column('#{x}')) with line,\\") }
f.puts("'#{fn}' using (column('e')) with line axes x1y2 lw 4")
f.close