NumRu::GAnalysis

# File ../../lib/numru/ganalysis/covariance.rb, line 60
def corelation(gphys0, gphys1, *dims)
  val0 = gphys0.val
  val1 = gphys1.val
  if val0.is_a?(NArrayMiss)
    mask = val0.get_mask
  else
    mask = NArray.byte(*(val0.shape)).fill!(1)
  end
  if val1.is_a?(NArrayMiss)
    mask2 = val1.get_mask
  else
    mask2 = NArray.byte(*(val1.shape)).fill!(1)
  end
  mask.mul!(mask2)
  val0 = NArrayMiss.to_nam(val0) unless val0.is_a?(NArrayMiss)
  val1 = NArrayMiss.to_nam(val1) unless val1.is_a?(NArrayMiss)
  val0 = val0.set_mask(mask)
  val1 = val1.set_mask(mask)
  p val0,val1 if $DEBUG
  gphys0 = gphys0.copy.replace_val(val0)
  gphys1 = gphys1.copy.replace_val(val1)

  covariance, ndiv = gphys0.covariance(gphys1,*dims)
  return covariance/(gphys0.stddev(*dims)*gphys1.stddev(*dims)), mask.to_type(NArray::LINT).sum(*dims)
end

# File ../../lib/numru/ganalysis/covariance.rb, line 8
def covariance(gphys0, gphys1, *dims)
  unless GPhys===gphys0 && GPhys===gphys1
    raise "gphys0 and gphys1 must be GPhys"
  end
  unless gphys0.shape == gphys1.shape
    raise "gphys0 and gphys1 must have the same shape"
  end
  units = gphys0.units*gphys1.units
  if dims.length == 0
    dims = Array.new
    gphys0.rank.times{|i| dims.push i }
  else
    dims = dims.map{|dim| gphys0.dim_index(dim) }
  end
  val0 = gphys0.val
  val1 = gphys1.val
  if val0.is_a?(NArrayMiss)
    if val1.is_a?(NArrayMiss)
      mask = val0.get_mask * val1.get_mask
      ndiv = mask.to_type(NArray::LINT).accum(*dims)
      val0 = val0.set_mask(mask)
      val1 = val1.set_mask(mask)
    else
      ndiv = val0.get_mask.to_type(NArray::LINT).accum(*dims)
      val1 = NArrayMiss.to_nam(val1,val0.get_mask)
    end
  elsif val1.is_a?(NArrayMiss)
    ndiv = val1.get_mask.to_type(NArray::LINT).accum(*dims)
    val0 = NArrayMiss.to_nam(val0,val1.get_mask)
  else
    ndiv = 1
    gphys0.shape.each_with_index{|s,i|
      ndiv *= s if dims.include?(i)
    }
  end
  val0 -= val0.accum(*dims).div!(ndiv)
  val1 -= val1.accum(*dims).div!(ndiv)
  nary = val0.mul_add(val1,*dims)
  if Float === nary
    ndiv = ndiv[0] if ndiv.is_a?(NArray)
    nary /= (ndiv-1)
    return UNumeric.new(nary, units), ndiv
  else
    nary.div!(ndiv-1)
    vary = VArray.new(nary,
                      {"long_name"=>"covariance","units"=>units.to_s},
                      "covariance")
    new_grid = gphys0.grid.delete_axes(dims, "covariance").copy
    return GPhys.new(new_grid,vary), ndiv
  end
end

