require 'statsample'


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)
	asdm0 = stat2([ns.split(/1+/)[1..-1]].compact.flatten.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 pairmap(l, m)

	return l.map { |x| pair2(x, m) }
end

def dist_d(c, w, m)

	return (0...c).map { |i| dense(w, i.to_f / (c - 1)) }
end

def read(fn)
	l = (f = File.open(fn)).readlines.map { |x| Kernel.eval(x) }
	f.close
	return l
end

def sum(l)

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

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

def corr(l, y1)

	yp = "#{y1}_p"
	ye = "#{y1}_e"
	xav = av(l.map { |x| x[y1] })
	yav = av(l.map { |x| x[yp] })

	tx = ty = txy = e = 0.0
	m = nil
	l.each \
	{
		|z|

		x = z[y1]
		y = z[yp]

		txy += (x - xav) * (y - yav)
		tx += (x - xav) ** 2
		ty += (y - yav) ** 2

		z[ye] = (x - y)
		e1 = z[ye]
		e += e1
		m = [m, e1.abs].compact.max
	}
	r = txy / (Math.sqrt(tx) * Math.sqrt(ty))
	e /= l.size
	return {'r' => r, 'e_a' => e, 'e_m' => m}

end

def dot(x, z)

	t = z['c']

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

end

def predict(l, vy, z)
	l.each { |x| x[1]["#{vy}_p"] = dot(x[0], z) }

end

def fita(l, vx, vy)

	a = {}

	vx.each { |v| a[v] = l.map { |x| x[0][v] }.to_vector() }

	vy2 = "#{vy}2"

	a[vy2] = l.map { |x| x[1][vy] }.to_vector()

	ds = a.to_dataset()

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

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

	predict(l, vy, z)

	a = corr(l.map { |x| x[1] }, vy)

	return z.merge!(a)

end

def fitb(l, vx, vy)

	a = {}

	vx.each { |v| a[v] = l.map { |x| x[0][v] }.to_vector() }

	a[vy] = l.map { |x| x[1][vy] }.to_vector()

	ds = a.to_dataset()

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

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

	predict(l, vy, z)

	a = corr(l.map { |x| x[1] }, vy)

	return z.merge!(a)

end

def fmt(a)
	a2 = {}
	a.each \
	{
		|k, v|
		a2[k] = v.is_a?(Numeric) ? sprintf("%.3g", v).to_f : v
	}
	return a2

end

def fit1(v, x, l, f)

	a = method(f).call(l, v, x)

#	a = fit(l, v, x)
	l1 = l.map { |y| y[1][x] }

	a.merge!({'mn' => l1.min, 'mx' => l1.max})

	$stderr.puts(fmt({'x' => x}.merge(a.select { |k, v| ['r', 'e_m', 'e_a', 'mn', 'mx'].member?(k) })).inspect)
	return a
end

def model(a)

	$stderr.puts(a['t'])
	(a['v'] + ['wr']).each \
	{
		|x|
		a[x] = fit1(a['v'], x, a['l'], a['f'])
	}

	return a

end

def advc(n)
	n1 = n

	l = [n]
	while (n >= n1)

		n = f2(n)
		l << n
	end
	return l.size
end

def dist_j2(fn, c)
	l = read(fn)


	l1 = []
	l2 = []
	c.times \
	{
		z1 = rand(l.size)
		ns = l[z1]['n'].to_s(2)

		z2 = rand(ns.length)
		ns[z2, 1] = (ns[z2, 1].to_i ^ 1).to_s
		n = ns.to_i(2)
		l1 << n
		l2 << {'d' => d(ns), 'c' => advc(n)}
	}
#	$stderr.puts(l2.inspect)
	return l1

end

vs = ['a1', 'a0', 'dh', 'dl', 'sd0', 'sd1', 'mx1']

d1 = lambda { dist_d(100, 100, 1) }
d2 = lambda { dist_j2('mixdb.txt', 100) }

[d1, d2].each_with_index \
{
	|d, i|
	l = pairmap(d2.call, 1)
	['fita', 'fitb'].each \
	{
		|f|
		a = model({'l' => l, 'v' => vs, 'f' => f, 't' => "#{i + 1} / #{['density', 'disordered'][i]} / #{f}"})
# a['v'].each { |x| p([x, a[x]]) }
	}
}