fit6.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)

	n1 = n = x['n']
	l = [n]
	while (n >= n1 && n != 1)
		
		n = f2(n)
		l << n
	end
	

	x['l'] = l
	x['ls'] = l.size
	
	x['ns'] = x['n'].to_s(2).length
	x['h'] = x['ls'].to_f / x['ns']
	x['h2'] = (x['h'] * $g).to_i
	
	return x
	
end

def next2(z)
	
	l = [z]
	p = z['p'] + 1
	
	l << adv({'n'=>z['n'] + 2**p, 'p'=>p})
	l << z.merge({'p'=>p})
	return l
end

def insert(l, x)
	l << x
end

def delete(l, j)
	z = l.delete_at(j)
	return z
end

def sum(l)

	t = 0
	l.each { |x| t += x }
	return t
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 dist(l)
	
	l2 = []
	l.each_with_index \
	{
		|x, i|
		h = x[1]['h2']
		l2[h] = [] if (l2[h].nil?)
		l2[h] << i
	}
	l1 = (0...l2.size).sort_by { |i| l2[i].nil? ? 0 : l2[i].size }
	
	return l2, l1
end

def rank(l1, l2)
	
	l1h, l1s = dist(l1)
	l2h, l2s = dist(l2)
	
	j = l1s.find { |x| !l2h[x].nil? && (l1h[x].nil? || l1h[x].size < $n) }

	j = l2h.size - 1 if (j.nil?)
	k = l2s.find { |x| !l2h[x].nil? }
	
	l = (0...l2.size).to_a
	l.sort_by! { |x| l2[x][1]['ls'] }
	k = l.find { |x| x != j }
	
	
	return l2h[j][rand(l2h[j].size)], k

end

def opt(c)
	
	l = []
	l1 = []
	
	insert(l, next2({'n'=>1, 'p'=>0}))


	puts('# ' + Time.now.to_s)
	t = Time.now.to_i

	c.times \
	{
		|i|
		$stderr.puts([i, sprintf('%.1fm', (Time.now.to_i - t) / 60.0), Time.now.to_s].join("\t")) if (i % 100 == 0)
		
		j, k = rank(l1, l)

		if (l.size > 1000) then
			
			z2 = delete(l, [j, k].max)
			z1 = delete(l, [j, k].min)
			
			l1 += [z1, z2]
			z = j < k ? z1 : z2
		else 
			z = delete(l, j)
			l1 += [z]
		end

		insert(l, next2(z[1]))
		insert(l, next2(z[2]))
		
		$stdout.flush
	}

	puts('# ' + Time.now.to_s)
	return l1.map { |x| x[1] }
end

def stat2(l, t)
	return stat(l).map { |x| x / t }
end

def d(s)
        c = s.split('').select { |x| x == '1' }.size
	d = c.to_f / s.length
        return d
        
end

def fit(l1, c)

	l1 = l1.transpose

	a = {}
	a['y'] = l1[c].to_vector()

	($c...l1.size).each \
	{
		|i|
		a["d#{i}"] = l1[i].to_vector()
	}

	ds = a.to_dataset()


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

	return [r.constant] + r.coeffs.values_at(*(a.keys - ['y']))

end

def err(z, l1, fn)
	
	z = z.dup
	c1 = z.shift

	l = []
	l1.each_with_index \
	{
		|l2, j|
		nl, = l2
		
		t = c1
		z.size.times { |i| t += z[i] * l2[i + $c] }
		
		b = fn.call(t, nl)

		l << [l2[0...$c], t, b].flatten
	}
	return l
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)
	
	l1 = ns.split(/1+/)
	l1.shift
	asdm0 = stat2(l1.map { |x| x.length }, nl)
	
	return [nl, x['ls'], x['h2'], d(ns), d(nsh), d(nsl), asdm1].flatten
end

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

def out2(fn, l)

	$stderr.puts("#{fn}\t#{l.size} pts")
	out("#{fn}e", l)
	out("#{fn}e2", l.sort_by { |x| x[1] })

end

def samplefit(c)
	
	l = opt(c)
	
	a = {}
	h = {}
	l.each \
	{ 
		|x| 
		n = x['n']
		h2 = x['h2']
		
		h[h2] = h.fetch(h2, {}).merge!({ n => nil })
		a[n] = x
	}
#	h.sort.each { |k, v| p([k, v.size]) }; exit;
	
	
	l3 = []
	
	c1 = 10
	
	h1 = h.select { |k, v| v.size >= c1 }.sort
	
	h1.each_with_index \
	{
		|kv, n|
		
		k, v = kv
		l1 = v.keys.sort

		c1.times \
		{
			|i|
			n = l1[(i.to_f / (c1 - 1) * (l1.size - 1)).to_i]
			l3 << a[n]
		}
	}
	
	l2 = l3.map { |x| data(x) }
	
	$c = 3
	z = fit(l2, $c - 1)

	out('out', l2)
	out2('out', err(z, l2, lambda { |x, nl| x / $g * nl }))
	
	return z, l3
	
end

def avg(l)
	return nil if (l.empty?)
	return sum(l) / l.size
end

def improve(z, l)
	
	l4 = []
	
	l.sort_by! { |x| x['h2'] }
	c = 20
	
	(0..(l.size - c)).each \
	{
		|j|
		l2 = []
		c.times { |i| l2 << data(l[j + i]) }
		
		l3 = l2.transpose.map { |x| avg(x.compact) }
		
		l1 = err(z, [l3], lambda { |x, nl| x / $g * nl })
		
		nl, ls, h2 = l1[0]
		
		l4 << [nl, ls, h2] + l1[0][3..4]
	}
	out2('out2', l4)
end

def compare()
	
	$n = 40
	$g = 50
	z, l = samplefit(5000)
	
#	$stderr.puts(z.inspect)
	improve(z, l)
end


l2 = compare()