#  Histogram class  ver0.23
#       written by TAKAHASHI Hitoshi <thitoshi@ne.scphys.kyoto-u.ac.jp>
#     date: 1999/04/06

class Histogram
  def initialize(title, xch, xmin, xmax, data = nil)
    raise "invalid number of channels" if xch <= 0
    raise "upper limit is less than lower" if xmin >= xmax
    @title = title.dup
    @xch, @xmin, @xmax = xch, xmin, xmax
    @xbin = Float(@xmax - @xmin) / @xch
    @contents = []
    if data == nil
      clear
    else
      type_check(data, Array)
      raise "number of initial data must be equal to or more than that of channels" if data.size < @xch
      @contents[0] = 0
      @contents[1..@xch] = data[0..@xch-1]
      @contents[@xch + 1] = 0
      @entry = sum
    end
  end

  attr_reader :xch, :xmin, :xmax
  attr_accessor :title, :entry

  def [](i)
    @contents[i]
  end

  def []=(i, d)
    @contents[i] = d
  end

  protected :entry=, :[]=

  def dup
    hist = Histogram.new(@title, @xch, @xmin, @xmax, @contents[1..@xch])
    hist.entry = @entry
    hist[0] = @contents[0]
    hist[@xch + 1] = @contents[@xch + 1]
    hist
  end

  def ===(other)
    if other.kind_of?(Histogram) && other.xch == @xch
      true
    else
      false
    end
  end

  def +(other)
    type_check(other, self.type)
    if other.xch != @xch
      raise "channel unmatch"
    end
    hist = dup
    for i in 0..@xch+1
      hist[i] += other[i]
    end
    hist.entry += other.entry
    hist
  end

  def -(other)
    type_check(other, self.type)
    if other.xch != @xch
      raise "channel unmatch"
    end
    hist = dup
    for i in 0..@xch+1
      hist[i] -= other[i]
    end
    hist.entry += other.entry
    hist
  end

  def *(times)
    type_check(times, Numeric)
    hist = dup
    for i in 0..@xch+1
      hist[i] *= times
    end
    hist
  end

  def /(times)
    type_check(times, Numeric)
    hist = dup
    for i in 0..@xch+1
      hist[i] /= times
    end
    hist
  end

  def clear
    for i in 0..@xch+1
      @contents[i] = 0
    end
    @entry = 0
  end

  def contents(x)
    @contents[channel(x)]
  end

  def to_a
    @contents[1..@xch]
  end

  def fill(x, d = 1)
    @contents[channel(x)] += d
    @entry += 1
  end

  def channel(x)
    if x > @xmax  # overflow
      @xch + 1
    elsif x < @xmin  # underflow
      0
    else
      ((x - @xmin) / @xbin).to_i + 1
    end
  end

  def underflow
    @contents[0]
  end

  def overflow
    @contents[@xch + 1]
  end

  def sum
    _sum = 0
    for i in 1..@xch
      _sum += @contents[i]
    end
    _sum
  end

  def mean
    _mean = 0
    for i in 1..@xch
      x = @xmin + @xbin * (i - 0.5)
      _mean += x * @contents[i]
    end
    _mean / sum
  end

  def sigma
    _mean = mean
    _sigma = 0
    for i in 1..@xch
      x = @xmin + @xbin * (i - 0.5)
      _sigma += (x - _mean) ** 2 * @contents[i]
    end
    Math.sqrt(_sigma / sum)
  end

  def min
    _min = @contents[1]
    for i in 2..@xch
      _min = @contents[i] if @contents[i] < _min
    end
    _min
  end

  def max
    _max = @contents[1]
    for i in 2..@xch
      _max = @contents[i] if @contents[i] > _max
    end
    _max
  end

  def to_s
    str = "# " + @title + "\n"
    for i in 1..@xch
      x = @xmin + @xbin * (i - 0.5)
      str << x.to_s << " " << @contents[i].to_s << "\n"
    end
    str
  end

  def each
    for i in 1..@xch
      yield @contents[i]
    end
  end

  def each_ch
    (@xmin + @xbin / 2).step(@xmax, @xbin) do |x|
      yield x
    end
  end

  def type_check(arg, expected)
    if ! arg.kind_of? expected
      raise TypeError, "wrong argument type #{arg.type} (expected #{expected})"
    end
  end

  private :type_check