# File ../../lib/numru/ganalysis/eof.rb, line 51
def eof(gphys, *args)

  unless defined?(@@EOF_engin)
    raise "SSL2 (Ruby-SSL2) or LAPACK (Ruby-LAPACK) or GSL (Ruby/GSL) must have been installed. (SSL2 or LAPACK is recommended for large computation)"
  end

  if Hash === args[-1]
    dims = args[0..-2]
    opts = args[-1]
  else
    dims = args
    opts = Hash.new
  end
  dims = dims.map{|dim| gphys.dim_index(dim) }
  n = 1
  n_lost = 1
  dims1 = Array.new
  shape1 = Array.new
  gphys.shape.each_with_index{|s,i|
    if dims.include?(i)
      n_lost *= s
    else
      n *= s
      dims1.push i
      shape1.push s
    end
  }
  new_grid = gphys.instance_variable_get("@grid").delete_axes(dims, "covariance matrix").copy
  new_index = NArray.sint(*new_grid.shape).indgen
  index = NArray.object(gphys.rank)
  index[dims] = true

  if w = (opts[:weight] || opts["weight"])
    if GPhys === w
      w = w.val
    end
    unless NArray === w
      raise "weight must be NArray of GPhys"
    end
    unless w.shape == new_grid.shape
      raise "shape of weight is invalid"
    end
    w /= w.mean
    w.reshape!(n)
  else
    if !(opts[:disable_weight]||opts["disable_weight"]) && /^lon/ =~ new_grid.coord(0).name && /^lat/ =~ new_grid.coord(1).name
      rad = NumRu::Units.new("radian")
      nlon = new_grid.coord(0).length
      lat = new_grid.coord(1).convert_units(rad).val
      w = NArray.new(lat.typecode,nlon).fill!(1) * NMath::cos(lat).reshape(1,lat.length)
      w /= w.mean
      w.reshape!(n)
    else
      w = nil
    end
  end

  ary = NArrayMiss.new(gphys.typecode, n_lost, n)
  ind_rank = dims1.length
  ind = Array.new(ind_rank,0)
  n.times{|n1|
    index[dims1] = ind
    val = gphys[*index].val
    val.reshape!(n_lost)
    val -= val.mean
    ary[true,n1] = val
    break if n1==n-1
    ind[0] += 1
    ind_rank.times{|i|
      if ind[i] == shape1[i]
        ind[i] = 0
        ind[i+1] += 1
      else
        break
      end
    }
  }
  ary.mul!(w.reshape(1,n)) if w


  nmodes = opts[:nmodes] || opts["nmodes"] || n
  case @@EOF_engin
  when "ssl2"
    print "start calc covariance matrix\n" if $DEBUG
    nary = NArray.new(gphys.typecode,n*(n+1)/2)
    nn = 0
    total_var = 0
    n.times{|n0|
      for n1 in n0...n
        nary[nn] = ary[n0].mul_add(ary[n1],0)/(n_lost-1)
        if n1==n0
          total_var += nary[nn]
        end
        nn += 1
      end
    }
    ary = nil # for GC
    print "start calc eigen vector\n" if $DEBUG
    val, vec = SSL2.seig2(nary,nmodes)
  when "lapack"
    print "start calc covariance matrix\n" if $DEBUG
    nary = NArray.new(gphys.typecode,n,n)
    total_var = 0.0
    n.times{|n0|
      nary[n0...n,n0] = (ary[true,n0...n].mul_add(ary[true,n0],0)/(n_lost-1)).get_array!
      total_var += nary[n0,n0]
    }
    ary = nil # for GC
    print "start calc eigen vector\n" if $DEBUG
    case nary.typecode
    when NArray::DFLOAT
      m, val, vec, isuppz, work, iwork, info, = NumRu::Lapack.dsyevr("V", "I", "L", nary, 0, 0, n-nmodes+1, n, 0.0, -1, -1)
      m, val, vec, = NumRu::Lapack.dsyevr("V", "I", "L", nary, 0, 0, n-nmodes+1, n, 0.0, work[0], iwork[0])
    when NArray::SFLOAT
      m, val, vec, isuppz, work, iwork, info, = NumRu::Lapack.ssyevr("V", "I", "L", nary, 0, 0, n-nmodes+1, n, 0.0, -1, -1)
      m, val, vec, = NumRu::Lapack.ssyevr("V", "I", "L", nary, 0, 0, n-nmodes+1, n, 0.0, work[0], iwork[0])
    end
    val = val[-1..0]
    vec = vec[true,-1..0]
  when "gsl"
    print "start calc covariance matrix\n" if $DEBUG
    nary = NArray.new(gphys.typecode,n,n)
    n.times{|n0|
      nary[n0...n,n0] = (ary[true,n0...n].mul_add(ary[true,n0],0)/(n_lost-1)).get_array!
      nary[n0,n0...n] = nary[n0...n,n0]
    }
    ary = nil # for GC
    print "start calc eigen vector\n" if $DEBUG
    val, vec = GSL::Eigen::symmv(nary.to_gm)
    GSL::Eigen.symmv_sort(val, vec, GSL::Eigen::SORT_VAL_DESC)
    vec = vec.to_na[0...nmodes,true].transpose(1,0)
    val = val.to_na
    total_var = val.sum
    val = val[0...nmodes]
  end

  axes = new_grid.instance_variable_get('@axes')
  axis_order = Axis.new
  axis_order.pos = VArray.new(NArray.sint(nmodes).indgen(1),
                              {}, "mode")
  axes << axis_order
  new_grid = Grid.new(*axes)
  vec /= w if w
  vec.reshape!(*new_grid.shape)
  vec *= NMath::sqrt( val.reshape( *([1]*(axes.length-1)+[nmodes]) ) )
  va_eof = VArray.new(vec,
                      {"long_name"=>"EOF vector","units"=>gphys.units.to_s },
                      "EOF")
  eof = GPhys.new(new_grid, va_eof)

  va_rate = VArray.new(val.div!(total_var),
                       {"long_name"=>"EOF contribution rate", "units"=>"1" },
                       "rate")
  rate = GPhys.new(Grid.new(axis_order), va_rate)

  return [eof, rate]
