fit.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]
	l2 = []
	while (n >= n1 && n != 1)
		
		l2 << n if (!$seen.member?(n))
			
		n = f2(n)
		l << n
	end
	

	x['l'] = l
	x['ls'] = l.size
	
	x['l2'] = l2
	x['l2s'] = l2.size

	x['ns'] = x['n'].to_s(2).length
	
	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)
	x[1]['l2'].each { |x| $seen[x] = nil }
	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, l.min.to_f
end

def rank(l, w1)
	
	
	l2 = []
	
	l.size.times \
	{ 
		|x| 
		l1 = l[x][1].values_at(*w1)
		a = {}
		w1.size.times \
		{
			|i|
			a[w1[i]] = l1[i]
		}
		l2 << [x, a]
	}
	
	w1.each \
	{
		|k|
		l1 = l2.map { |x| x[1][k] }
		
		a, sd, = stat(l1)
		
		sd = 1 if (sd == 0)
		l2.each \
		{
			|x|
			x[1][k + 'b'] = (x[1][k] - a) / sd
		}
	}

	b = w1.map { |x| x + 'b' }
	l2.each \
	{
		|x|
		x[1]['z'] = sum(x[1].values_at(*b))
	}
	
	l2.sort_by! { |x| [x[1]['z'], x[1]['d']] }
	l2.reverse!
		
	j = l2[rand([l2.size, 100].min)][0]
	k = l2[-1][0]

	return j, k
	
end

def opt(w1, c)
	
	l = []

	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(l, w1)
			
		if (l.size > 1000) then
			
			z2 = delete(l, [j, k].max)
			z1 = delete(l, [j, k].min)
			
			z = j < k ? z1 : z2
		else 
			z = delete(l, j)
		end

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

	puts('# ' + Time.now.to_s)
	return l.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 data(l, ls)
	l2 = []
	l.each \
	{
		|x|
		
		ns = x.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)
		
		z = ls.to_f / nl
		
		l2 << [nl, ls, z, d(ns), d(nsh), d(nsl), asdm0].flatten
	}
	return l2
end

def sample()
	
	$seen = {}

	l = opt(['l2s'], 2000)

	h = {}
	l.each \
	{ 
		|x| 
		x['l'].reverse.each_with_index { |y, i| h[i] = h.fetch(i, []) + [y] }
	}

	l2 = []
	h1 = h.select { |k, v| v.size >= 4 }.each \
	{
		|k, v|
		l1 = v.sort
		l2 += data(l1[0..1] + l1[-2..-1], k)
	}

	return l2
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, c, fn)
	
	e1 = e2 = 0
	z = z.dup
	c1 = z.shift
	i2 = 0
	l = []
	l1.each_with_index \
	{
		|l2, j|
		nl, ls = l2
		
		t = c1
		z.size.times { |i| t += z[i] * l2[i + $c] }
		
		a = fn.call(l2[c], nl)
		b = fn.call(t, nl)
		
		e1 += (l2[c] - t).abs
		e2 += (a - b).abs
		
		l << [nl, ls, l2[c], t, a, b, e1 / (j + 1), e2 / (j + 1)]
	}
	return l
end


l2 = sample()

$c = 3
z = fit(l2, $c - 1)
l1 = err(z, l2, $c - 1, lambda { |x, nl| x * nl })
l1.each { |x| puts(x.join("\t")) }