end


class IntHisto < Histogram
  def initialize(title, xmin, xmax, xstep, data = nil)
    type_check(xmin, Integer)
    type_check(xmax, Integer)
    type_check(xstep, Integer)
    raise "upper limit is less than lower" if xmin >= xmax
    raise "invalid step" if xstep <= 0 || xstep > xmax - xmin
    @title = title.dup
    @xbin, @xmin = xstep, xmin
    @xch = (xmax - @xmin) / @xbin + 1
    @xmax = (@xch - 1) * @xbin + @xmin
    @contents = []
    if data == nil
      clear
    else
      type_check(data, Array)
      raise "number of initial data must be equal to or more than that of channels" if data.size < @xch
      @contents[0] = 0
      @contents[1..@xch] = data[0..@xch-1]
      @contents[@xch + 1] = 0
      @entry = sum
    end
  end

  def dup
    hist = IntHisto.new(@title, @xmin, @xmax, @xbin, @contents[1..@xch])
    hist.entry = @entry
    hist[0] = @contents[0]
    hist[@xch + 1] = @contents[@xch + 1]
    hist
  end

  def channel(x)
    type_check(x, Integer)
    if x > @xmax  # overflow
      @xch + 1
    elsif x < @xmin  # underflow
      0
    else
      (x - @xmin) / @xbin + 1
    end
  end

  def to_s
    str = "# " + @title + "\n"
    for i in 1..@xch
      x = @xmin + @xbin * (i - 1)
      str << x.to_s << " " << @contents[i].to_s << "\n"
    end
    str
  end

  def mean
    _mean = 0
    for i in 1..@xch
      x = @xmin + @xbin * (i - 1)
      _mean += x * @contents[i]
    end
    _mean / sum
  end

  def sigma
    _mean = mean
    _sigma = 0
    for i in 1..@xch
      x = @xmin + @xbin * (i - 1)
      _sigma += (x - _mean) ** 2 * @contents[i]
    end
    Math.sqrt(_sigma / sum)
  end

  def each_ch
    @xmin.step(@xmax, @xbin) do |x|
      yield x
    end
  end
end


