lib/numru/gphys/gphys.rb in gphys-1.2.2.1 vs lib/numru/gphys/gphys.rb in gphys-1.4.3
- old
+ new
@@ -68,10 +68,76 @@
USAGE
See the manual of ((|GPhys::IO.each_along_dims_write|)).
+---GPhys.join_md_nocheck(gpnary)
+ Join multiple GPhys objects that are ordered perfectly in a NArray.
+
+ LIMITATION (as of 2013-03-04)
+ * joining assoc_coords is yet to be supported; Currently
+ assoc_coords are ignored if any.
+
+ ARGUMENT
+ * gpnarray [NArray of GPhys] having the same rank with that of
+ its component GPhys objects. multiple GPhys objects are joined
+ along the dimension with multiple elements (the order is kept).
+
+ RETURN VALUE
+ * a GPhys
+
+---GPhys.join_md(gpnary)
+ Join multiple GPhys objects (ordered in a NArray).
+
+ Like GPhys.join_md_nocheck but it supports insersion
+ of omitted 1-element dimensions and transpose for
+ gpnary (the input NArray). It means that the rank of gpnary
+ can be smaller than that of its compoent GPhys objects, and
+ the order of dimensions can be arbitrary. Also,
+ the order of coordinate values along each dimension does not
+ have to be monotonic; the method supports sorting and spliting
+ along dimensions. For example, if gpnary == NArray.object(2):[gp0, gp1],
+ where the first object gp0 has the 1st coordinate [0,1,7,8] and
+ the second object gp1 has the 1st coordinate [3,4,5,6],
+ gpnary is restructured as [ gp0[0..1,false], gp1, gp0[2..3,false] ],
+ and join is made by using GPhys.join_md_nocheck.
+
+ This method is generally faster than GPhys.join unless the split
+ is one-dimensional.
+
+ ARGUMENT
+ * gpnarray [NArray of GPhys]
+
+ RETURN VALUE
+ * a GPhys
+
+---GPhys.join(gpary)
+ Join multiple GPhys objects (no need for any pre-ordering).
+
+ ARGUMENT
+ * gpnarray [Array (or 1D NArray) of GPhys]
+
+ RETURN VALUE
+ * a GPhys
+
+---GPhys.concat(gpary, axis_or_ary, name=nil, attr=nil)
+ Concatenate an Array (or 1D NArray) of GPhys objects
+ along the new dimension specified by the 2nd to 4th arguments.
+ The rank of the result (a GPhys) is one plus the rank of
+ the GPhys objects.
+
+ ARGUMENTS
+ * gpary [1D NArray or Array of GPhys]
+ * axis_or_ary [an Axis or a 1D NArray or Array of floats]
+ * name [String; optional] name of the coordinate;
+ needed if axis_or_ary is not an Axis.
+ * attr [Hash; optional] attributes of the coordinate;
+ used if axis_or_ary is not an Axis.
+
+ RETURN VALUE
+ * a GPhys
+
==Instance Methods
---data
Returns the data object
RETURN VALUE
@@ -85,13 +151,21 @@
RETURN VALUE
* a Grid
NOTE
- * There is a PROTECTED method (('grid')), which returns
- the grid object without duplicating.
+ * (('grid')) does not make a copy.
+---grid
+ Returns the grid object without copying.
+
+ RETURN VALUE
+ * a Grid
+
+ NOTE
+ * Use (('grid_copy')) to avoid side effects
+
---copy
Make a deep clone onto memory
RETURN VALUE
* a GPhys
@@ -434,21 +508,20 @@
* an Array of Integer
---shape
Aliased to ((<shape_current>))
----cyclic_ext(dim_or_dimname, modulo)
+---cyclic_ext(dim_or_dimname)
Extend a dimension cyclically.
The extension is done only when adding one grid point makes a full circle.
Thus, data at coordinate values [0,90,180,270] with modulo 360 are extended
(to at [0,90,180,270,360]), but data at [0,90,180] are not extended with
the same modulo: in this case, self is returned.
ARGUMENTS
* dim_or_dimname (String or Integer)
- * modulo (Numeric)
RETURN VALUE
* a GPhys (possibly self)
@@ -479,42 +552,43 @@
=end
require "numru/gphys/grid"
require "numru/misc/md_iterators"
require "numru/gphys/narray_ext"
+require "numru/gphys/mdstorage" # for GPhys.join
module NumRu
class GPhys
include NumRu::Misc::MD_Iterators
def initialize(grid, data)
- raise ArgumentError,"1st arg not a Grid" if ! grid.is_a?(Grid)
- raise ArgumentError,"2nd arg not a VArray" if ! data.is_a?(VArray)
+ #raise ArgumentError,"1st arg not a Grid" if ! grid.is_a?(Grid)
+ #raise ArgumentError,"2nd arg not a VArray" if ! data.is_a?(VArray)
if ( grid.shape_current != data.shape_current )
raise ArgumentError, "Shapes of grid and data do not agree. " +
"#{grid.shape_current.inspect} vs #{data.shape_current.inspect}"
end
@grid = grid
@data = data
end
attr_reader :grid, :data
- protected :grid
+ ###protected :grid # protection is lifted
def grid_copy
# deep clone of the grid
@grid.copy
end
def copy
# deep clone onto memory
- GPhys.new( @grid.copy, @data.copy )
+ self.class.new( @grid.copy, @data.copy )
end
def inspect
- "<GPhys grid=#{@grid.inspect}\n data=#{@data.inspect}>"
+ "<#{self.class} grid=#{@grid.inspect}\n data=#{@data.inspect}>"
end
def name
data.name
end
@@ -569,11 +643,11 @@
# ==NOTE:
# * VArray#convert_units does not copy data if to == @data.units
# * @grid is shared with self (no duplication)
# Thus, use GPhys#copy to separate all sub-objects (deep clone).
data = @data.convert_units(to)
- GPhys.new(@grid, data)
+ self.class.new(@grid, data)
end
def long_name
@data.long_name
end
@@ -586,11 +660,11 @@
slicer[0].keys[0].is_a?(String)
slicer = __process_hash_slicer(slicer[0])
else
slicer = __rubber_expansion( slicer )
end
- GPhys.new( @grid[*slicer], @data[*slicer] )
+ self.class.new( @grid[*slicer], @data[*slicer] )
end
def []=(*args)
val = args.pop
slicer = args
@@ -619,29 +693,29 @@
if has_assoccoord? && args.length==1 && ((spec=args[0]).is_a?(Hash)) &&
( acnms = (spec.keys & assoccoordnames ) ).length > 0
acspec = Hash.new
acnms.each{|nm| acspec[nm] = spec.delete(nm)}
grid, sl = @grid.cut_assoccoord(acspec)
- gphys = GPhys.new( grid, self.data[*sl] )
+ gphys = self.class.new( grid, self.data[*sl] )
else
gphys = self
end
newgrid, slicer = gphys.grid.cut( *args )
- GPhys.new( newgrid, gphys.data[ *slicer ] )
+ self.class.new( newgrid, gphys.data[ *slicer ] )
end
def cut_rank_conserving( *args )
newgrid, slicer = @grid.cut_rank_conserving( *args )
- GPhys.new( newgrid, @data[ *slicer ] )
+ self.class.new( newgrid, @data[ *slicer ] )
end
Axis.defined_operations.each do |method|
eval <<-EOS, nil, __FILE__, __LINE__+1
def #{method}(dim_or_dimname, *extra_args)
vary, grid = @grid.#{method}(@data, dim_or_dimname, *extra_args)
if grid
- GPhys.new( grid, vary )
+ self.class.new( grid, vary )
else
vary # scalar
end
end
EOS
@@ -723,25 +797,15 @@
end
## <-- For graphics
def coerce(other)
case other
- when Numeric
- ##na_other = self.data.val.fill(other) # Not efficient!
- va_other, = self.data.coerce(other)
- c_other = GPhys.new( @grid[ *([0..0]*self.rank) ],
- va_other.reshape!( *([1]*self.rank) ) )
- c_other.put_att('units',nil) # should be treated as such, not 1
- when Array, NArray
- va_other, = self.data.coerce(other)
- c_other = GPhys.new( @grid, va_other )
- c_other.put_att('units',nil) # should be treated as such, not 1
- when VArray
- c_other = GPhys.new( @grid, other )
- else
- raise "Cannot coerse #{other.class}"
- end
+ when Numeric, Array, NArrayMiss
+ c_other = UNumeric::Num2Coerce.new( other )
+ else
+ raise "Cannot coerse #{other.class}"
+ end
[c_other, self]
end
def shape_coerce(other)
#
@@ -774,31 +838,31 @@
end
def shape_coerce_full(other)
o, s = shape_coerce(other)
if o.length < s.length
- o = GPhys.new( s.grid, o + NArray.new(o.typecode,*s.shape) )
+ o = self.class.new( s.grid, o + NArray.new(o.typecode,*s.shape) )
elsif o.length > s.length
- s = GPhys.new( o.grid, s + NArray.new(s.typecode,*o.shape) )
+ s = self.class.new( o.grid, s + NArray.new(s.typecode,*o.shape) )
end
[o, s]
end
def transpose(*dims)
grid = @grid.transpose(*dims)
data = @data.transpose(*dims)
- GPhys.new( grid, data )
+ self.class.new( grid, data )
end
for f in VArray::Math_funcs
eval <<-EOS, nil, __FILE__, __LINE__+1
#def GPhys.#{f}(gphys)
# raise ArgumentError, "Not a GPhys" if !gphys.is_a?(GPhys)
- # GPhys.new( gphys.grid, VArray.#{f}(gphys.data) )
+ # self.class.new( gphys.grid, VArray.#{f}(gphys.data) )
#end
def #{f}(*arg)
- GPhys.new( self.grid, self.data.#{f}(*arg) )
+ self.class.new( self.grid, self.data.#{f}(*arg) )
end
EOS
end
for f in VArray::Binary_operators
eval <<-EOS, nil, __FILE__, __LINE__+1
@@ -811,22 +875,22 @@
newgrid = myself.grid.merge(other.grid)
else
vary = myself.data#{f} other
newgrid = myself.grid_copy
end
- GPhys.new( newgrid, vary )
+ self.class.new( newgrid, vary )
else
if other.respond_to?(:grid) #.is_a?(GPhys)
vary = myself#{f} other.data
else
vary = myself#{f} other
end
- GPhys.new( other.grid.copy, vary )
+ self.class.new( other.grid.copy, vary )
end
else
vary = self.data#{f} other
- GPhys.new( @grid.copy, vary )
+ self.class.new( @grid.copy, vary )
end
end
EOS
end
for f in VArray::Binary_operatorsL
@@ -839,18 +903,18 @@
end
for f in VArray::Unary_operators
eval <<-EOS, nil, __FILE__, __LINE__+1
def #{f}
vary = #{f.delete("@")} self.data
- GPhys.new( @grid.copy, vary )
+ self.class.new( @grid.copy, vary )
end
EOS
end
for f in VArray::NArray_type1_methods
eval <<-EOS, nil, __FILE__, __LINE__+1
def #{f}(*args)
- GPhys.new( self.grid.copy, self.data.#{f}(*args) )
+ self.class.new( self.grid.copy, self.data.#{f}(*args) )
end
EOS
end
for f in VArray::NArray_type2_methods
eval <<-EOS, nil, __FILE__, __LINE__+1
@@ -860,33 +924,35 @@
EOS
end
for f in VArray::NArray_type3_methods
eval <<-EOS, nil, __FILE__, __LINE__+1
def #{f}(*args)
+ arg_hash = args.pop if args[-1].is_a?(Hash)
+ # arg_hash: to support "min_count" in NArrayMiss
args = args.collect{|i| @grid.dim_index(i)}
- result = self.data.#{f}(*args)
+ args.push( arg_hash ) if arg_hash
+ result = self.data.#{f}(*args)
if Numeric===result || UNumeric===result
result
else
- GPhys.new( self.grid.delete_axes(args, "#{f}"), result )
+ args.pop if args[-1].is_a?(Hash)
+ self.class.new( self.grid.delete_axes(args, "#{f}"), result )
end
end
EOS
end
def shape_current
@data.shape_current
end
alias shape shape_current
- def cyclic_ext(dim_or_dimname, modulo)
+ # Old version of cyclic_ext
+ def cyclic_ext_with_modulo(dim_or_dimname, modulo)
# Cyclic extention to push the first element after the last element
# if appropriate.
- # <<developper's memo by horinout, 2005/01>>
- # in future modulo should be read based on conventions if nil
-
vx = coord(dim_or_dimname)
return self if vx.length <= 1
vvx = vx.val
width = (vvx[-1] - vvx[0]).abs
@@ -905,10 +971,31 @@
else
return self
end
end
+ def cyclic_ext(dim_or_dimname)
+ # Cyclic extention to push the first element after the last element
+ # if appropriate. (by using the cut method)
+
+ ax = axis(dim_or_dimname)
+ if ax.cyclic_extendible?
+ modulo = ax.modulo
+ v = ax.pos.val
+ v0 = v[0]
+ eps = 1e-2/ax.length
+ if v0 < v[-1] # increasing
+ v1 = v0 + modulo*(1+eps)
+ else
+ v1 = v0 - modulo*(1+eps)
+ end
+ return self.cut(ax.name=>v0..v1)
+ else
+ return self
+ end
+ end
+
def self.each_along_dims(gphyses, loopdims)
if !gphyses.is_a?(Array)
gphyses = [gphyses] # put in an Array (if a single GPhys)
end
gp = gphyses[0]
@@ -998,11 +1085,11 @@
vary = VArray.new(NArrayMiss.new(vtst.typecode, *grid.shape),
rs.data)
else
vary = VArray.new(NArray.new(vtst.typecode, *grid.shape), rs.data)
end
- results_whole.push( GPhys.new( grid, vary ) )
+ results_whole.push( self.new( grid, vary ) )
end
end
for j in 0...results.length
rs = results[j]
results_whole[j][idx_hash] = rs.data
@@ -1022,10 +1109,314 @@
def marshal_load(ary)
@data = ary[0]
@grid = ary[1]
end
+ ####### join multiple GPhys objects #######
+
+ def GPhys.join_md_nocheck(gpnary)
+ #< check >
+ if !gpnary.is_a?(NArray)
+ raise(ArgumentError,"Input must be an NArray of GPhys")
+ end
+ rank = gpnary.rank
+
+ #< axes >
+ gp0 = gpnary[0]
+
+ axes = Array.new
+ for d in 0...rank
+ if gpnary.shape[d] > 1 # --> join axes
+ sel = [0]*d + [true] + [0]*(rank-d-1) # [0,..0,true,0,...0]
+ axs = gpnary[*sel].collect{|gp| gp.axis(d)} # axes along d-th dim
+ ax = Axis.join(axs)
+ else
+ ax = gp0.axis(d)
+ end
+ axes.push(ax)
+ end
+
+ #< grid >
+ grid = Grid.new(*axes)
+ if gp0.assoc_coords
+ assoc_coords = gp0.assoccoordnames.collect do |aname|
+ GPhys.join( gpnary.collect{|gp| gp.assoc_coord_gphys(aname)}, true )
+ end
+ grid.set_assoc_coords(assoc_coords)
+ end
+
+ #< data >
+ data = VArrayComposite.new( gpnary.collect{|gp| gp.data} )
+
+ #< new gphys >
+ GPhys.new(grid, data)
+ end
+
+ def GPhys.join_md(gpnary)
+ #< Check >
+
+ if !gpnary.is_a?(NArray)
+ raise(ArgumentError,"Input must be an NArray of GPhys")
+ end
+ arank = gpnary.rank # rank of the input NArray
+ ashape = gpnary.shape
+ rank = gpnary[0].rank
+
+ #< Reshape and transpose gpnary if needed >
+
+ # / find dimmensions to join /
+ dimmap = Array.new
+ for i in 0...arank
+ if ashape[i] > 1 # join needed
+ sel = Array.new(arank, 0) # [0,0,...,0]
+ gp0 = gpnary[ *sel ]
+ sel[i] = 1
+ gp1 = gpnary[ *sel ] # [0,..,0,1,0,...,0]
+ for d in 0...rank
+ c0 = gp0.coord(d)[0]
+ c1 = gp1.coord(d)[0]
+ if c0.val != c1.convert_units(c0.units).val
+ dimmap[i] = d # dimension to join (found)
+ break
+ end
+ raise("Corresponding dim is not found for #{i}") if d==rank-1
+ end
+ else
+ dimmap[i] = nil # no need to join this dimension
+ end
+ end
+ if (x=dimmap-[nil]).length != x.uniq.length
+ raise "Dimensions to join cannot be determined uniquely"
+ end
+
+ # / "solo" dimensions (dimensions no need to join) /
+ sdims = (0...rank).collect{|d| d} - dimmap
+ for i in 0...arank
+ if dimmap[i].nil?
+ dimmap[i] = sdims.shift # assign dimensions orderly to
+ # minimize the need to transpose
+ end
+ end
+ sdims.each do |d|
+ dimmap.insert(d,d) # assign dimensions orderly to
+ gpnary = gpnary.newdim(d) # minimize the need to transpose
+ end
+ # now, gpnary.rank == rank
+
+ # / transpose gpnary if needed /
+ if dimmap != (0...rank).collect{|d| d}
+ imap = Array.new
+ dimmap.each_with_index do |d,j|
+ imap[j] = d
+ end
+ gpnary = gpnary.transpose(*imap)
+ end
+
+ #< Sort along dimensions to join >
+ gpnary = __sort_gpnary(gpnary)
+
+ #< Join! >
+ self.join_md_nocheck(gpnary)
+ end
+
+ # Join multiple GPhys objects (not need for any pre-ordering).
+ #
+ # ARGUMENT
+ # * gpnarray [Array (or 1D NArray) of GPhys]
+ #
+ def GPhys.join(gpary, ignore_overlap=false)
+
+ #< initialization with the first GPhys object >
+
+ gp = gpary[0]
+ rank = gp.rank
+ gpstore = MDStorage.new(rank)
+ gpstore[ *Array.new(rank, 0) ] = gp # first element
+ x0s = (0...rank).collect{|d|
+ pos = gp.axis(d).pos
+ x0 = UNumeric[ pos.val[0], pos.units ]
+ [ x0 ] # first values of each coordinate
+ }
+
+ #< scan the coordiantes of the remaining GPhys objects >
+ for k in 1...gpary.length
+ gp = gpary[k]
+ idx = Array.new
+ for d in 0...rank
+ pos = gp.axis(d).pos
+ x0 = UNumeric[ pos.val[0], pos.units ]
+ i = x0s[d].index(x0)
+ if i.nil?
+ x0s[d].push(x0)
+ i = x0s[d].length-1
+ end
+ idx.push(i)
+ end
+ gpstore[*idx] = gp
+ end
+
+ if !ignore_overlap && gpstore.count_non_nil != gpary.length
+ raise(ArgumentError,"Cannot uniquely locate one or more objects; some overlap in the grids?")
+ end
+
+ gpnary = gpstore.to_na
+
+ #< Sort along dimensions to join >
+ gpnary = __sort_gpnary(gpnary)
+
+ #< Join! >
+ self.join_md_nocheck(gpnary)
+ end
+
+ def GPhys.concat(gpary, axis_or_ary, name=nil, attr=nil)
+ #< interpret gpary (1st arg) >
+ gpary = NArray.to_na(gpary) if gpary.is_a?(Array)
+ if !gpary.is_a?(NArray) || gpary.rank != 1
+ raise(ArgumentError,"1st arg must be a 1D NArray or Array of GPhys")
+ end
+ len = gpary.length
+
+ #< interpret axis_or_ary (2nd arg) and make an Axis if not >
+ case axis_or_ary
+ when Axis
+ ax = axis_or_ary
+ if ax.length != len
+ raise(ArgumentError,"length mismatch #{len} vs #{ax.length}")
+ end
+ else
+ ary = axis_or_ary # must be an NArray or Array
+ ary = NArray.to_na(ary).to_f if ary.is_a?(Array)
+ if !ary.is_a?(NArray) || ary.rank != 1
+ raise(ArgumentError,
+ "If not an Axis, 2nd arg must be 1D NArray or Array of float")
+ end
+ if ary.length != len
+ raise(ArgumentError,"length mismatch #{len} vs #{ary.length}")
+ end
+ if name.nil?
+ raise(ArgumentError,
+ "3rd arg (name) is needed if the 2nd arg is not an Axis")
+ end
+ va = VArray.new(ary, attr, name)
+ ax = Axis.new().set_pos(va)
+ end
+
+ #< new grid >
+ grid = gpary[0].grid.insert_axis(-1,ax) # insert_axis: non-destructive
+
+ #< join VArrays >
+ ds = gpary.collect{|gp| gp.data}
+ gpary[0].rank.times{ds.newdim!(0)} # for VArrayComposite.new
+ data = VArrayComposite.new(ds)
+
+ #< result >
+ GPhys.new(grid, data)
+ end
+
+ ############## < private class methods > ##############
+ def GPhys::__sort_gpnary(gpnary)
+ #< Sort along dimensions to join >
+
+ shape = gpnary.shape
+ rank = shape.length
+ for d in 0...rank
+ n = shape[d]
+ if n > 1 # --> d is a dimesnion to join; possibly need to sort
+ sel = Array.new(rank, 0) # [0,0,...,0]
+ xs = Array.new # will be [ [pos, k-th obj, j-th elem], ... ]
+ xunits = gpnary[0].axis(d).pos.units
+ for k in 0...n
+ sel[d] = k
+ axis = gpnary[*sel].axis(d)
+ if axis.pos.units != xunits
+ pos = axis.pos.convert_units(xunits)
+ if !axis.cell?
+ axis.set_pos(pos)
+ else
+ if axis.cell_center?
+ bds = axis.cell_bounds.convert_units(xunits)
+ axis.set_cell(pos,bds)
+ else
+ ctrs = axis.cell_center.convert_units(xunits)
+ axis.set_cell(ctr,pos)
+ end
+ end
+ end
+ x = axis.pos.val
+ for i in 0...x.length
+ xs.push( {:pos=>x[i], :kobj=>k, :idx=>i} )
+ end
+ end
+
+ if ( gpnary[0].axis(d).length == 1 or # -> cannot determine order
+ xs[0][:pos] <= xs[1][:pos] ) # -> increasing
+ xs.sort!{|a,b| a[:pos] <=> b[:pos]} # increasing order
+ else
+ xs.sort!{|a,b| b[:pos] <=> a[:pos]} # decreasing order
+ end
+
+ x0 = xs.shift
+ xconf = [ { :kobj=>x0[:kobj], :ids=>[x0[:idx]] } ]
+ ic = 0
+ xs.each do |x|
+ if x[:kobj] == xconf[ic][:kobj] # same obj
+ xconf[ic][:ids].push(x[:idx])
+ else # -> new obj
+ ic += 1
+ xconf[ic] = {:kobj=>x[:kobj], :ids=>[x[:idx]] }
+ end
+ end
+ if xconf.length == n
+ # At most, only reodering is needed
+ gpnary2 = gpnary.clone # to replace gpnary
+ selj = Array.new(rank, true) # [true, true,...]
+ selk = Array.new(rank, true) # [true, true,...]
+ for j in 0...xconf.length
+ if j != (k=xconf[j][:kobj]) # need to reorder the dimension d
+ selj[d] = j
+ selk[d] = k
+ gpnary2[*selj] = gpnary[*selk]
+ end
+ end
+ gpnary = gpnary2
+ else
+ # Need to divide some object(s)
+ shape2 = shape.clone
+ n2 = xconf.length
+ shape2[d] = n2
+ gpnary2 = NArray.object(*shape2) # to replace gpnary
+ selj = Array.new(rank, true) # [true, true,...]
+ selk = Array.new(rank, true) # [true, true,...]
+ for j in 0...xconf.length
+ selj[d] = j
+ selk[d] = xconf[j][:kobj]
+ ids = xconf[j][:ids] # indices to subset d-th D of gpnary elem
+ regular = true
+ for b in 0...ids.length-1
+ if ids[b+1] - ids[b+1] != 1
+ regular = false
+ break
+ end
+ end
+ if regular
+ ids = ids[0]..ids[-1]
+ else
+ #ids = NArray.to_na(ids)
+ end
+ seld = Array.new(rank, true) # [true, true,...]
+ seld[d] = ids
+ gpnary2[*selj] = gpnary[*selk].collect{|g| g[*seld]}
+ end
+ gpnary = gpnary2
+ end
+ end
+
+ end
+ gpnary
+ end
+ private_class_method :__sort_gpnary
+
############## < private methods > ##############
private
def __rubber_expansion( args )
if (id = args.index(false)) # substitution into id
@@ -1180,12 +1571,14 @@
p(sp=gpz3.shape)
gpz3_cext = gpz3[ [0...sp[0],0], false ]
p gpz3_cext.coord(0).val, gpz3_cext.val
print "cyclic extention if appropriate\n"
- gpz3_cext = gpz3.cyclic_ext(0,10.0)
+ gpz3.axis(0).pos.set_att("modulo",[10.0])
+ gpz3_cext = gpz3.cyclic_ext(0)
p gpz3_cext.coord(0).val, gpz3_cext.val
+ gpz3.axis(0).pos.del_att("modulo")
print "axis to gphys\n"
ax = gpz3.axis(1)
p ax.to_gphys
p ax.to_gphys( ax.cell_bounds[0..-2].copy )
@@ -1218,6 +1611,104 @@
mar = Marshal.dump(gpz)
p mar.class, mar.length
g = Marshal.load(mar)
p g
+ # test join
+ p "## testing join_md_nocheck..."
+ nx = 4
+ ny = 3
+ nz = 2
+ vx0 = VArray.new( NArray.float(nx).indgen! ).rename("x")
+ vy0 = VArray.new( NArray.float(ny).indgen! + 0.5 ).rename("y")
+ vyb0 = VArray.new( NArray.float(ny+1).indgen! ).rename("yb")
+ vx1 = VArray.new( NArray.float(nx).indgen! + nx ).rename("x")
+ vy1 = VArray.new( NArray.float(ny).indgen! + (ny+0.5) ).rename("y")
+ vyb1 = VArray.new( NArray.float(ny+1).indgen! + ny ).rename("yb")
+ xax0 = Axis.new().set_pos(vx0)
+ yax0 = Axis.new(true).set_cell(vy0,vyb0).set_pos_to_center
+ xax1 = Axis.new().set_pos(vx1)
+ yax1 = Axis.new(true).set_cell(vy1,vyb1).set_pos_to_center
+ vz = VArray.new( NArray.float(nz).indgen! ).rename("z")
+ zax = Axis.new().set_pos(vz)
+ v = VArray.new( NArray.float(nx,ny,nz).indgen!, nil, "V" )
+ gp00 = GPhys.new( Grid.new(xax0, yax0, zax), v )
+ gp10 = GPhys.new( Grid.new(xax1, yax0, zax), v+100 )
+ gp01 = GPhys.new( Grid.new(xax0, yax1, zax), v+200 )
+ gp11 = GPhys.new( Grid.new(xax1, yax1, zax), v+300 )
+ x = NArray.float(nx,ny).random!
+=begin
+ gp00.set_assoc_coords([GPhys.new(Grid.new(xax0,yax0),VArray.new(x,nil,"A"))])
+ gp10.set_assoc_coords([GPhys.new(Grid.new(xax1,yax0),VArray.new(x,nil,"A"))])
+ gp01.set_assoc_coords([GPhys.new(Grid.new(xax0,yax1),VArray.new(x,nil,"A"))])
+ gp11.set_assoc_coords([GPhys.new(Grid.new(xax1,yax1),VArray.new(x,nil,"A"))])
+=end
+ gpnary = NArray.to_na( [[ [gp00,gp10], [gp01,gp11] ]] )
+ gpu = GPhys.join_md_nocheck(gpnary)
+ p "#",gpu, gpu.coord(0).val, gpu.coord(1).val, gpu.axis(1).cell_bounds.val
+ p "#",gpu.val
+
+ p "## testing join_md..."
+ gpnary = NArray.to_na( [ [gp00,gp01], [gp10,gp11] ] )
+ gpu = GPhys.join_md(gpnary)
+ p gpu.copy, gpu.val
+
+ gpnary = NArray.to_na( [gp00,gp01] )
+ gpu = GPhys.join_md(gpnary)
+ p gpu.copy, gpu.val
+
+ gpnary = NArray.to_na( [gp10,gp00] )
+ gpu = GPhys.join_md(gpnary)
+ p gpu.copy, gpu.val
+
+ gp00.coord(0)[0..1] = NArray.to_na([1, 2])
+ gp00.coord(0)[2..3] = NArray.to_na([7, 8])
+ gp10.coord(0)[0..3] = NArray.to_na([3,4,5,6])
+ gpnary = NArray.to_na( [gp00,gp10] )
+ gpu = GPhys.join_md(gpnary)
+ p gpu.copy, gpu.val
+
+ ans = NArray.to_na(
+ [ [ [ 0.0, 1.0, 100.0, 101.0, 102.0, 103.0, 2.0, 3.0 ],
+ [ 4.0, 5.0, 104.0, 105.0, 106.0, 107.0, 6.0, 7.0 ],
+ [ 8.0, 9.0, 108.0, 109.0, 110.0, 111.0, 10.0, 11.0 ] ],
+ [ [ 12.0, 13.0, 112.0, 113.0, 114.0, 115.0, 14.0, 15.0 ],
+ [ 16.0, 17.0, 116.0, 117.0, 118.0, 119.0, 18.0, 19.0 ],
+ [ 20.0, 21.0, 120.0, 121.0, 122.0, 123.0, 22.0, 23.0 ] ] ] )
+ if gpu.val == ans
+ p "test ok"
+ else
+ raise "test failed"
+ end
+
+ p "## testing join..."
+ gpu = GPhys.join( [gp11,gp10,gp01,gp00] )
+ p gpu.copy, gpu.val
+ if gpu.val[true,0..2,true] == ans and gpu.val[true,3..5,true] == ans+200
+ p "test ok"
+ else
+ raise "test failed"
+ end
+
+ p "## testing to confirm the failure of join..."
+ begin
+ gpu = GPhys.join( [gp11,gp11+0.1,gp01,gp00] )
+ # overlap in grids: cannot not uniquely join
+ rescue ArgumentError
+ puts "test ok: confirmed the expected failure."
+ end
+
+ p "## testing concat..."
+ gpary = [gp00[0..2,0..1,0], gp00[0..2,0..1,0]+100]
+ gpu = GPhys.concat( gpary, [0.0, 1.0], "time", {"units"=>"days"} )
+ ans = NArray.to_na( [ [ [ 0.0, 1.0, 2.0 ],
+ [ 4.0, 5.0, 6.0 ] ],
+ [ [ 100.0, 101.0, 102.0 ],
+ [ 104.0, 105.0, 106.0 ] ] ] )
+ if gpu.coord(2).name == "time" && gpu.val == ans
+ p "test ok"
+ else
+ raise "test failed"
+ end
+
end
+