outline23.rb

require 'statsample'

def f2(n)
	return n.odd? ? (n * 3 + 1) / 2 : n / 2

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["#{vy}_y"] = dot(x, z) }

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}_y")
	
	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
	}
	d = Math.sqrt(tx) * Math.sqrt(ty)
	r = txy / (d.abs < $z ? $z : d)
	e /= l.size
	return {'r' => r, 'e_a' => e, 'e_m' => m}
	
end

def fit(l, vx, vy)

	a = {}

	(vx + [vy]).each { |v| a[v] = l.map { |x| x[v] }.to_vector() }
	
	begin
		r = Statsample::Regression.multiple(a.to_dataset(), vy)
	rescue Statsample::Regression::LinearDependency
		return nil
	end
	
#	$stderr.puts(r.summary)

	z = r.coeffs.merge({'c' => r.constant})
	
	predict(l, vy, z)
	a = corr(l, vy)

	return a['r'], z

end

def read(fn)
	l = (f = File.open(fn)).readlines.map { |x| Kernel.eval(x) }
	f.close
	$stderr.puts(['read', fn, l.size].inspect)
	return l
end

def adj(l, z, ks, wt)
	ks.each \
	{
		|k|
		l.each { |x| x[k] = (x[k] - z["#{k}_a"]) * wt[k] }
	}
end

def stat(k, l)
	
	t = l.inject { |a, x| a + x }
	t2 = l.inject(0) { |a, x| a + (x ** 2) }
	
        c = l.size
        a = t.to_f / c
        z = t2.to_f / c - a ** 2
        sd = Math.sqrt(z < 0 ? 0 : z)
        raise if (sd.nan?)
	
        return {"#{k}_a" => a, "#{k}_s" => sd, 
		"#{k}_mn" => l.min, "#{k}_mx" => l.max}
end

def outd(f, l)
	f.puts('$dat << eof')
	k = l[0].keys
	
	f.puts(k.join("\t"))
	l.each { |x| f.puts(x.values.join("\t")) }

	f.puts('eof')
	return k
end

def hide(a, x)
	return (a.member?(x) && a[x].nil?)
end

def plot(f, d, k, a, t = nil)
	f.puts("# #{t}") if (!t.nil?)
	f.puts("set colors classic; set title '#{t}'; ")
#	f.puts("set key top right opaque; ")
	f.puts("set ytics nomirror; set y2tics;")
	f.puts("plot \\")
	k.each \
	{ 
		|x| 
		next if (hide(a, x))
		opt = a.fetch(x, '')
		opt = ' with line ' + opt if (!opt.include?('with') && !opt.include?('pt'))
		opt += ' lw 2 ' if (!opt.include?('lw') && !opt.include?('pt'))
		f.puts("#{d} using (column('#{x}')) #{opt} title '#{x}',\\") 
	}
	f.puts
#	f.puts("reset; pause -1;")
end

def outa(f, l, a = {}, t = nil)
	
	k = outd(f, l)
	plot(f, '$dat', k, a, t)
end

def outafn(l, a = {}, fno = nil, t = '')
	
	outa(f = File.open(fn = "gnuplot#{fno}.cmd", 'w'), l, a, t)
	f.close
	$stderr.puts([fn, l.size, t, l[0].keys].inspect)
end

def streamafn(a = nil, a1 = {}, fno = nil, t = nil)
	$f, $n = [{}, 1] if ($f.nil?)
	
	fno = $n if (fno.nil?)
	
	if (a.nil?) then
		$f[$f.keys.last]['f'].close
		$n += 1
		return
	end
	
	if (!$f.member?(fno)) then
		$f[fno] = {'fn' => (fn2 = "out#{fno}.txt"), 
			'f' => File.open(fn2, 'w')}
		$f[fno]['f'].puts(a.keys.join("\t"))
		plot(f = File.open(fn = "gnuplot#{fno}.cmd", 'w'), "'#{fn2}'", a.keys, a1, t)
		f.close
		$stderr.puts([fn, fn2, t, a.keys].inspect)
	end
	
	$f[fno]['f'].puts(a.values.join("\t"))
	$f[fno]['f'].flush
end

def keys(l, ks)
	return l.map { |x| Hash[[ks, x.values_at(*ks)].transpose] }
end

def seq(n)
	l = [n]
	
	while (n != 1)
		n = f2(n)
		
		l << n
		
	end
	return l
end

def refit(l1, l2, ks, n)
	
	k = 'hc'
	ky = "#{k}_y"
	ke = "#{k}_e"
	
	r, z= fit(l2, ks, k)
	p(r, z)
	outafn(keys(l2.sort_by { |x| x[k] }, [k, ky, ke]), {ke => 'axes x1y2 with impulse'}, n)


	predict(l1, ky, z)
	
	l1.each { |x| x['x'] = 2 }
	l2.each { |x| x['x'] = 1 }
	f = File.open("out#{n}.txt", 'w')
	
	(l3 = (l1 + l2).sort_by { |x| x[k] }).each_with_index \
	{
		|x, i|
		f.puts([i, x[k], x[ky], x[ky] - x[k], x['x']].join("\t"))
	}
	f.close
	return l3
end

	
def model(l1, ks)
	
	k = 'hc'	

	z = {}
	wt = {}
	ks.each \
	{
		|k|
		z.merge!(stat(k, l1.map { |x| x[k] }))
		wt[k] = 1.0 / z["#{k}_s"]
	}
	adj(l1, z, ks, wt)

	c = 250

	fit(l1, ks, k)
	ky = "#{k}_y"
	ke = "#{k}_e"
	l1.each { |x| x[ke] = (x[k] - x[ky]).abs }

	l1.sort_by! { |x| x[ke] }
	l2 = l1.slice!(0...c)
	
	return refit(l1, l2, ks, 1)

end

def linear1(l1)

	l1 = l1.map { |x| x.dup }

	ks = ['elo', 'a0', 'a1', 'a01', 'a1m']
	l1 = model(l1, ks)
	k = 'hc'
	ky = "#{k}_y"
	ke = "#{k}_e"
	l1.each { |x| x[ke] = (x[k] - x[ky]).abs }

	d = l1.size / 250
	p(d)
	
	l2 = []
	i = 0
	while (i < l1.size)
		l2 << l1[i...(i + d)].min_by { |x| x[ke] }
		i += d
	end
	
	
	refit(l1, l2, ks, 2)
	
end

$z = 1e-4

l1 = read('outline2.txt')
k = 'hc'
l1.each \
{
	|x|
	c = seq(x['ns'].to_i(2)).size
	x[k] = c.to_f / x['ns'].length
}
l1.sort_by! { |x| x[k] }
l1.each_with_index { |x, j| x[k] = j }

linear1(l1)

# plot 'out1.txt' using 1:3:5 with point linecolor variable pt 5 ps 0.5 title 'hc_y','' using 1:2:5 with line lt 3 lw 2 dt 3 title 'hc'