class Histogram2
  def initialize(title, xch, xmin, xmax, ych, ymin, ymax)
    raise "invalid number of x channel" if xch <= 0
    raise "upper limit of x is less than lower" if xmin >= xmax
    raise "invalid number of y channel" if ych <= 0
    raise "upper limit of y is less than lower" if ymin >= ymax
    @title = title.dup
    @xch, @xmin, @xmax = xch, xmin, xmax
    @ych, @ymin, @ymax = ych, ymin, ymax
    @xbin = Float(@xmax - @xmin) / @xch
    @ybin = Float(@ymax - @ymin) / @ych
    @contents = []
    (@ych + 2).times do |i|
      @contents[i] = []
    end
    clear
  end

  attr_reader :xch, :xmin, :xmax, :ych, :ymin, :ymax
  attr_accessor :title, :entry

  def [](i, j)
    @contents[i][j]
  end

  def []=(i, j, d)
    @contents[i][j] = d
  end

  protected :entry=, :[]=

  def dup
    hist = Histogram2.new(@title, @xch, @xmin, @xmax, @ych, @ymin, @ymax)
    for i in 0..@ych+1
      for j in 0..@xch+1
	hist[i, j] = @contents[i][j]
      end
    end
    hist.entry = @entry
    hist
  end

  def ===(other)
    if other.kind_of?(Histogram2) && other.xch == @xch && other.ych == @ych
      true
    else
      false
    end
  end

  def +(other)
    type_check(other, self.type)
    if other.xch != @xch || other.ych != @ych
      raise "channel unmatch"
    end
    hist = dup
    for i in 0..@ych+1
      for j in 0..@xch+1
	hist[i, j] += other[i, j]
      end
    end
    hist.entry += other.entry
    hist
  end

  def -(other)
    type_check(other, self.type)
    if other.xch != @xch || other.ych != @ych
      raise "channel unmatch"
    end
    hist = dup
    for i in 0..@ych+1
      for j in 0..@xch+1
	hist[i, j] -= other[i, j]
      end
    end
    hist.entry += other.entry
    hist
  end

  def *(times)
    type_check(times, Numeric)
    hist = dup
    for i in 0..@ych+1
      for j in 0..@xch+1
	hist[i, j] *= times
      end
    end
    hist
  end

  def /(times)
    type_check(times, Numeric)
    hist = dup
    for i in 0..@ych+1
      for j in 0..@xch+1
	hist[i, j] /= times
      end
    end
    hist
  end

  def clear
    for i in 0..@ych+1
      for j in 0..@xch+1
	@contents[i][j] = 0
      end
    end
    @entry = 0
  end

  def contents(x, y)
    xch, ych = xchannel(x), ychannel(y)
    @contents[ych][xch]
  end

  def to_a
    array = []
    for i in 0..@ych-1
      array[i] = @contents[i][1..@xch]
    end
    array
  end

  def fill(x, y, d = 1)
    xch, ych = xchannel(x), ychannel(y)
    @contents[ych][xch] += d
    @entry += 1
  end

  def xchannel(x)
    if x > @xmax  # overflow
      @xch + 1
    elsif x < @xmin  # underflow
      0
    else
      ((x - @xmin) / @xbin).to_i + 1
    end
  end

  def ychannel(y)
    if y > @ymax  # overflow
      @ych + 1
    elsif y < @ymin  # underflow
      0
    else
      ((y - @ymin) / @ybin).to_i + 1
    end
  end

  def projection_x(ch1, ch2, title)
    data = []
    for i in 0..@xch-1
      data[i] = 0
      ch1.upto(ch2) do |ch|
	data[i] += @contents[ch][i + 1]
      end
    end
    Histogram.new(title, @xch, @xmin, @xmax, data)
  end

  def projection_y(ch1, ch2, title)
    data = []
    for i in 0..@ych-1
      data[i] = 0
      ch1.upto(ch2) do |ch|
	data[i] += @contents[i + 1][ch]
      end
    end
    Histogram.new(title, @ych, @ymin, @ymax, data)
  end

  private :projection_x, :projection_y

  def xproj
    projection_x(1, @ych, @title + " x projection")
  end

  def yproj
    projection_y(1, @xch, @title + " y projection")
  end

  def xslice(y)
    ch = ychannel(y)
    projection_x(ch, ch, @title + " x slice")
  end

  def yslice(x)
    ch = xchannel(x)
    projection_y(ch, ch, @title + " y slice")
  end

  def xband(y1, y2)
    if y1 < y2
      projection_x(ychannel(y1), ychannel(y2), @title + " x band")
    else
      projection_x(ychannel(y2), ychannel(y1), @title + " x band")
    end
  end

  def yband(x1, x2)
    if x1 < x2
      projection_y(ychannel(x1), ychannel(x2), @title + " y band")
    else
      projection_y(ychannel(x2), ychannel(x1), @title + " y band")
    end
  end

  def sum
    _sum = 0
    for i in 1..@ych
      for j in 1..@xch
	_sum += @contents[i][j]
      end
    end
    _sum
  end

  def min
    _min = @contents[1][1]
    for i in 1..@ych
      for j in 1..@xch
	_min = @contents[i][j] if @contents[i][j] < _min
      end
    end
    _min
  end

  def max
    _max = @contents[1][1]
    for i in 1..@ych
      for j in 1..@xch
	_max = @contents[i][j] if @contents[i][j] > _max
      end
    end
    _max
  end

  def to_s
    str = "# " + @title + "\n"
    for i in 1..@ych
      y = @ymin + @ybin * (i - 0.5)
      for j in 1..@xch
	x = @xmin + @xbin * (j - 0.5)
	str << x.to_s << " " << y.to_s << " " << @contents[i][j].to_s << "\n"
      end
    end
    str
  end

  def each
    for i in 1..@ych
      for j in 1..@xch
	yield @contents[i][j]
      end
    end
  end

  def each_ch
    (@ymin + @ybin / 2).step(@ymax, @ybin) do |y|
      (@xmin + @xbin / 2).step(@xmax, @xbin) do |x|
	yield x, y
      end
    end
  end

  def type_check(arg, expected)
    if ! arg.kind_of? expected
      raise TypeError, "wrong argument type #{arg.type} (expected #{expected})"
    end
  end

  private :type_check