end

# File ../../lib/numru/ganalysis/eof.rb, line 209
def eof2(gphys1, gphys2, *args)
  if Hash === args[-1]
    opts = args[-1]
  else
    opts = Hash.new
  end
  raise ArgumentError, "The 1st arg must be a GPhys of rank 1: arg1 = #{gphys1.inspect}" unless gphys1.rank==1
  raise ArgumentError, "The 2nd arg must be a GPhys of rank 1: arg2 = #{gphys2.inspect}" unless gphys2.rank==1
  raise ArgumentError, "The 1st and 2nd args must have the same length: #{gphys1.length}!=#{gphys2.length}" unless gphys2.rank==1
  nam = NArrayMiss.new(gphys1.typecode, 2, gphys1.length)
  nam[0,true] = gphys1.val
  nam[1,true] = gphys2.val
  gphys = GPhys.new(Grid.new(Axis.new.set_pos(VArray.new(NArray[0,1],{},"var")),
                             gphys1.axis(0)),
                    VArray.new(nam,gphys1.data.attr_copy,gphys1.name))
  eof(gphys, 1, opts)
end

# File ../../lib/numru/ganalysis/histogram.rb, line 15
def histogram(gphys0,opts=Hash.new)
  unless HistogramGSL
    raise "gsl is necessary to use this method"
  end
  unless GPhys === gphys0
    raise "gphys0 (1st arg) must be GPhys"
  end
  unless Hash === opts
    raise "opts (2nd arg) must be Hash"
  end
  if opts["bins"]
    bins = opts["bins"]
    unless (bins.is_a?(NArray) || bins.is_a?(Array))
      raise(TypeError, "option 'bins' must be Array or NArray")
    end
    bins = bins.to_gslv if bins.is_a?(NArray)
    hist = GSL::Histogram.alloc(bins)
  else
    nbins = opts["nbins"] || gphys0.total/500
    nbins = 10 if nbins < 10
    min = opts["min"] || gphys0.min.val
    max = opts["max"] || gphys0.max.val
    if log_bins = (opts["log_bins"] && (min > 0))
      min = Math.log10(min)
      max = Math.log10(max)
    end
    hist = GSL::Histogram.alloc(nbins,[min,max])
  end
  val = gphys0.val
  val = val.get_array![val.get_mask!] if NArrayMiss === val
  val = NMath.log10(val) if log_bins
  hist.increment(val)

  bounds = hist.range.to_na
  bounds = 10 ** bounds if log_bins
  center = (bounds[0..-2]+bounds[1..-1])/2
  cell_width = (bounds[1..-1] - bounds[0..-2]) / 2
  name = gphys0.name
  attr = gphys0.data.attr_copy
  bounds = VArray.new(bounds, attr, name)
  center = VArray.new(center, attr, name)
  axis = Axis.new(true)
  axis.set_cell(center, bounds, name)
  axis.set_pos_to_center

  bin = hist.bin.to_na
  bin /= cell_width if opts["log_bins"]
  bin = VArray.new(bin,
                   {"long_name" => (log_bins ? "number per unit bin width" : "number in bins"), "units"=>"1"},
                   "bin")
  new_gphys = GPhys.new(Grid.new(axis), bin)
  new_gphys.set_att("mean",[hist.mean])
  new_gphys.set_att("standard_deviation",[hist.sigma])
  return new_gphys
