dist2b.rb

require 'statsample'


def f2(n)

	n = (n * 3 + 1) / 2 while (n.odd?)
	n /= 2 while (n.even?)
	return n

end

def adv(x, l1)

	x['n'] = x['nb'].to_i(2)
	n1 = n = x['n']
	l = [n]
	while (n != 1)
		
		n = f2(n)
		l << n
	end
	

	x['l'] = l
	x['ls'] = l.size
	
	x['ns'] = x['n'].to_s(2).length
	x['m'] = (0...l.size).max_by { |x| l[x] }

	x['y'] = (l.max.to_s(2).length.to_f / x['ns'] * 100).to_i
	
	l1 << x
	
	return x
	
end

def adv1(x, l1, z)
	
	x1 = adv(x, l1)
	
	x = data(x1)
	
	x = predict1(x, z)

	x1['y'] = (x['y'] - x['y_p']).abs
	
	x['ns'] = x1['ns']
	out1(x, x.keys) if (rand(50) == 0)

	return x1

end

def out1(a, l)
	if ($f.nil?) then
		$f = File.open('data.txt', 'w')
		$f.puts(l.join("\t"))
	end
	$f.puts(a.values_at(*l).join("\t"))
	$f.flush
	
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
end

def stat2(l, t, n)
	return Hash[[["a#{n}", "sd#{n}", "mx#{n}"], stat(l).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(x)
	
	n = x['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' => x['n'], 'y' => x['y'], 'd' => d(ns), 'dh' => d(nsh), 'dl' => d(nsl)}.merge(asdm0).merge(asdm1)
end

def out(fn, l)
	
	f = File.open("#{fn}.txt", 'w')
	f.puts(l[0].keys.join("\t"))
	l.each { |x| f.puts(x.values.join("\t")) }
	f.close
	
end

def dist(w, advx, c)
	
	$stdout.puts("#{w} width")
	l = []
	l1 = []
	
	20.times { l << advx.call({'nb' => (0...w).map { rand(2).to_s }.join }, l1) }

	seen = {0 => nil}

	i = n = 0
	while (n < c)

		i += 1
		l.sort_by! { |x| -x['y'] }

		t = l.size / 10
		case i % 3
			when 0
				x = l[rand(t)]
				s = x['nb']
				r = rand(s.length)
				s[r, 1] = (s[r, 1].to_i  ^ 1).to_s
				
			when 1
				x = l[rand(t)]
				y = l[rand(t)]
				sx = x['nb']
				sy = y['nb']
				
				s = ''
				w.times \
				{
					|j|
					s[j, 1] = (rand(2) == 0) ? sx[j, 1] : sy[j, 1]
				}
				
			when 2
				x = l[rand(t)]
				y = l[rand(t)]
				sx = x['nb']
				sy = y['nb']
				s = ''
				r = rand(w + 1)
				w.times \
				{
					|j|
					s[j, 1] = (j < r) ? sx[j, 1] : sy[j, 1]
				}
		end
		
		next if (seen.member?(s))
		seen[s] = nil
		
		l.pop if (l.size > 500)
		l << advx.call({'nb' => s}, l1)
		n += 1
	end

	$stderr.puts("#{i - n} dups")
	return l1
	

end

def hist(l1)
	
	$h = {}
	l1.each \
	{
		|x|
		
		$h[x['y']] = [] if (!$h.member?(x['y']))
		$h[x['y']] << x
	}
	
	l = []


	$h.sort.each \
	{
		|k, v|
		p([k, v.size])
	}

	$c = 10
	$h.sort.select { |k, v| v.size >= $c }.each \
	{ 
		|k, v|
		$c.times { l << data(v.delete_at(rand(v.size))) }
	}
	$stderr.puts("#{l.size} pts")

	return l
	
end

def fit(l, y, lx)

	a = {}

	(lx + [y]).each { |x| a[x] = l.map { |b| b[x] }.to_vector() }

	ds = a.to_dataset()

	r = Statsample::Regression.multiple(ds, y)
	$stderr.puts(r.summary)

	return r.coeffs.merge({'c' => r.constant})

end

def adv2(n, c)

	l = []

	c.times \
	{
		n = f2(n)
		l << n
	}
	return l
	
end

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

def sum(l)

	t = 0
	l.each { |x| t += x }
	return t
end
	
def av(l)
	return nil if (l.empty?)
	return sum(l) / l.size
end

def predict1(x, z)
	t = dot(x, z)
	
	ns1 = x['n'].to_i.to_s(2).length
	
	t = 100 if (t < 100)
	l2 = [t.to_f * ns1]
	
	adv2(x['n'].to_i, 10).each \
	{
		|n|
		x2 = data({'n' => n})
		ns = n.to_s(2).length
		
		t = dot(x2, z)
		t = 100 if (t < 100)
		
		t *= ns.to_f
		l2 << t
	}
	x['y_p'] = av(l2) / ns1
	x['y_p'] = 100 if (x['y_p'] < 100)
	return x
end

def predict(z, l1)
	l1.each_with_index \
	{
		|x, i|
		predict1(x, z)
	}
end

l2 = (1..3).map { hist(dist(50, method(:adv), 50 * 1e3)) }.flatten

z = fit(l2, 'y', l2[0].keys - ['y', 'n'])

predict(z, l2)

out('out1', l2.sort_by { |x| x['y'] })


w = 50
loop { dist(w += 5, lambda { |x, l1| adv1(x, l1, z) }, 2 * 1e3) }