matrix6.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 count(n)
	
	c = 0
	n1 = n
	l = [n]
	begin
		n = f2(n) 
		l << n
		c += 1
	end while n > n1
	j = (0...l.size).max_by { |x| l[x] }
	
	return {'c' => c, 'n' => n1}
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(200, i.to_f / (c - 1))) }

end

def sum(l)

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

def real(x)
	raise x.imag.to_s if (x.imag > 1e-15)
	return x.real
end

def fmt(n)
	
	
	n.real? ? sprintf("%.3g", n) : sprintf("(%.3g %.3gi)", n.real, n.imag)
end

def detect2(n, v, d, vi, xz, vx)
	
	x = data(n)
	
	x['wr'] = 1.0
	nl = 1.0
	
	y0 = Matrix[vx.map { |k| x[k] }].transpose - xz
	t = 1
	y = nil
	
	begin
	
		at = v * (d ** t) * vi
		y = at * y0

		yn = y.column(0).norm
		yd = y.column(0) / yn - v.column(0)
		
		yh = Hash[[vx, yd.to_a].transpose].merge({'y_n' => yn, 'yd_n' => yd.norm})
		
		y2 = y + xz
		x1 = Hash[[vx, y2.column(0).to_a].transpose]
		
		nl *= x1['wr']
		x1['nl'] = nl
		break if (nl.abs < 1.0)
		
		t += 1
	end while (t < 100)
	
	return {'t' => t, 'yn' => (y - y0).column(0).norm}
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

def sample(n)

	l2 = []
	n.times \
	{
		l = dist(500)
		l.sort_by! { |x| -x['c'] }
		l2 += l[0...50]
	}
		
	l = []
	l1 = l2.map { |x| x['c'] }
	mn = l1.min
	mx = l1.max
	(0..100).each \
	{
		|i|
		d = i.to_f / 100

		x1 = nil
		begin 
			j = (0...l2.size).min_by { |x| ((l2[x]['c'].to_f - mn) / (mx - mn) - d).abs }
			break if (j.nil?)
			x1 = l2.delete_at(j)
		end while (!l.select { |x| x['n'] == x1['n'] }.empty?)
		l << x1 if (!x1.nil?)
	}
	
	return l.sort_by { |x| -x['c'] }
	
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

File.open(fn = 'out.txt', 'w').close

sample(10).each_with_index \
{
	|x, i|
	
	out(fn, {'i' => i, 'c' => x['c'], 'd' => d(x['n'].to_s(2))}.merge(detect2(x['n'], v, d, vi, xz, vx)))
}