require "numo/narray"

p Numo::NArray::VERSION

z = Numo::DFloat.zeros(10)
p z

z = Numo::DFloat.zeros(10,10)
printf "%d bytes", z.byte_size

ri 'Numo::DFloat#+'

z = Numo::DFloat.zeros(10)
z[4] = 1
p z

z = Numo::DFloat[10..49]
p z

z = Numo::Int32.new(50).seq
z = z.reverse

z = Numo::Int32.new(3,3).seq
p z

nz = Numo::NArray[1,2,0,0,4,0].ne(0).where
p nz

z = Numo::DFloat.eye(3)
p z

z = Numo::DFloat.new(3,3,3).rand
p z

z = Numo::DFloat.new(10,10).rand
zmin, zmax = z.minmax
p zmin, zmax

z = Numo::DFloat.new(30).rand
m = z.mean
p m

z = Numo::DFloat.ones(10,10)
z[1..-2,1..-2] = 0
p z

0 * Float::NAN
Float::NAN == Float::NAN
Float::INFINITY > Float::NAN
Float::NAN - Float::NAN
0.3 == 3 * 0.1

z = Numo::Int32.zeros(5,5)
z.diagonal(-1)[] = Numo::Int32[1..4]
p z

# todo: rangewithstep
x = Numo::Int32.new(1,8).seq
y = Numo::Int32.new(8,1).seq
z = (x+y)%2
p z

# NArray allows unraveled index access
z = Numo::Int32.new(6,7,8).seq
p z[100]

z = Numo::DFloat.new(5,5).rand
zmin, zmax = z.minmax
z = (z - zmin)/(zmax - zmin)
p z

# todo: record
color = Numo::Struct.new do
  uint8 "r"
  uint8 "g"
  uint8 "b"
  uint8 "a"
end

x = Numo::DFloat.ones(5,3)
y = Numo::DFloat.ones(3,2)
z = x.dot y
p z

z = Numo::Int32.new(11).seq
z[(3 < z) & (z <= 8)] *= -1
p z

p [*0...5,-1].inject(:+)
p Numo::Int32[0...5].sum(-1)

z = Numo::Int32.new(5).seq
z**z
2 << z >> 2
z <- z
1i*z
z/1/1
z<z>z

p Numo::Int32[0] / Numo::Int32[0]
p Numo::DFloat[Float::NAN].cast_to(Numo::Int32).cast_to(Numo::DFloat)

z = Numo::DFloat.new(10).rand(-10,+10)
p (z + (0.5*z.sign)).trunc
# todo: copysign

# todo: intersect1d

a = Numo::DFloat.new(3).fill(1)
b = Numo::DFloat.new(3).fill(2)
p (a+b)*(-a/2)
(a+b.inplace)*(-a.inplace/2)
p b

z = Numo::DFloat.new(10).rand(10)

p z - z%1
p z.floor
p z.ceil - 1
p z.cast_to(Numo::Int32)
p z.trunc

z = Numo::DFloat.zeros(5,5)
z += Numo::Int32.new(5).seq
p z

z = Numo::DFloat.linspace(0,1,12)[1..-2]
p z

z = Numo::DFloat.new(10).rand
z = z.sort
p z

z = Numo::DFloat.zeros(10)
z.freeze
z[0] = 1

z = Numo::DFloat.new(10,2).rand
x,y = z[true,0], z[true,1]
r = Numo::NMath.sqrt(x**2+y**2)
t = Numo::NMath.atan2(y,x)
p r
p t

z = Numo::DFloat.new(10).rand
z[z.max_index] = 0
p z

z = Numo::Int32.new(100).seq
v = rand*100
index = (z-v).abs.min_index
p z[index]