end


class IntHisto2 < Histogram2
  def initialize(title, xmin, xmax, xstep, ymin, ymax, ystep)
    type_check(xmin, Integer)
    type_check(xmax, Integer)
    type_check(xstep, Integer)
    type_check(ymin, Integer)
    type_check(ymax, Integer)
    type_check(ystep, Integer)
    raise "upper limit of x is less than lower" if xmin >= xmax
    raise "upper limit of y is less than lower" if ymin >= ymax
    raise "invalid step" if xstep <= 0 || xstep > xmax - xmin
    raise "invalid step" if ystep <= 0 || ystep > ymax - ymin
    @title = title.dup
    @xbin, @xmin = xstep, xmin
    @ybin, @ymin = ystep, ymin
    @xch = (xmax - @xmin) / @xbin + 1
    @ych = (ymax - @ymin) / @ybin + 1
    @xmax = (@xch - 1) * @xbin + @xmin
    @ymax = (@ych - 1) * @ybin + @ymin
    @contents = []
    (@ych + 2).times do |i|
      @contents[i] = []
    end
    clear
  end

  def dup
    hist = IntHisto2.new(@title, @xmin, @xmax, @xbin, @ymin, @ymax, @ybin)
    for i in 0..@ych+1
      for j in 0..@xch+1
	hist[i, j] = @contents[i][j]
      end
    end
    hist.entry = @entry
    hist
  end

  def xchannel(x)
    type_check(x, Integer)
    if x > @xmax  # overflow
      @xch + 1
    elsif x < @xmin  # underflow
      0
    else
      (x - @xmin) / @xbin + 1
    end
  end

  def ychannel(y)
    type_check(y, Integer)
    if y > @ymax  # overflow
      @ych + 1
    elsif y < @ymin  # underflow
      0
    else
      (y - @ymin) / @ybin + 1
    end
  end

  def projection_x(ch1, ch2, title)
    data = []
    for i in 0..@xch-1
      data[i] = 0
      ch1.upto(ch2) do |ch|
	data[i] += @contents[ch][i + 1]
      end
    end
    IntHisto.new(title, @xmin, @xmax, @xbin, data)
  end

  def projection_y(ch1, ch2, title)
    data = []
    for i in 0..@ych-1
      data[i] = 0
      ch1.upto(ch2) do |ch|
	data[i] += @contents[i + 1][ch]
      end
    end
    IntHisto.new(title, @ymin, @ymax, @ybin, data)
  end

  def to_s
    str = "# " + @title + "\n"
    for i in 1..@ych
      y = @ymin + @ybin * (i - 1)
      for j in 1..@xch
	x = @xmin + @xbin * (j - 1)
	str << x.to_s << " " << y.to_s << " " << @contents[i][j].to_s << "\n"
      end
    end
    str
  end

  def each_ch
    @ymin.step(@ymax, @ybin) do |y|
      @xmin.step(@xmax, @xbin) do |x|
	yield x, y
      end
    end
  end
end