=begin
= NArrayMiss Class

NArrayMiss is a additional class processing of missing value with to
((<NArray|URL:http://www.ruby-lang.org/en/raa-list.rhtml?name=NArray>))
for Ruby.

To use NArrayMiss class, you need invoking "require 'narray_miss.rb'" in your script.

== Index

* ((<Class Constants>))
* ((<Class Methods>))
* ((<Class Instance Methods>))
  * ((<NArrayMiss information>))
  * ((<Slicing Array>))
  * ((<Filling values>))
  * ((<Arithmetic operator>))
  * ((<Bitwise operator|Bitwise operator (only for byte, sint and int)>))
  * ((<Comparison>))
  * ((<Statistics>))
  * ((<Sort>))
  * ((<Transpose>))
  * ((<Changing Shapes of indices>))
  * ((<Type conversion>))
  * ((<Iteration>))
  * ((<Boolean and mask related|Boolean and mask related (only for byte, sint and int)>))
  * ((<Complex compound number|Complex compound number (only for scomplex and complex)>))
  * ((<Byte swap>))
  * ((<Mask and missing value>))
  * ((<Others>))


=end

require 'narray'


class NArrayMiss

=begin
== Class Constants
--- NArrayMiss::BYTE
     type code for 1 byte unsigned integer.
--- NArrayMiss::SINT
     type code for 2 byte signed integer.
--- NArrayMiss::INT
     type code for 4 byte signed integer.
--- NArrayMiss::SFLOAT
     type code for single precision float.
--- NArrayMiss::FLOAT
     type code for double precision float.
--- NArrayMiss::SCOMPLEX
     type code for single precision complex.
--- NArrayMiss::COMPLEX
     type code for double precision complex.
--- NArrayMiss::OBJECT
     type code for Ruby object.

go back to ((<Index>))
=end

  BYTE = NArray::BYTE
  SINT = NArray::SINT
  INT = NArray::INT
  SFLOAT = NArray::SFLOAT
  FLOAT = NArray::FLOAT
  SCOMPLEX = NArray::SCOMPLEX
  COMPLEX = NArray::COMPLEX
  OBJECT = NArray::OBJECT

  class << self
    alias :__new__ :new
    private :__new__
  end

  def initialize(array, mask)
    if array.shape!=mask.shape
      raise "array and mask must have the same shape"
    end
    @array = array
    @mask = mask
  end
  private :initialize

=begin
== Class Methods
--- NArrayMiss.new(typecode, size, ...)
     create (({NArrayMiss})) of ((|typecode|)).
     All elements are initialized with 0.
--- NArrayMiss.byte(size, ...)
     same as NArrayMiss.new(NArrayMiss::BYTE, ((|size|)), ...).
--- NArrayMiss.sint(size, ...)
     same as NArrayMiss.new(NArrayMiss::SINT, ((|size|)), ...).
--- NArrayMiss.int(size, ...)
     same as NArrayMiss.new(NArrayMiss::INT, ((|size|)), ...).
--- NArrayMiss.sfloat(size, ...)
     same as NArrayMiss.new(NArrayMiss::SFLOAT, ((|size|)), ...).
--- NArrayMiss.float(size, ...)
     same as NArrayMiss.new(NArrayMiss::FLOAT, ((|size|)), ...).
--- NArrayMiss.scomplex(size, ...)
     same as NArrayMiss.new(NArrayMiss::SCOMPLEX, ((|size|)), ...).
--- NArrayMiss.complex(size, ...)
     same as NArrayMiss.new(NArrayMiss::COMPLEX, ((|size|)), ...).
--- NArrayMiss.object(size, ...)
     same as NArrayMiss.new(NArrayMiss::OBJECT, ((|size|)), ...).
--- NArrayMiss[](value, ...)
     create (({NArrayMiss})) form [((|value|)), ...].
--- NArrayMiss.to_nam(array [,mask])
     create (({NArrayMiss})) from ((|array|)).
     ((|array|)) must be (({Array})) or (({NArray})).
--- NArrayMiss.to_nam_no_dup(array [,mask])
     convert from ((|array|)) to (({NArrayMiss})).

go back to ((<Index>))
=end

  def self.new(*arg)
    array = NArray.new(*arg)
    mask = NArray.byte(*arg[1..-1])
    __new__(array, mask)
  end
  def self.byte(*arg)
    NArrayMiss.new(BYTE,*arg)
  end
  def self.sint(*arg)
    NArrayMiss.new(SINT,*arg)
  end
  def self.int(*arg)
    NArrayMiss.new(INT,*arg)
  end
  def self.sfloat(*arg)
    NArrayMiss.new(SFLOAT,*arg)
  end
  def self.float(*arg)
    NArrayMiss.new(FLOAT,*arg)
  end
  def self.scomplex(*arg)
    NArrayMiss.new(SCOMPLEX,*arg)
  end
  def self.complex(*arg)
    NArrayMiss.new(COMPLEX,*arg)
  end
  def self.object(*arg)
    NArrayMiss.new(OBJECT,*arg)
  end
  def self.[](*arg)
    NArrayMiss.to_nam(NArray[*arg])
  end
  def self.to_nam_no_dup(*arg)
    if arg.length > 2 || arg.length==0 then
      raise("NArrayMiss.to_nar( array [,mask]] )")
    end

    array = arg[0]
    if Numeric===array then array = NArray[array] end
    if Array===array then array = NArray.to_na(array) end
    if !array.is_a?(NArray) then
      raise("argument must be Numeric, NArray or Array")
    end

    if arg.length==2 then
      mask = arg[1]
      if Numeric===mask then mask = [mask] end
      if Array===mask then
	mask = NArray.to_na(mask).ne(0)
      end
      if mask.class == FalseClass then
	mask = NArray.byte(*array.shape)
      end
      if mask.class == TrueClass then
	mask = NArray.byte(*array.shape).fill(1)
      end
      if !(NArray===mask && mask.typecode==BYTE) then
	  raise("mask must be Numeric, Array, true, false or NArray(byte)")
      end
      if mask.length!=array.length
        raise "mask.length must be same as array.length"
      end
    else
      mask = NArray.byte(*array.shape).fill(1)
    end
    __new__(array,mask)
  end
  def self.to_nam(*arg)
    if !(Numeric===arg[0]) && Array===arg[0] && !arg[0].is_a?(NArray)
      raise "first argument must be Numeric, NArray or Array"
    end
    arg[0] = arg[0].dup if !(Numeric===arg[0])
    if arg.length==2 && !(Numeric===arg[1]) && arg[1].class!=TrueClass && arg[1].class!=FalseClass then
      arg[1] = arg[1].dup
    end
    NArrayMiss.to_nam_no_dup(*arg)
  end


=begin
== Class Instance Methods
=end


=begin
=== NArrayMiss information
--- NArrayMiss#dim
     return the dimension which is the number of indices.
--- NArrayMiss#rank
     same as (({NArrayMiss#dim})).
--- NArrayMiss#shape
     return the (({Array})) of sizes of each index.
--- NArrayMiss#size
     return the number of total elements.
--- NArrayMiss#total
     alias to size
--- NArrayMiss#length
     alias to size
--- NArrayMiss#rank_total
     return the number of total of the shape.
--- NArrayMiss#typecode
     return the typecode.
=end

  def dim
    @array.dim
  end
  def rank
    @array.rank
  end
  def shape
    @array.shape
  end
  def size
    @array.size
  end
  alias :total :size
  alias :length :size

  def rank_total(*arg)
    @array.rank_total(*arg)
  end

  def typecode
    @array.typecode
  end


=begin
=== Slicing Array
--- NArrayMiss#[](index)
     return the value at [((|index|))].
     ((|index|)) must be (({Integer, Range, Array, true})).
     Index order is FORTRAN type.
--- NArrayMiss#slice(index)
     same as (({NArrayMiss#[]})) but keeps the rank of original array by not elimiting dimensions whose length became equal to 1 (which (({NArrayMiss#[]})) dose).
     This is not the case with the 1-dimensional indexing and masking.
--- NArrayMiss#set_without_validation(index,value)
     replace elements at ((|index|)) by ((|value|)).
--- NArrayMiss#[]=(index, value)
     replace elements at ((|index|)) by ((|value|)) and
     make replaced elements valid.     
=end

  def [](*arg)
    obj = @array[*arg]
    if Numeric===obj
      return obj
    else
      return NArrayMiss.to_nam_no_dup(obj,@mask[*arg])
    end
  end
  def slice(*arg)
    NArrayMiss.to_nam_no_dup(@array.slice(*arg),@mask.slice(*arg))
  end

  def set_without_validation(*arg)
    if arg.length==1 then
      if !arg[0] then
	@mask[] = 0
      elsif arg[0].class == NArrayMiss then
	@array[] = arg[0].get_array!
	@mask[] = arg[0].get_mask!
      else
	@array[] = arg[0]
      end
    else
      if !arg[-1] then
	@mask[*arg[0..-2]] = 0
      elsif arg[-1].class == NArrayMiss then
	@array[*arg[0..-2]] = arg[-1].get_array!
	@mask[*arg[0..-2]] = arg[-1].get_mask!
      else
	@array[*arg[0..-2]] = arg[-1]
      end
    end
    return self
  end

  def []=(*arg)
    self.set_without_validation(*arg)
    if arg[-1].class != NArrayMiss && arg[-1] then
      if arg.length==1 then
        @mask=1
      else
        @mask[*arg[0..-2]] = 1
      end
    end
    return self
  end


=begin
=== Filling values
--- NArrayMiss#indgen!([start[,step]])
     set values from ((|start|)) with ((|step|)) increment.
--- NArrayMiss#indgen([start[,step]])
     same as (({NArrayMiss#indgen!})) but create new object.
--- NArrayMiss#fill!(value)
     fill elements with ((|value|)).
--- NArrayMiss#fill(value)
     same as (({NArrayMiss#fill!})) but create new object.
--- NArrayMiss#random!(max)
     set random values between 0<=x<((|max|)).
--- NArrayMiss#random(max)
     same as (({NArrayMiss#random!})) but create new object.
--- NArrayMiss#randomn(max)
     set normally distributed random values with mean=0, dispersion=1 (Box-Muller)
=end

  for operator in ["indgen","fill","random"]
    eval(<<-EOL,nil,__FILE__,__LINE__+1)
    def #{operator}(*arg)
      obj = self.dup
      obj.#{operator}!(*arg)
      obj
    end
    def #{operator}!(*arg)
      @array = @array.#{operator}!(*arg)
      @mask[true] = 1
      self
    end
    EOL
  end
  def randomn
    obj = self.dup
    obj[@mask] = @array[@mask].randomn
    obj
  end


=begin
=== Arithmetic operator
--- NArrayMiss#-@
--- NArrayMiss#+(other)
--- NArrayMiss#-(other)
--- NArrayMiss#*(other)
--- NArrayMiss#/(other)
--- NArrayMiss#%(other)
--- NArrayMiss#**(other)
--- NArrayMiss#abs
--- NArrayMiss#add!(other)
--- NArrayMiss#sbt!(other)
--- NArrayMiss#mul!(other)
--- NArrayMiss#div!(other)
--- NArrayMiss#mod!(other)
--- NArrayMiss#mul_add(other, dim, ...)
=end

  def -@
    array = @array.dup
    array[@mask] = -@array[@mask]
    NArrayMiss.to_nam_no_dup(array, @mask.dup)
  end
  for operator in ["+","-","*","/","%","**"]
    dummy = {"+"=>0,"-"=>0,"*"=>1,"/"=>1,"%"=>1,"**"=>1}[operator]
    eval(<<-EOL,nil,__FILE__,__LINE__+1)
    def #{operator}(arg)
      if !arg then
	@mask = 0
	return self
      else
        term1,term2,mask,flag = routine1(arg,#{dummy})
        result = term1 #{operator} term2
        routine2(result,mask,flag)
      end
    end
    EOL
  end
  def abs
    array = @array.dup
    array[@mask] = @array[@mask].abs
    NArrayMiss.to_nam_no_dup(array, @mask.dup)
  end

  for operator in ["add!","sbt!","mul!","div!","mod!"]
    dummy = {"add!"=>0,"sbt!"=>0,"mul!"=>1,"div!"=>1,"mod!"=>1}[operator]
    eval(<<-EOL,nil,__FILE__,__LINE__+1)
    def #{operator}(arg)
      term1,term2,mask,flag = routine1(arg,#{dummy})
      result = term1.#{operator} term2
      routine2(result,mask,flag)
    end
    EOL
  end
  def mul_add(*arg)
    if arg.length==1 then
      return (self*arg[0]).sum
    else
      return (self*arg[0]).sum(*arg[1..-1])
    end
  end


=begin
=== Bitwise operator (only for byte, sint and int)
--- NArrayMiss#~@
--- NArrayMiss#&(other)
--- NArrayMiss#|(other)
--- NArrayMiss#^(other)
=end

  def ~@
    NArrayMiss.to_nam_to_dup(~@array, @mask.dup)
  end
  for operator in ["&","|","^"]
    dummy = {"&"=>1,"|"=>0,"^"=>1}[operator]
    eval(<<-EOL,nil,__FILE__,__LINE__+1)
    def #{operator}(arg)
      term1,term2,mask,flag = routine1(arg,#{dummy})
      result = term1 #{operator} term2
      routine2(result,mask,flag)
    end
    EOL
  end


=begin
=== Comparison
--- NArrayMiss#eq(other)
--- NArrayMiss#ne(other)
--- NArrayMiss#gt(other)
--- NArrayMiss#ge(other)
--- NArrayMiss#lt(other)
--- NArrayMiss#le(other)
--- NArrayMiss#>(other)
--- NArrayMiss#>=(other)
--- NArrayMiss#<(other)
--- NArrayMiss#<=(other)
--- NArrayMiss#and(other)
--- NArrayMiss#or(other)
--- NArrayMiss#xor(other)
--- NArrayMiss#not(other)
=end

  for operator in ["eq","ne","gt","ge","lt","le"]
    eval(<<-EOL,nil,__FILE__,__LINE__+1)
    def #{operator}(arg)
      term1,term2,mask,flag = routine1(arg,0)
      result = term1.#{operator} term2
      routine2(result,mask,flag)
    end
    EOL
  end
  for operator in [">",">=","<","<="]
    eval(<<-EOL,nil,__FILE__,__LINE__+1)
    def #{operator}(arg)
      term1,term2,mask,flag = routine1(arg,0)
      result = term1 #{operator} term2
      routine2(result,mask,flag)
    end
    EOL
  end

  for operator in ["and","or","xor"]
    dummy = {"and"=>1,"or"=>0,"xor"=>1}[operator]
    eval(<<-EOL,nil,__FILE__,__LINE__+1)
    def #{operator}(arg)
      term1,term2,mask,flag = routine1(arg,#{dummy})
      result = term1.#{operator} term2
      routine2(result,mask,flag)
    end
    EOL
  end
  def not
    NArrayMiss.to_nam_no_dup(@array.not, @mask.dup)
  end

#  def ==(arg)
#    if art.kind_of?(NArrayMiss) then
#      @array==arg.get_array! && @mask==arg.get_mask!
#    else
#      false
#    end
#  end


=begin
=== Statistics
--- NArrayMiss#sum(dim, ... ["min_count"=>i])
     return summation of elements in specified dimensions.
     Elements at which the number of elements for summation is less than ((|i|)) is invalid.
--- NArrayMiss#accum(dim, ...)
     same as (({NArrayMiss#sum})) but not elimiting dimensions whose length became equal to 1.
--- NArrayMiss#min(dim, ...)
     return minimum in specified dimensions.
     Elements at which the number of valid elements in the dimension is less than ((|i|)) is invalid.
--- NArrayMiss#max(dim, ...)
     return maximum in specified dimensions.
     Elements at which the number of valid elements in the dimension is less than ((|i|)) is invalid.
--- NArrayMiss#median(dim, ...)
     return median in specified dimensions.
     Elements at which the number of valid elements in the dimension is less than ((|i|)) is invalid.
--- NArrayMiss#mean(dim, ...)
     return mean of elements in specified dimensions.
     Elements at which the number of elements for then mean is less than ((|i|)) is invalid.
--- NArrayMiss#stddev(dim, ...)
     return standard deviation of elements in specified dimensions.
     Elements at which the number of elements for calculation the standard deviation is less than ((|i|)) is invalid.
=end

  def accum(*arg)
    if @mask.count_true == 0 then
      return nil
    else
      array = @array.dup
      array[@mask.not] = 0
      return NArrayMiss.to_nam_no_dup(array.accum(*arg),
			       @mask.to_type(NArray::INT).accum(*arg).ne(0))
    end
  end

  for operator in ["sum","min","max"]
    str = {"sum"=>"0","min"=>"array.max","max"=>"array.min",}[operator]
    eval(<<-EOL,nil,__FILE__,__LINE__+1)
    def #{operator}(*arg)
      if @mask.count_true == 0 then
        return nil
      end
      min_count=1
      options = ["min_count"]
      if arg.length!=0 && arg[-1].class==Hash then
	option = arg[-1]
	arg = arg[0..-2]
	option.each_key{|key|
	  if !options.index(key) then
	    raise(ArgumentError,key+" option is not exist")
	  end
	}
	min_count = option["min_count"]
      end
      mask = @mask.to_type(NArray::INT)
      array = @array.dup
      array[@mask.not] = #{str}
      mask = mask.sum(*arg)
      if NArray===mask then
	if mask.ge(min_count).count_true == 0 then
	  return nil
        else
          return NArrayMiss.to_nam_no_dup(array.#{operator}(*arg),
                                          mask.ge(min_count))
        end
      else
        if mask < min_count
          return nil
        else
          return array.#{operator}
        end
      end
    end
    EOL
  end
  def mean(*arg)
    if @mask.count_true == 0 then
      return nil
    end
    min_count = 1
    options=["min_count"]
    if arg.length!=0 && arg[-1].class==Hash then
      option = arg[-1]
      arg2=arg[0..-2]
      option.each_key{|key|
	if !options.index(key) then
	  raise(ArgumentError,key+" option is not exist")
	end
      }
      min_count = option["min_count"]
    else
      arg2=arg
    end
    count = @mask.to_type(NArray::INT).sum(*arg2)
    if count.class == NArray then
      count = NArrayMiss.to_nam(count,count.ge(min_count))
      if count.ge(min_count).count_true == 0
        return nil
      else
        return self.sum(*arg2)/count
      end
    else
      if count < min_count then
        return nil
      else
	return self.sum(*arg2)/count
      end
    end
  end
  def stddev(*arg)
    if @mask.count_true == 0 then
      return nil
    end
    min_count=2
    options=["min_count"]
    if arg.length!=0 && arg[-1].class==Hash then
      option = arg[-1]
      arg = arg[0..-2]
      option.each_key{|key|
	if !options.index(key) then
	  raise(ArgumentError,key+" option is not exist")
	end
      }
      min_count = option["min_count"].to_i
      if min_count<2
        raise(ArgumentError, "min_count must be >= 2")
      end
    end
    count = @mask.to_type(NArray::INT).sum(*arg)
    count2 = @mask.to_type(NArray::INT).accum(*arg)
    if count.class==NArray then
      count = NArrayMiss.to_nam_no_dup(count,count.ge(min_count))
      if count.get_mask!.count_true == 0
        return nil
      end
    else
      if count < min_count then
	return nil
      end
    end
    if self.integer? then
      a = self.to_f
    else
      a = self
    end
    var = ( (a-a.accum(*arg)/count2)**2 ).sum(*arg)/(count-1)
    obj = NMMath::sqrt(var)
    return obj
  end

  def median(*arg)
    if arg.length==0 then
      return @array[@mask].median
    else
      nshape = NArray.to_na(@array.shape)
      nshape[arg]=1
      nslice = nshape[nshape.ne(1).where]
      index = NArray.object(@mask.rank)
      index[nshape.eq(1).where] = true
      obj = NArrayMiss.new(@array.typecode,*nslice.to_a)
      total = 1
      nslice.each{|n| total *= n}
      for i in 0...total
        index[nshape.ne(1).where] = pos(i,nslice)
        mask = NArray.byte(*@array.shape).fill(0)
        mask[*index] = 1
        mask = @mask&mask
        if mask.count_true != 0 then
          obj[*pos(i,nslice)] = @array[mask].median
        end
      end
      return obj
    end
  end


=begin
=== Sort
--- NArrayMiss#sort(dim)
     sort in the 0..((|dim|)) (All dimensions if omitted)
--- NArrayMiss#sort_index(dim)
     return index of sort result.
=end

  for operator in ["sort","sort_index"]
    eval(<<-EOL,nil,__FILE__,__LINE__+1)
    def #{operator}(*arg)
      obj=NArrayMiss.new(@array.typecode,*@array.shape)
      if arg.length==0 then
	obj[@mask] = @array[@mask].#{operator}
        return obj
      else
        nshape = NArray.to_na(@array.shape)
        nshape[arg]=1
        nslice = nshape[nshape.ne(1).where]
        index = NArray.object(@mask.rank)
        index[nshape.eq(1).where] = true
        obj = NArrayMiss.new(@array.typecode,*@array.shape)
        total = 1
        nslice.each{|n| total *= n}
        for i in 0...total
          index[nshape.ne(1).where] = pos(i,nslice)
          mask = NArray.byte(*@array.shape).fill(0)
          mask[*index] = 1
          mask = @mask&mask
          if mask.count_true != 0 then
	    obj[mask] = @array[mask].#{operator}
	  end
	end
	return obj
      end
    end
    EOL
  end


=begin
=== Transpose
--- NArrayMiss#transpose(dim0, dim1, ...)
     transpose array. The 0-th dimension goes to the ((|dim0|))-th dimension of new array.
=end

  def transpose(*arg)
    obj = self.dup
    shape = arg.collect{|i| obj.shape[i]}
    obj.reshape!(*shape)
    obj.set_without_validation( @array.transpose(*arg) )
    obj.set_mask(@mask.transpose(*arg))
    obj
  end


=begin
=== Changing Shapes of indices
--- NArrayMiss#reshape!(size, ...)
     change shape of array.
--- NArrayMiss#reshape(size, ...)
     same as (({NArrayMiss#reshape!})) but create new object.
--- NArrayMiss#shape=(size, ...)
     same as (({NArrayMiss#reshape!})).
--- NArrayMiss#newdim!(dim)
     insert new dimension with size=1
--- NArrayMiss#newdim(dim)
     same as (({NArrayMiss#newdim!})) but create new object.
--- NArrayMiss#rewrank!(dim)
     same as (({NArrayMiss#newdim!})).
--- NArrayMiss#rewrank=(dim)
     same as (({NArrayMiss#newdim!})).
=end

  def reshape!(*arg)
    @array = @array.reshape!(*arg)
    @mask = @mask.reshape!(*arg)
    self
  end
  def reshape(*arg)
    obj = self.dup
    obj.reshape!(*arg)
  end
  alias :shape= :reshape!
  def newdim!(*arg)
    @array = @array.newdim!(*arg)
    @mask = @mask.newdim!(*arg)
    self
  end
  alias :rewrank! :newdim!
  alias :rewrank= :newdim!
  def newdim(*arg)
    obj = self.dup
    obj.newdim!(*arg)
  end
  alias :rewrank :newdim


=begin
=== Type conversion
--- NArrayMiss#floor
     return (({NArrayMiss})) of integer whose elements processed (({floor})).
--- NArrayMiss#ceil
     return (({NArrayMiss})) of integer whose elements processed (({ceil})).
--- NArrayMiss#round
     return (({NArrayMiss})) of integer whose elements processed (({round})).
--- NArrayMiss#to_i
     return (({NArrayMiss})) of integer whose elements processed (({to_i})).
--- NArrayMiss#to_f
     return (({NArrayMiss})) of float whose elements processed (({to_f})).
--- NArrayMiss#to_type(typecode)
     return (({NArrayMiss})) of ((|typecode|)).
--- NArrayMiss#to_a
     convert (({NArrayMiss})) to (({Array})).
--- NArrayMiss#to_na!([missing_value])
     convert (({NArrayMiss})) to (({NArray})).
     if there is argument, set missing_value for invalid elements.
--- NArrayMiss#to_na([missing_value])
     convert (({NArrayMiss})) to (({NArray})).
     if there is argument, set missing_value for invalid elements.
--- NArrayMiss#to_s
     convert (({NArrayMiss})) to (({String})) as a binary data.
--- NArrayMiss#to_string
     create (({NArrayMiss})) of object whose elements are processed (({to_s}))
=end

  for operator in ["floor","ceil","round","to_i","to_f"]
    eval(<<-EOL,nil,__FILE__,__LINE__+1)
    def #{operator}
      NArrayMiss.to_nam_no_dup(@array.#{operator}, @mask.dup)
    end
    EOL
  end
  def to_type(typecode)
    NArrayMiss.to_nam_no_dup(@array.to_type(typecode), @mask.dup)
  end
  def to_a
    @array.to_a
  end
  def to_na!(*arg)
    if arg.length==0
      return @array
    elsif arg.length==1 then
      self.set_missing_value!(arg[0])
      return @array
    else
      raise(ArgumentError, "Usage: NArray#to_na([missing_value])")
    end
  end
  def to_na(*arg)
    return self.dup.to_na!(*arg)
  end
  def to_s
    @array.to_s
  end
  def to_string
    obj = NArrayMiss.obj(*@array.shape)
    obj.set_without_validation( @array.to_string )
    obh.set_mask(@mask)
    obj
  end


=begin
=== Iteration
--- NArrayMiss#each{|x| ...}
--- NArrayMiss#each_valid{|x| ...}
--- NArrayMiss#each_valid_with_index{|x,i| ...}
--- NArrayMiss#collect{|x| ...}
--- NArrayMiss#collect!{|x| ...}
=end

  def each
    for i in 0..self.total-1
      yield(@array[i])
    end
  end
  def each_valid
    for i in 0..self.total-1
      yield(@array[i]) if @mask[i]
    end
  end
  def each_valid_with_index
    for i in 0..self.total-1
      yield(@array[i],i) if @mask[i]
    end
  end
  def collect!
    for i in 0..self.total-1
      self[i] = yield(self[i])
    end
    self
  end
  def collect(&blk)
    self.dup.collect!(&blk)
  end


=begin
=== Boolean and mask related (only for byte, sint and int)
--- NArrayMiss#count_false
     return the number of elements whose value==0 and valid.
--- NArrayMiss#count_true
     return the number of elements whose value!=0 and valid.
--- NArrayMiss#mask(mask)
     return (({NArrayMiss#get_mask!&((|mask|))})).
--- NArrayMiss#all?
     return true if all the valid elements are not 0, else false.
--- NArrayMiss#any?
     return true if any valid element is not 0, else false.
--- NArrayMiss#none?
     return true if none of the valid elements is not 0, else false.
--- NArrayMiss#where
     return (({NArray})) of indices where valid elements are not 0.
--- NArrayMiss#where2
     return (({Array})) including two (({NArray}))s of indices,
     where valid elements are not 0, and 0, respectively.
=end

  def count_false
    if @array.typecode==BYTE then
      return @array.count_false-@mask.count_false
    else
      raise("cannot count_true NArrayMiss except BYTE type")
    end
  end
  def count_true
    if @array.typecode==BYTE then
      return (@array&@mask).count_true
    else
      raise("cannot count_true NArrayMiss except BYTE type")
    end
  end
  def mask(arg)
    obj = self.dup
    if arg.class==NArrayMiss then
      arg = arg.get_array!&arg.get_mask!
    end
    obj.set_mask(@mask&arg)
  end

  def all?
    @array[@mask].all?
  end
  def any?
    @array[@mask].any?
  end
  def none?
    @array[@mask].none?
  end

  def where
    (@array&@mask).where
  end
  def where2
    self.where-@mask.where
  end




=begin
=== Complex compound number (only for scomplex and complex)
--- NArrayMiss#real
--- NArrayMiss#imag
--- NArrayMiss#conj
--- NArrayMiss#angle
--- NArrayMiss#imag=(other)
--- NArrayMiss#im
=end

  def real
    NArrayMiss.to_nam_no_dup(@array.real,@mask)
  end
  def imag
    NArrayMiss.to_nam_no_dup(@array.imag,@mask)
  end
  def conj
    NArrayMiss.to_nam_no_dup(@array.conj,@mask)
  end
  def angle
    NArrayMiss.to_nam_no_dup(@array.angle,@mask)
  end
  def imag=(arg)
    @array.image=(arg)
    self
  end
  def im
    NArrayMiss.to_nam_no_dup(@array.im,@mask)
  end


=begin
=== Byte swap
--- NArrayMiss#swap_byte
     swap byte order.
--- NArrayMiss#hton
     convert to network byte order.
--- NArrayMiss#ntoh
     convert from network byte order.
--- NArrayMiss#htov
     convert to VAX byte order.
--- NArrayMiss#vtoh
     convert from VAX byte order.
=end

  def swap_byte
    obj = self.dup
    obj.set_without_validation(@array.swap_byte)
    obj
  end
  def hton
    NArray.to_nam(@array.hton,@mask.hton)
  end
  alias :ntoh :hton
  def htov
    NArray.to_nam(@array.htov,@mask.htov)
  end
  alias :vtoh :htov


=begin
=== Mask and missing value
--- NArrayMiss#set_valid(index)
     validate element at ((|index|)).
     ((|index|)) must be (({Integer, Range, Array, or ture})).
--- NArrayMiss#validation(index)
     alias to set_valid
--- NArrayMiss#set_invalid(index)
     invaliadate element at ((|index|)).
     ((|index|)) must be (({Integer, Range, Array, or ture})).
--- NArrayMiss#invalidation(index)
     alias to set_invalid
--- NArrayMiss#all_valid
     set all elements valid
--- NArrayMiss#all_invalid
     set all elements invalid
--- NArrayMiss#set_mask(mask)
     masking by ((|mask|))
--- NArrayMiss#set_missing_value(value)
     replace invalid elements by ((|value|)).
--- NArrayMiss#get_mask!
     return (({NArray})) of byte as mask.
--- NArrayMiss#get_mask
     return (({NArray})) of byte as mask.
--- NArrayMiss#get_array!
     return (({NArray})) as data.
--- NArrayMiss#get_array
     return (({NArray})) as data.
--- NArrayMiss#valid?
     return (({Array})) whose elements are true or false corresponding to valid or invalid of elements, respectively.
--- NArrayMiss#all_valid?
     return true if all elements are valid, else false.
--- NArrayMiss#none_valid?
     return true if all elements are invalid, else false.
--- NArrayMiss#all_invalid?
     alias to none_valid?
--- NArrayMiss#any_valid?
     return true if any elements are valid, else false.
     
--- NArrayMiss#count_valid
     return the number of valid elements.
--- NArrayMiss#count_invalid
     return the number of invalid elements.
=end

  def set_valid(*pos)
    @mask[*pos] = 1
  end
  alias validation set_valid
  def set_invalid(*pos)
    @mask[*pos] = 0
  end
  alias invalidation set_invalid
  def all_valid
    @mask[true]=1
    self
  end
  def all_invalid
    @mask[true]=0
    self
  end
  def set_mask(mask)
    if mask.class == Array then
      tmp = NArray.byte(*@mask.shape)
      tmp[true] = mask
      mask = tmp
    end
    if mask.class == NArrayMiss then
      mask = mask.to_na(0)
    end
    if mask.class == NArray then
      if mask.typecode != BYTE then
	raise("mask must be NArrayMiss.byte, NArray.byte or Array")
      end
      if @array.shape != mask.shape then
	raise("mask.shape must be same as array")
      end
      @mask = mask.dup
      return self
    else
      raise("mask must be NArray.byte or Array")
    end
  end

  def set_missing_value!(val)
    @array[@mask.not] = val
    self
  end
  def set_missing_value(val)
    obj = self.dup
    obj.set_missing_value!(val)
  end

  def get_mask!
    @mask
  end
  def get_mask
    @mask.dup
  end
  def get_array!
    @array
  end
  def get_array
    @array.dup
  end

  def valid?
    where = self.get_mask!.where2
    tf = Array.new(self.total)
    for i in where[0]
      tf[i] = true
    end
    for i in where[1]
      tf[i] = false
    end
    tf
  end
  def all_valid?
    @mask.all?
  end
  def none_valid?
    @mask.none?
  end
  alias :all_invalid? :none_valid?
  def any_valid?
    @mask.any?
  end

  def count_valid(*arg)
    if arg.length==0 then
      return @mask.count_true
    else
      return @mask.to_type(NArray::INT).sum(*arg)
    end    
  end
  def count_invalid(*arg)
    if arg.length==0 then
      return @mask.count_false
    else
      return NArray.int(*@mask.shape).fill(1).sum(*arg)-
	     @mask.to_type(NArray::INT).sum(*arg)
    end
  end


=begin
=== Others
--- NArrayMiss#integer?
     return true if (({NArrayMiss})) is byte, sint or int, else false.
--- NArrayMiss#complex?
     return true if (({NArrayMiss})) is scomplex or complex, else false.
--- NArrayMiss#dup
--- NArrayMiss#coerce(object)
--- NArrayMiss#inspect

go back to ((<Index>))
=end

  def integer?
    @array.integer?
  end
  def complex?
    @array.complex?
  end


  def dup
    NArrayMiss.to_nam(@array,@mask)
  end

  alias __clone__ clone
  def clone
    obj = __clone__
    obj.set_array(@array.clone)
    obj.set_mask(@mask.clone)
    return obj
  end

  def coerce(x)
    if Numeric===x then
      return [NArrayMiss.new(NArray[x].typecode,*self.shape).fill(x),self]
    elsif x.class==Array || x.class==NArray then
      return [NArrayMiss.to_nam(x), self]
    else
      raise("donnot know how to cange #{x.class} to NArrayMiss")
    end
  end


  def inspect
#    "array -> " + @array.inspect + "\nmask  -> " + @mask.inspect
    count_line = 0
    max_line = 10
    max_col = 80
    sep = ", "
    const = Hash.new
    NArray.constants.each{|c| const[NArray.const_get(c)] = c}
    str_ret = "NArrayMiss."+const[typecode].downcase+"("+shape.join(",")+"):"
    if rank == 0 then
      str_ret += " []"
      return str_ret
    else
       str_ret += "\n"
    end
    str = ""
    index = Array.new(rank,0)
    index[0] = true
    i = 1
    (rank-1).times{ str_ret += "[ " }
    while(true)
      i.times{ str_ret += "[ " }

      str = @array[*index].inspect
      ary = str[str.index("[")+1..str.index("]")-1].strip.split(/\s*,\s*/)
      miss = @mask[*index].where2[1]
      miss = miss[miss<ary.length].to_a
      if ary[-1]=="..." && miss[-1]==ary.length-1 then miss.pop end
      for j in miss
	ary[j] = "-"
      end
      while ( rank*4+ary.join(", ").length > max_col )
        ary.pop
        ary[-1] = "..."
      end
      str_ret += ary.join(", ")
      i = 1
      while (i<rank)
	if index[i]<shape[i]-1 then
	  str_ret += " ]"+sep+"\n"
	  count_line += 1
	  index[i] += 1
	  break
	else
	  str_ret += " ]"
	  index[i] = 0
	  i += 1
	end
      end

      if i>=rank then
	str_ret += " ]"
	return str_ret
      elsif count_line>=max_line then
	str_ret += " ..."
	return str_ret
      end

      (rank-i).times{ print("  ") }
    end
    return str_ret
  end


  def _dump(limit)
    Marshal::dump([@array._dump(nil),@mask._dump(nil)])
  end
  def self._load(o)
    ary, mask = Marshal::load(o)
    ary = NArray._load(ary)
    mask = NArray._load(mask)
    NArrayMiss.to_nam_no_dup(ary,mask)
  end




#  private
  private
  def pos(n,shape)
    rank = shape.length
    result = NArray.int(rank)
    m=n
    for i in 0..rank-2
      j = rank-1-i
      result[j] = m/shape[j-1]
      m = m%shape[j-1]
    end
    result[0] = m
    result
  end
  def routine1(arg,dummy)
    flag=true
    if Numeric===arg then
      term1 = @array
      term2 = arg
      mask = true
    elsif arg.class == Array then
      term1 = @array
      term1[@mask.not] = dummy
      term2 = NArray.to_na(arg)
      mask = NArray.byte(*term2.shape).fill(1)
    elsif arg.class == NArray then
      term1 = @array
      term1[@mask.not] = dummy
      term2 = arg
      mask = NArray.byte(*term2.shape).fill(1)
    elsif arg.class == NArrayMiss then
      mask = arg.get_mask
      term1 = @array
      term1[@mask.not] = dummy
      term2 = arg.to_na
      term2[arg.get_mask!.not] = dummy
    else
      term1,term2 = arg.coerce(self)
      flag=false
    end
    [term1,term2,mask,flag]
  end
  def routine2(result,mask,flag)
    if flag then
      obj = NArrayMiss.to_nam_no_dup(result)
      if mask==true then
	obj.set_mask(@mask)
      else
	mask = @mask+mask
	obj.set_mask(mask.eq(2))
      end
#      obj.set_without_validation(@mask.not,@array[@mask.not])
      if Numeric===obj && obj.get_mask![0]==1 then
	return obj.get_array[0]
      else
	return obj
      end
    else
      result
    end
  end
    
end



module NMMath


  func1 =  ["sqrt","exp","log","log10","log2",
            "sin","cos","tan","sinh","cosh","tanh",
            "asin","acos","atan","asinh","acosh","atanh"]
  func2 = ["atan2"]

  for operator in func1
    eval <<-EOL,nil,__FILE__,__LINE__+1
    def #{operator}(x)
      if Numeric===x || x.class==Array || x.class==NArray then
	NMath::#{operator}(x)
      elsif x.class == NArrayMiss then
	obj = NArrayMiss.new(x.typecode,*x.shape)
	mask = x.get_mask!
	obj[mask] = NMath::#{operator}(x.get_array![mask])
	obj[mask.not] = x[mask.not]
	obj.set_mask(mask)
        obj
      end
    end
    module_function :#{operator}
    EOL
  end

  for operator in func2
    eval <<-EOL,nil,__FILE__,__LINE__+1
    def #{operator}(x,y)
      if Numeric===x || x.class==Array || x.class==NArray then
	mask1 = nil
      elsif x.class == NArrayMiss then
	obj = NArrayMiss.new(x.typecode,*x.shape)
	mask1 = x.get_mask!
      end
      if Numeric===y || y.class==Array || y.class==NArray then
	mask2 = nil
      elsif y.class == NArrayMiss then
	obj = NArrayMiss.new(y.typecode,*y.shape)
	mask2 = y.get_mask!
      end
      if mask2.nil? then
	if mask1.nil? then
	  return NMath::#{operator}(x,y)
	else
	  obj[mask1] = NMath::#{operator}(x.get_array![mask1],y)
	  obj[mask1.not] = x[mask1.not]
	  obj.set_mask(mask1)
	  return obj
	end
      else
        if mask1.nil? then
          obj[mask2] = NMath::#{operator}(x,y.get_array![mask2])
          obj[mask2.not] = y[mask2.not]
          obj.set_mask(mask2)
          return obj
        else
          obj[mask1&mask2] = NMath::#{operator}(x.get_array![mask1],y.get_array![mask2])
          obj[(mask1&mask2).not] = y[(mask1&mask2).not]
          obj[(mask1&mask2).not] = x[(mask1&mask2).not]
          return obj
        end
      end
    end
    module_function :#{operator}
    EOL
  end

  for operator in func1+func2
    eval <<-EOL,nil,__FILE__,__LINE__+1
    def #{operator}(*x)
      x = [self]+x if NArrayMiss===self
      NMMath::#{operator}(*x)
    end
    EOL
  end
end

class NArrayMiss
  include NMMath
end