require "stringio"
s = StringIO.new("1, 2, 3, 4, 5
                  6,  ,  , 7, 8
                   ,  , 9,10,11")
z = Numo::NArray[*s.readlines.map{|l| l.split(",").map{|x| x.strip.empty? ? Float::NAN : x.to_f}}]

z = Numo::Int32.new(3,3).seq
z.each_with_index{|x,*i| p [i,x]}

x = Numo::DFloat.linspace(-1,1,10)
y = Numo::DFloat.linspace(-1,1,10).expand_dims(1)
d = Numo::NMath.sqrt(x*x+y*y)
sigma, mu = 1.0, 0.0
g = Numo::NMath.exp(-( (d-mu)**2 / ( 2.0 * sigma**2 ) ) )
p g

x = Numo::DFloat.new(5, 10).rand
y = x - x.mean(1).expand_dims(1)

z = Numo::Int32.new(3,3).rand(10)
p z
p z[z[true,1].sort_index,true]

z = Numo::Int32.new(3,10).rand(3)
(~z.ne(0).any?(0)).any?

z = Numo::DFloat.new(10).rand
x = 0.5
m = z[(z - x).abs.min_index]
p m

# todo: unique
w,h = 16,16
i = Numo::UInt32.new(h,w,3).rand(2)
f = i[false,0]*256*256 + i[false,1]*256 +i[false,2]
p f.flatten.sort.to_a.uniq

a = Numo::Int32.new(3,4,3,4).rand(10)
sum = a.sum(-2,-1)
p sum

a = Numo::DFloat.new(3,3).seq
b = Numo::DFloat.new(3,3).seq
p a.mulsum(b.transpose,1)
# speed?

z = Numo::NArray[1,2,3,4,5]
nz = 3
z0 = Numo::Int32.zeros(z.size + (z.size-1)*(nz))
# todo: rangewithstep
# z0[(0..-1).step(nz+1)] = z
# p z0

a = Numo::Int32.ones(5,5,3)
b = Numo::Int32.new(5,5).fill(2)
p a * b[:*,:*,:-]

a = Numo::Int32.new(5,5).seq
a[[0,1],true] = a[[1,0],true].copy
p a
# todo: identity check between read/write array

def moving_average(a, n=3)
  ret = a.cumsum
  ret[n..-1] = ret[n..-1] - ret[0..-n-1]
  ret[n-1..-1] / n
end
z = Numo::DFloat.new(20).seq
p moving_average(z, 3)

# todo: logical_not
z = Numo::Int32.new(100).rand(2)
p z
z.inplace ^ 1
p z

z = Numo::DFloat.new(100).rand(-1,1)
p z
-z.inplace
p z

def distance(p0, p1, p)
  t = p1 - p0
  l = (t**2).sum(1)
  u = -((p0[true,0]-p[false,0])*t[true,0] + (p0[true,1]-p[false,1])*t[true,1]) / l
  u = u.reshape(u.size,1)
  d = p0 + u*t - p
  return Numo::NMath.sqrt((d**2).sum(1))
end

p0 = Numo::DFloat.new(10,2).rand(-10,10)
p1 = Numo::DFloat.new(10,2).rand(-10,10)
p  = Numo::DFloat.new( 1,2).rand(-10,10)
p distance(p0, p1, p)

p0 = Numo::DFloat.new(10,2).rand(-10,10)
p1 = Numo::DFloat.new(10,2).rand(-10,10)
p  = Numo::DFloat.new(10,2).rand(-10,10)
p a = p.shape[0].times.map{|i| distance(p0, p1, p.slice(i,true))}
# todo: concat narray

module Symetric
  def []=(i,j,value)
    super(i,j,value)
    super(j,i,value) if i != j
  end
end

def symetric(z)
  y = z + z.transpose
  y.diagonal.store(z.diagonal)
  y.extend(Symetric)
end

s = symetric(Numo::Int32.new(5,5).rand(10))
s[2,3] = 42
p s

p, n = 10, 20
m = Numo::DFloat.ones(p,n,n)
v = Numo::DFloat.ones(p,n,1)
s = m.transpose(0,2,1).mulsum(v,0,1)
p s
# todo: tensordot?

n, k = 16, 4
z = Numo::DFloat.ones(n,n)
s = z.reshape(n/k,k,n/k,k).sum(1,3)
# todo: reduceat?

def iterate(z)
  # Count neighbours
  n = z[0..-3,0..-3] + z[0..-3,1..-2] + z[0..-3,2..-1] +
      z[1..-2,0..-3]                  + z[1..-2,2..-1] +
      z[2..-1,0..-3] + z[2..-1,1..-2] + z[2..-1,2..-1]

  # Apply rules
  birth = n.eq(3) & z[1..-2,1..-2].eq(0)
  survive = (n.eq(2) | n.eq(3)) & z[1..-2,1..-2].eq(1)
  z[] = 0
  #z[1..-2,1..-2][birth | survive] = 1
  y = z[0..-3,0..-3].copy
  y[birth | survive] = 1
  z[1..-2,1..-2] = y
end

z = Numo::Int32.new(50,50).rand(2)
100.times{ iterate(z) }
p z

z = Numo::DFloat.new(10000).rand
n = 5
p z[z.sort_index[-n..-1]]
# todo: shuffle, argpartition

def cartesian(*arrays)
  arrays = arrays.map{|a| Numo::Int32.cast(a)}
  shape = arrays.map{|x| x.size}
  asz = arrays.size

  ix = Numo::Int32.zeros(*shape, asz)
  arrays.each_with_index do |arr,n|
    s = [1]*asz
    s[n] = arr.size
    ix[false,n] = arr.reshape(*s)
  end
  return ix.reshape(ix.size/asz,asz)
end

p cartesian([1, 2, 3], [4, 5], [6, 7])

x = Numo::DFloat.new(5e7).rand
x**3 # probably fast

a = Numo::Int32.new(8,3).rand(5)
b = Numo::Int32.new(2,2).rand(5)
c = a[false,:new,:new].eq b
rows = (c.count_true(1,2,3) >= b.shape[1]).where
p rows

z = Numo::Int32.new(10,3).rand(5)
e = (z[true,1..-1].eq z[true,0..-2]).all?(1)
u = z[(~e).where,true]
p z
p u

# no method: einsum
a = Numo::DFloat.new(10).rand(0,1)
b = Numo::DFloat.new(10).rand(0,1)

a.sum           # np.sum(A)
a*b             # A * B
a.mulsum(b)     # np.inner(A, B)
a[false,:new]*b # np.outer(A, B)

x = Numo::DFloat.new(100).rand
n = 1000 # number of bootstrap samples
idx = Numo::Int32.new(n, x.size).rand(x.size)
means = x[idx].mean(1)
confint = means[means.sort_index[means.size/100.0*Numo::DFloat[2.5, 97.5]]]
p confint
# todo: percentile, rand_norm