bug: fixed transpose function and fixed reciprocal-cell calculation

Now the transpose function works for 2D and 3D arrays.

The reciprocal cell calculation was wrong for skewed lattices.

Signed-off-by: Nick Papior <[email protected]>

flos/num/array.lua

@@ -62,7 +62,7 @@ end
 
 
 function Array:initialize(...)
- 
+
  local sh = nil
 
  -- One may also initialize an Array by passing a
@@ -75,7 +75,7 @@ function Array:initialize(...)
  if sh == nil then
  sh = shape.Shape(...)
  end
- 
+
  -- Create the shape container
  rawset(self, "shape", sh)
 
@@ -200,7 +200,7 @@ function Array.range(i1, i2, step)
  elseif is < 0 and i1 < i2 then
  error("flos.Array range with negative step-length and i1 < i2 is not allowed.")
  end
- 
+
  local new = Array.empty( 1 )
  local j = 0
  for i = i1, i2, is do
@@ -266,7 +266,7 @@ function Array:set_linear(i, v)
 
  -- If we are at the last dimension, return immediately.
  if #self.shape == 1 then
- 
+
  self[i] = v
 
  return
@@ -278,7 +278,7 @@ function Array:set_linear(i, v)
  local j = m.tointeger( m.ceil(i / n_dim) )
  -- Transform i into the linear index in the underlying array
  self[j]:set_linear( m.tointeger(i - (j-1) * n_dim), v)
- 
+
 end
 
 
@@ -323,7 +323,7 @@ function Array:reshape(...)
  if #arg == 0 then
  arg[1] = 0
  end
- 
+
  -- In case a shape is passed
  local sh
  if shape.isShape(arg[1]) then
@@ -334,7 +334,7 @@ function Array:reshape(...)
  if sh == nil then
  error("flos.Array cannot align shapes, incompatible dimensions")
  end
- 
+
  -- Create the new array
  local new = Array( sh )
 
@@ -522,9 +522,9 @@ function Array:norm(axis)
  norm = norm + self[i]:norm(0) ^ 2
  end
  norm = m.sqrt(norm)
- 
+
  elseif #self.shape == 1 then
- 
+
  norm = 0.
  for i = 1, #self do
  norm = norm + self[i] * self[i]
@@ -539,9 +539,9 @@ function Array:norm(axis)
  for i = 1, #norm do
  norm[i] = self[i]:norm()
  end
- 
+
  end
- 
+
  return norm
 end
 
@@ -560,14 +560,14 @@ function Array:scalar_project(P, axis)
  local ax = ax_(axis)
 
  if ax == 0 then
- 
+
  -- Calculate norm of the projection vector
  return self:flatten():dot( P:flatten() ) / P:norm(0)
- 
+
  else
- 
+
  error("flos.Array could not project on anything but flattened array")
- 
+
  end
 
 end
@@ -587,15 +587,15 @@ function Array:project(P, axis)
  local ax = ax_(axis)
 
  if ax == 0 then
- 
+
  -- Calculate norm of the projection vector
  local dnorm2 = P:norm(0) ^ 2
  return self:flatten():dot( P:flatten() ) / dnorm2 * P
- 
+
  else
- 
+
  error("flos.Array could not project on anything but flattened array")
- 
+
  end
 
 end
@@ -631,7 +631,7 @@ function Array:min(axis)
 
  error("NotimplementedYet")
  end
- 
+
  return min
 end
 
@@ -666,7 +666,7 @@ function Array:max(axis)
 
  error("NotimplementedYet")
  end
- 
+
  return max
 end
 
@@ -679,7 +679,7 @@ function Array:sum(axis)
 
  local sum
  if ax == 0 then
- 
+
  -- Special case for the 1D case
  if #self.shape == 1 then
  sum = self[1]
@@ -753,7 +753,7 @@ function Array.cross(lhs, rhs)
  local cross = Array( sh )
 
  if #cross.shape == 1 then
- 
+
  cross[1] = lhs[2] * rhs[3] - lhs[3] * rhs[2]
  cross[2] = lhs[3] * rhs[1] - lhs[1] * rhs[3]
  cross[3] = lhs[1] * rhs[2] - lhs[2] * rhs[1]