end

# File ../../lib/numru/ganalysis/histogram.rb, line 72
def histogram2D(gphys0, gphys1, opts=Hash.new)
  unless HistogramGSL
    raise "gsl is necessary to use this method"
  end
  unless GPhys === gphys0
    raise "gphys0 (1st arg) must be GPhys"
  end
  unless GPhys === gphys1
    raise "gphys1 (2nd arg) must be GPhys"
  end
  unless Hash === opts
    raise "opts (3nd arg) must be Hash"
  end

  nbins0 = opts["nbins0"] || gphys0.total/500
  nbins0 = 10 if nbins0 < 10
  nbins1 = opts["nbins1"] || gphys1.total/500
  nbins1 = 10 if nbins1 < 10

  min0 = opts["min0"] || gphys0.min.val
  max0 = opts["max0"] || gphys0.max.val
  min1 = opts["min1"] || gphys1.min.val
  max1 = opts["max1"] || gphys1.max.val

  hist = GSL::Histogram2d.alloc(nbins0,[min0,max0],nbins1,[min1,max1])
  val0 = gphys0.val
  val1 = gphys1.val
  mask = nil
  if NArrayMiss === val0
    mask = val0.get_mask!
    val0 = val0.get_array!
  end
  if NArrayMiss === val1
    if mask
      mask = mask & val1.get_mask!
    else
      mask = val1.get_mask!
    end
    val1 = val1.get_array!
  end
  if mask
    val0 = val0[mask]
    val1 = val1[mask]
  end
  hist.increment(val0.to_gslv, val1.to_gslv)

  bounds0 = hist.xrange.to_na
  center0 = (bounds0[0..-2]+bounds0[1..-1])/2
  name = gphys0.name
  attr = gphys0.data.attr_copy
  bounds0 = VArray.new(bounds0, attr, name)
  center0 = VArray.new(center0, attr, name)
  axis0 = Axis.new(true)
  axis0.set_cell(center0, bounds0, name)
  axis0.set_pos_to_center

  bounds1 = hist.yrange.to_na
  center1 = (bounds1[0..-2]+bounds1[1..-1])/2
  name = gphys1.name
  attr = gphys1.data.attr_copy
  bounds1 = VArray.new(bounds1, attr, name)
  center1 = VArray.new(center1, attr, name)
  axis1 = Axis.new(true)
  axis1.set_cell(center1, bounds1, name)
  axis1.set_pos_to_center

  bin = hist.bin.to_na.reshape!(nbins1,nbins0).transpose(1,0)
  bin = VArray.new(bin,
                   {"long_name"=>"number in bins", "units"=>"1"},
                   "bin")
  new_gphys = GPhys.new(Grid.new(axis0,axis1), bin)
  new_gphys.set_att("mean0",[hist.xmean])
  new_gphys.set_att("standard_deviation0",[hist.xsigma])
  new_gphys.set_att("mean1",[hist.ymean])
  new_gphys.set_att("standard_deviation1",[hist.ysigma])
  new_gphys.set_att("covariance",[hist.cov])
  return new_gphys
end