outline19.rb

require 'statsample'

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}_p"] = 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}_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
	}
	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
		raise 'lineardependent'
		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)
	ks.each \
	{
		|k|
		l.each { |x| x[k] -= z["#{k}_a"] }
	}
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 += ' lw 2 ' if (!opt.include?('lw'))
		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 keys(l, ks)
	return l.map { |x| Hash[[ks, x.values_at(*ks)].transpose] }
end

def linear(fn, n)
	ks = ['d', 'e', 'da', 'ea', 'dlo', 'dhi', 'elo', 'ehi', 'mp1', 'mp2',
		'mx0', 'mx1', 'mx01', 'a0', 'a1', 'a01', 'a1m']

	l1 = read(fn)

	l1.sort_by! { |x| x['hg'] }

	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)

	r, z = fit(l1, ks - ['mx01', 'a01', 'ea', 'ehi', 'elo'], 'hg')

	p([r, z])
	outafn(keys(l1, ['hg', 'hg_p']), {}, n)

end

$z = 1e-4

linear('outline2.txt', 1)
linear('outline3.txt', 2)