@@ -814,7 +814,7 @@ function Array.dot(lhs, rhs)
  if lhs.shape ~= rhs.shape then
  error("flos.Array dot dimensions for 1D dot product are not the same")
  end
- 
+
  -- This is a element wise product and sum
  dot = 0.
  for i = 1, #lhs do
@@ -872,11 +872,11 @@ function Array.dot(lhs, rhs)
  else
 
  error("flos.Array dot for arrays with anything but 1 or 2 dimensions is not implemented yet")
- 
+
  end
 
  return dot
- 
+
 end
 
 
@@ -885,24 +885,43 @@ end
 function Array:transpose()
 
  -- Check dimensions, we cannot transpose a 1D array
- if #self.shape == 1 then
- error("flos.Array cannot transpose a vector, reshape, then transpose")
- end
+ local nd = #self.shape
+ local new = nil
 
- -- First create the reversed shape
- local sh = self.shape:reverse()
+ if nd == 1 then
+  return self:copy()
 
- -- Create return array
- local new = Array( sh )
+ elseif nd == 2 then
+
+ -- Create return array
+ new = Array( self.shape:reverse() )
+
+ for i = 1 , self.shape[1] do
+ for j = 1 , self.shape[2] do
+ new[j][i] = self[i][j]
+ end
+ end
+
+ elseif nd == 3 then
+
+ -- Create return array
+ new = Array( self.shape:reverse() )
+
+ for i = 1 , self.shape[1] do
+ for j = 1 , self.shape[2] do
+ for k = 1 , self.shape[3] do
+ new[k][j][i] = self[i][j][k]
+ end
+ end
+ end
+
+ else
+
+ error("flos.Array transpose only works up to 3D arrays")
 
- -- Perform transpose
- local size = self:size()
- for i = 1, size do
- new:set_linear(size-i+1, self:get_linear(i))
  end
 
  return new
-
 end
 
 
@@ -925,14 +944,14 @@ function Array.__add(lhs, rhs)
  if sh == nil then
  error("flos.Array + requires the same shape for two different Arrays")
  end
- 
+
  -- Create the return array
  ret = Array( lhs.shape )
  -- Element-wise additions
  for i = 1, #lhs do
  ret[i] = lhs[i] + rhs[i]
  end
- 
+
  elseif isArray(lhs) then
 
  ret = Array( lhs.shape )
@@ -943,7 +962,7 @@ function Array.__add(lhs, rhs)
  end
 
  elseif isArray(rhs) then
- 
+
  -- Add scalar to all RHS
  ret = Array( rhs.shape )
 
@@ -1135,7 +1154,7 @@ function Array.__pow(lhs, rhs)
  if sh == nil then
  error("flos.Array ^ requires the same shape for two different Arrays")
  end
- 
+
  ret = Array( lhs.shape )
  for i = 1, #lhs do
  ret[i] = lhs[i] ^ rhs[i]
@@ -1145,7 +1164,7 @@ function Array.__pow(lhs, rhs)
  if type(rhs) == "string" and rhs == "T" then
  return lhs:transpose()
  end
- 
+
  ret = Array( lhs.shape )
  for i = 1, #lhs do
  ret[i] = lhs[i] ^ rhs
@@ -1155,7 +1174,7 @@ function Array.__pow(lhs, rhs)
  if type(lhs) == "string" and lhs == "T" then
  return rhs:transpose()
  end
- 
+
  ret = Array( rhs.shape )
  for i = 1, #rhs do
  ret[i] = lhs ^ rhs[i]

flos/optima/lattice.lua

@@ -67,11 +67,7 @@ function Lattice:update_reciprocal()
  self.rcell[3][3] = c[1][1]*c[2][2] - c[1][2]*c[2][1]
 
  for i = 1, 3 do
- local n = c[i][1]*self.rcell[i][1] +
- c[i][2]*self.rcell[i][2] +
- c[i][3]*self.rcell[i][3]
-
- self.rcell[i] = self.rcell[i] / n
+ self.rcell[i] = self.rcell[i] / c[i]:dot(self.rcell[i])
  end
 
 end