From f2ca3c596964fbea42a2a6d9b681b81de3b0c08b Mon Sep 17 00:00:00 2001
From: James Osborn <osborn@alcf.anl.gov>
Date: Wed, 13 Nov 2024 23:28:56 -0600
Subject: [PATCH] initial implementation for Color

---
 .github/workflows/test.yml  |   3 +-
 src/io/qioInternal.nim      |  32 +++-
 src/maths/matrixConcept.nim |  25 ++-
 src/physics/color.nim       | 292 +-----------------------------
 src/physics/colorOld.nim    | 289 ++++++++++++++++++++++++++++++
 src/physics/colorTensor.nim |  57 ++++++
 src/physics/tensorwrap.nim  | 347 ++++++++++++++++++++++++++++++++++++
 7 files changed, 746 insertions(+), 299 deletions(-)
 create mode 100644 src/physics/colorOld.nim
 create mode 100644 src/physics/colorTensor.nim
 create mode 100644 src/physics/tensorwrap.nim

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 8cd40fd..4afcadc 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -5,13 +5,14 @@ on:
     branches:
       - 'devel'
       - 'master'
+      - 'tensorwrap'
 
 jobs:
   build:
     strategy:
       matrix:
         mpi-impl: [openmpi, mpich]
-        nim-branch: [version-1-6, version-2-0, devel]
+        nim-branch: [version-2-0, version-2-2, devel]
         fuel-compat: [0, 1]
       fail-fast: false
     name: nim-${{ matrix.nim-branch }}-${{ matrix.mpi-impl }}-FUELCompat:${{ matrix.fuel-compat }}
diff --git a/src/io/qioInternal.nim b/src/io/qioInternal.nim
index 5a5a1b8..139c34c 100644
--- a/src/io/qioInternal.nim
+++ b/src/io/qioInternal.nim
@@ -44,10 +44,30 @@ template IOname*[T](x:typedesc[T]):string =
 #  "QDP_F" & $N & "_ColorMatrix"
 #template IOname*[N:static int](x:typedesc[Color[MatrixArray[N,N,DComplex]]]):string =
 #  "QDP_D" & $N & "_ColorMatrix"
-template IOname*[N:static int,T](x:typedesc[Color[MatrixArray[N,N,T]]]):string =
-  when T is SComplex:
-    "QDP_F" & $N & "_ColorMatrix"
-  elif T is DComplex:
-    "QDP_D" & $N & "_ColorMatrix"
+# template IOname*[N:static int,T](x:typedesc[Color[MatrixArray[N,N,T]]]):string =
+#   static: echo $x.type
+#   when T is SComplex:
+#     "QDP_F" & $N & "_ColorMatrix"
+#   elif T is DComplex:
+#     "QDP_D" & $N & "_ColorMatrix"
+#   else:
+#     IOnameDefault type T
+# #template IOname*(x:typedesc[Color[AsMatrix]]):string =
+#  static: echo $x.type
+#  when x.index(int,int) is SComplex:
+#    "QDP_F" & $x.Nc & "_ColorMatrix"
+#  elif x.index(int,int) is DComplex:
+#    "QDP_D" & $x.Nc & "_ColorMatrix"
+#  else:
+#    IOnameDefault x
+template IOname*[T](x:typedesc[Color[T]]):string =
+  mixin Nc
+  when T is AsMatrix:
+    when x.index(int,int) is SComplex:
+      "QDP_F" & $x.getNc & "_ColorMatrix"
+    elif x.index(int,int) is DComplex:
+      "QDP_D" & $x.getNc & "_ColorMatrix"
+    else:
+      IOnameDefault x
   else:
-    IOnameDefault type T
+    IOnameDefault x
diff --git a/src/maths/matrixConcept.nim b/src/maths/matrixConcept.nim
index 21355db..b5be25c 100644
--- a/src/maths/matrixConcept.nim
+++ b/src/maths/matrixConcept.nim
@@ -41,6 +41,9 @@ template createAsType2(t,c:untyped) =
       c(index(X.type[], type T))
     else:
       index(X[], T)
+  template index*[X:t](x: typedesc[X]; i,j: typedesc[int]): typedesc =
+    mixin index
+    index(X[], int, int)
   template `[]`*(x: t; i: Scalar): untyped = c(x[][i])
   template `[]`*(x:t; i,j:SomeInteger):untyped =
     #echoType: x
@@ -65,10 +68,12 @@ template createAsType2(t,c:untyped) =
       x[][i] = y
   template `[]=`*(x: t; i,j: SomeInteger; y: typed) =
     x[][i,j] = y
-  template len*(x:t):untyped = getConst(x[].len)
-  template nrows*(x:t):untyped = getConst(x[].nrows)
-  template ncols*(x:t):untyped = getConst(x[].ncols)
-  template len*(x:typedesc[t]):auto = getConst(x[].len)
+  template len*(x:t):int = getConst(x[].len)
+  template nrows*(x:t):int = getConst(x[].nrows)
+  template ncols*(x:t):int = getConst(x[].ncols)
+  template len*(x:typedesc[t]):int = getConst(x[].len)
+  template nrows*(x:typedesc[t]):int = getConst(x[].nrows)
+  template ncols*(x:typedesc[t]):int = getConst(x[].ncols)
   #template mvLevel*(x:t):untyped =
   #  mixin mvLevel
   #  mvLevel(x[])
@@ -99,6 +104,9 @@ template has*[T:AsMatrix](x: typedesc[T], y: typedesc): bool =
   when y is AsMatrix: true
   else: has(T.type[], y)
 
+template `*`*(x: typedesc[AsMatrix], y: typedesc[AsVector]): typedesc =
+  asVector(eval(x[]) * eval(y[]))
+
 #declareScalar(AsScalar)
 #declareScalar(AsVarScalar)
 #declareVector(AsVector)
@@ -226,7 +234,12 @@ template index*[I,J:static[int],T,K](x: typedesc[MatrixArrayObj[I,J,T]];
   else:
     false # error
 
+template index*[I,J:static[int],T](x: typedesc[MatrixArrayObj[I,J,T]];
+                                   k,l: typedesc[int]): typedesc =
+  T
+
 template `len`*(x:MatrixArrayObj):untyped = x.I
+template nrows*(x:typedesc[MatrixArrayObj]):int = x.I
 template nrows*(x:MatrixArrayObj):untyped = x.I
 template ncols*(x:MatrixArrayObj):untyped = x.J
 template `[]`*(x:MatrixArrayObj):untyped = x.mat
@@ -408,6 +421,10 @@ proc setColumn*(r:var AsMatrix; x:AsVector; i:int) {.inline.} =
   for j in 0..<nr:
     assign(r[j,i], x[j])
 
+template `*`*[R,C:static[int],T1,T2](x: typedesc[MatrixArrayObj[R,C,T1]],
+                                     y: typedesc[VectorArrayObj[C,T2]]): typedesc =
+  VectorArray[R,eval(T1)*eval(T2)]
+
 import matrixOps
 export matrixOps
 
diff --git a/src/physics/color.nim b/src/physics/color.nim
index 91dc40d..6170453 100644
--- a/src/physics/color.nim
+++ b/src/physics/color.nim
@@ -1,289 +1,5 @@
-import base/basicOps
-import base/wrapperTypes
-export wrapperTypes
-import maths/types
-import maths
-import simd/simdWrap
+#import colorOld
+#export colorOld
 
-makeWrapperType(Color):
-  ## wrapper type for colored objects
-
-type
-  Color2*[T] = Color[T]
-  Color3*[T] = Color[T]
-  Color4*[T] = Color[T]
-
-template asVarWrapper*(x: Color, y: typed): untyped =
-  #static: echo "asVarWrapper Color"
-  #var cy = asColor(y)
-  #cy
-  asVar(asColor(y))
-
-template index*[T,I](x: typedesc[Color[T]], i: typedesc[I]): typedesc =
-  when I is Color:
-    index(T.type, I.type[])
-  elif I.isWrapper:
-    Color[index(T.type, I.type)]
-  else:
-    index(T.type, I.type)
-
-template `[]`*[T](x: Color, i: T): untyped =
-  when T is Color:
-    x[][i[]]
-  elif T.isWrapper:
-    #indexed(x, i)
-    var tColorBracket = asColor(x[][i])
-    tColorBracket
-  else:
-    x[][i]
-template `[]`*(x: Color, i,j: SomeInteger): auto = x[][i,j]
-
-#[
-template `[]=`*[T](x: Color, i:T, y: typed): untyped =
-  when T.isWrapper and T isnot Color:
-    x[][asScalar(i)] = y
-  else:
-    x[][i] = y
-]#
-template `[]=`*[T](x: Color, i: T; y: auto) =
-  when T is Color2:
-    x[][i[]] = y
-  elif T.isWrapper:
-    #indexed(x, i)
-    var tColorBracket = asColor(x[][i])
-    tColorBracket := y
-  else:
-    x[][i] = y
-template `[]=`*(x: Color, i,j: SomeInteger, y: auto) =
-  x[][i,j] = y
-
-# forward from value to value
-template forwardVV(f: untyped) {.dirty.} =
-  template f*(x: Color): auto =
-    mixin f
-    f(x[])
-# forward from value to type
-#template forwardVT(f: untyped) {.dirty.} =
-#  template f*[T](x: Color[T]): untyped =
-#    mixin f
-#    f(type T)
-# forward from type to type
-template forwardTT(f: untyped) {.dirty.} =
-  template f*[T](x: typedesc[Color[T]]): auto =
-    mixin f
-    f(type T)
-# forward from type to type and wrap
-template forwardTTW(f: untyped) {.dirty.} =
-  template f*[T](x: typedesc[Color[T]]): auto =
-    mixin f
-    Color[f(type T)]
-
-forwardVV(len)
-forwardVV(nrows)
-forwardVV(ncols)
-forwardVV(nVectors)
-forwardVV(simdType)
-#forwardVV(simdLength)
-forwardVV(getNs)
-forwardVV(numNumbers)
-
-forwardTT(len)
-forwardTT(nrows)
-forwardTT(ncols)
-forwardTT(nVectors)
-forwardTT(simdType)
-forwardTT(simdLength)
-forwardTT(getNs)
-forwardTT(numberType)
-
-forwardTTW(toSingle)
-forwardTTW(toDouble)
-
-#template eval*[T](x: typedesc[Color[T]]): typedesc = asColor(eval(type T))
-template eval*[T:Color](x: typedesc[T]): typedesc = asColor(eval((type T)[]))
-
-template has*[T:Color](x: typedesc[T], y: typedesc): bool =
-  mixin has
-  when y is Color: true
-  else: has(T.type[], y)
-
-template row*(x: Color, i: untyped): untyped =
-  mixin row
-  asColor(row(x[],i))
-template setRow*(r: Color; x: Color2; i: untyped): untyped =
-  setRow(r[], x[], i)
-
-template getNc*[T](x: Color[T]): untyped =
-  when T is Mat1:
-    x[].nrows
-  elif T is Vec1:
-    x[].len
-  else:
-    static:
-      echo "error: unknown Nc"
-      echo x.repr
-      echo type(x).name
-      qexExit 1
-    0
-
-template binDDRet(fn,wr,T1,T2) =
-  template fn*(x: T1, y: T2): untyped =
-    wr(fn(x[], y[]))
-  #template fn*(x: T1, y: T2): auto =
-  #  var tmp {.noInit.}: wr(evalType(fn(x[],y[]))
-  #  wr(fn(x[], y[]))
-
-#binDDRet(`+`, asColor, Color, Color2)
-binDDRet(`-`, asColor, Color, Color2)
-#binDDRet(`*`, asColor, Color, Color2)
-binDDRet(`/`, asColor, Color, Color2)
-
-setBinop(`+`, add, Color, Color2, asColor(evalType(x[]+y[])))
-setBinop(`*`, mul, Color, Color2, asColor(evalType(x[]*y[])))
-
-template load1*(x: Color): untyped = asColor(load1(x[]))
-template `-`*(x: Color): untyped = asColor(-(x[]))
-template assign*(r: var Color, x: SomeNumber) =
-  assign(r[], x)
-template assign*(r: var Color, x: AsComplex) =
-  assign(r[], x)
-template assign*(r: var Color, x: Color2) =
-  assign(r[], x[])
-template `:=`*(r: var Color, x: SomeNumber) =
-  r[] := x
-template `:=`*(r: var Color, x: AsComplex) =
-  r[] := x
-template `:=`*(r: var Color, x: Color2) =
-  r[] := x[]
-template `+=`*(r: var Color, x: Color2) =
-  r[] += x[]
-template `-=`*(r: var Color, x: Color2) =
-  r[] -= x[]
-template `*=`*(r: var Color, x: SomeNumber) =
-  r[] *= x
-template iadd*(r: var Color, x: SomeNumber) =
-  iadd(r[], x)
-template iadd*(r: var Color, x: AsComplex) =
-  iadd(r[], x)
-template iadd*(r: var Color, x: Color2) =
-  iadd(r[], x[])
-template isub*(r: var Color, x: SomeNumber) =
-  isub(r[], x)
-template isub*(r: var Color, x: Color2) =
-  isub(r[], x[])
-template imul*(r: var Color, x: SomeNumber) =
-  imul(r[], x)
-template imadd*(r: var Color, x: Color2, y: Color3) =
-  imadd(r[], x[], y[])
-template imadd*(r: var Color, x: AsComplex, y: Color3) =
-  imadd(r[], x, y[])
-template imsub*(r: var Color, x: Color2, y: Color3) =
-  imsub(r[], x[], y[])
-template peqOuter*(r: var Color, x: Color2, y: Color3) =
-  peqOuter(r[], x[], y[])
-template meqOuter*(r: var Color, x: Color2, y: Color3) =
-  meqOuter(r[], x[], y[])
-template `+`*(r: Color, x: SomeNumber): untyped =
-  asColor(r[] + x)
-template `+`*(r: SomeNumber, x: Color): untyped =
-  asColor(r + x[])
-template `-`*(r: Color, x: SomeNumber): untyped =
-  asColor(r[] - x)
-template `+`*(r: Color, x: AsComplex): untyped =
-  asColor(r[] + x)
-template `-`*(r: Color, x: AsComplex): untyped =
-  asColor(r[] - x)
-template add*(r: var Color, x: Color2, y: Color3) =
-  add(r[], x[], y[])
-template sub*(r: var Color, x: Color2, y: Color3) =
-  sub(r[], x[], y[])
-template `*`*(x: Color, y: SomeNumber): untyped =
-  asColor(x[] * y)
-#template `*`*(x: SomeNumber, y: Color): auto =
-  #asColor(x * y[])
-template `*`*[X:SomeNumber,Y:Color](x: X, y: Y): auto =
-  #static: echo "SomeNumber * Color"
-  #var tmp {.noInit.}: asColor(type(X)*type(Y)[])
-  var tmp {.noInit.}: asColor(evalType(x*y[]))
-  #static: echo $type(tmp)
-  mul(tmp[], x, y[])
-  tmp
-template `*`*(x: Simd, y: Color2): untyped =
-  asColor(x * y[])
-#template `*`*(x: AsReal, y: Color2): untyped =
-#  asColor(x * y[])
-template `*`*[X:AsReal,Y:Color](x: X, y: Y): auto =
-  #var tmp {.noInit.}: asColor(type(X)*type(Y)[])
-  var tmp {.noInit.}: asColor(evalType(x*y[]))
-  mul(tmp[], x, y[])
-  tmp
-template `*`*(x: AsImag, y: Color2): untyped =
-  asColor(x * y[])
-template `*`*(x: AsComplex, y: Color2): untyped =
-  asColor(x * y[])
-template `*`*(x: Color, y: AsComplex): untyped =
-  asColor(x[] * y)
-#template mul*(x: SomeNumber, y: Color2): untyped =
-#  asColor(`*`(x, y[]))
-template mul*[X:SomeNumber,Y:Color](x: X, y: Y): auto =
-  #var tmp {.noInit.}: asColor(type(X)*type(Y)[])
-  var tmp {.noInit.}: asColor(evalType(x*y[]))
-  mul(tmp[], x, y[])
-  tmp
-template mul*(x: Color, y: Color2): untyped =
-  asColor(`*`(x[], y[]))
-template mul*(r: var Color, x: Color2, y: Color3) =
-  mul(r[], x[], y[])
-template mul*(r: var Color, x: SomeNumber, y: Color3) =
-  mul(r[], x, y[])
-template mul*(r: var Color, x: Color2, y: SomeNumber) =
-  mul(r[], x[], y)
-template mul*(r: var Color, x: AsComplex, y: Color3) =
-  mul(r[], x, y[])
-template mul*(x: AsComplex, y: Color2): untyped =
-  asColor(mul(x, y[]))
-template random*(x: var Color) =
-  gaussian(x[], r)
-template gaussian*(x: Color, r: untyped) =
-  gaussian(x[], r)
-template uniform*(x: var Color, r: var untyped) =
-  uniform(x[], r)
-template z4*(x: var Color, r: var untyped) =
-  z4(x[], r)
-template z2*(x: var Color, r: var untyped) =
-  z2(x[], r)
-template u1*(x: var Color, r: var untyped) =
-  u1(x[], r)
-template projectU*(r: var Color) =
-  projectU(r[])
-template projectU*(r: var Color, x: Color2) =
-  projectU(r[], x[])
-template projectUderiv*(r: var Color, u: Color2, x: Color3, chain: Color4) =
-  projectUderiv(r[], u[], x[], chain[])
-template projectUderiv*(r: var Color, x: Color3, chain: Color4) =
-  projectUderiv(r[], x[], chain[])
-template projectSU*(r: var Color) =
-  projectSU(r[])
-template projectSU*(r: var Color, x: Color2) =
-  projectSU(r[], x[])
-template projectTAH*(r: var Color) =
-  projectTAH(r[])
-template projectTAH*(r: var Color, x: Color2) =
-  projectTAH(r[], x[])
-template checkU*(x: Color):untyped = checkU(x[])
-template checkSU*(x: Color):untyped = checkSU(x[])
-template norm2*(x: Color): untyped = norm2(x[])
-template norm2*(r: var auto, x: Color): untyped = norm2(r, x[])
-template inorm2*(r: var auto, x: Color2) = inorm2(r, x[])
-template dot*(x: Color, y: Color2): untyped =
-  dot(x[], y[])
-template idot*(r: var auto, x: Color2, y: Color3) = idot(r, x[], y[])
-template redot*(x: Color, y: Color2): untyped =
-  redot(x[], y[])
-template trace*(x: Color): untyped = trace(x[])
-template simdSum*(x: Color): untyped = asColor(simdSum(x[]))
-template re*(x: Color): untyped = asColor(re(x[]))
-template im*(x: Color): untyped = asColor(im(x[]))
-template exp*(x: Color): untyped = asColor(exp(x[]))
-template expDeriv*(x: Color, c: Color2): untyped = asColor(expDeriv(x[], c[]))
-template ln*(x: Color): untyped = asColor(ln(x[]))
+import colorTensor
+export colorTensor
diff --git a/src/physics/colorOld.nim b/src/physics/colorOld.nim
new file mode 100644
index 0000000..91dc40d
--- /dev/null
+++ b/src/physics/colorOld.nim
@@ -0,0 +1,289 @@
+import base/basicOps
+import base/wrapperTypes
+export wrapperTypes
+import maths/types
+import maths
+import simd/simdWrap
+
+makeWrapperType(Color):
+  ## wrapper type for colored objects
+
+type
+  Color2*[T] = Color[T]
+  Color3*[T] = Color[T]
+  Color4*[T] = Color[T]
+
+template asVarWrapper*(x: Color, y: typed): untyped =
+  #static: echo "asVarWrapper Color"
+  #var cy = asColor(y)
+  #cy
+  asVar(asColor(y))
+
+template index*[T,I](x: typedesc[Color[T]], i: typedesc[I]): typedesc =
+  when I is Color:
+    index(T.type, I.type[])
+  elif I.isWrapper:
+    Color[index(T.type, I.type)]
+  else:
+    index(T.type, I.type)
+
+template `[]`*[T](x: Color, i: T): untyped =
+  when T is Color:
+    x[][i[]]
+  elif T.isWrapper:
+    #indexed(x, i)
+    var tColorBracket = asColor(x[][i])
+    tColorBracket
+  else:
+    x[][i]
+template `[]`*(x: Color, i,j: SomeInteger): auto = x[][i,j]
+
+#[
+template `[]=`*[T](x: Color, i:T, y: typed): untyped =
+  when T.isWrapper and T isnot Color:
+    x[][asScalar(i)] = y
+  else:
+    x[][i] = y
+]#
+template `[]=`*[T](x: Color, i: T; y: auto) =
+  when T is Color2:
+    x[][i[]] = y
+  elif T.isWrapper:
+    #indexed(x, i)
+    var tColorBracket = asColor(x[][i])
+    tColorBracket := y
+  else:
+    x[][i] = y
+template `[]=`*(x: Color, i,j: SomeInteger, y: auto) =
+  x[][i,j] = y
+
+# forward from value to value
+template forwardVV(f: untyped) {.dirty.} =
+  template f*(x: Color): auto =
+    mixin f
+    f(x[])
+# forward from value to type
+#template forwardVT(f: untyped) {.dirty.} =
+#  template f*[T](x: Color[T]): untyped =
+#    mixin f
+#    f(type T)
+# forward from type to type
+template forwardTT(f: untyped) {.dirty.} =
+  template f*[T](x: typedesc[Color[T]]): auto =
+    mixin f
+    f(type T)
+# forward from type to type and wrap
+template forwardTTW(f: untyped) {.dirty.} =
+  template f*[T](x: typedesc[Color[T]]): auto =
+    mixin f
+    Color[f(type T)]
+
+forwardVV(len)
+forwardVV(nrows)
+forwardVV(ncols)
+forwardVV(nVectors)
+forwardVV(simdType)
+#forwardVV(simdLength)
+forwardVV(getNs)
+forwardVV(numNumbers)
+
+forwardTT(len)
+forwardTT(nrows)
+forwardTT(ncols)
+forwardTT(nVectors)
+forwardTT(simdType)
+forwardTT(simdLength)
+forwardTT(getNs)
+forwardTT(numberType)
+
+forwardTTW(toSingle)
+forwardTTW(toDouble)
+
+#template eval*[T](x: typedesc[Color[T]]): typedesc = asColor(eval(type T))
+template eval*[T:Color](x: typedesc[T]): typedesc = asColor(eval((type T)[]))
+
+template has*[T:Color](x: typedesc[T], y: typedesc): bool =
+  mixin has
+  when y is Color: true
+  else: has(T.type[], y)
+
+template row*(x: Color, i: untyped): untyped =
+  mixin row
+  asColor(row(x[],i))
+template setRow*(r: Color; x: Color2; i: untyped): untyped =
+  setRow(r[], x[], i)
+
+template getNc*[T](x: Color[T]): untyped =
+  when T is Mat1:
+    x[].nrows
+  elif T is Vec1:
+    x[].len
+  else:
+    static:
+      echo "error: unknown Nc"
+      echo x.repr
+      echo type(x).name
+      qexExit 1
+    0
+
+template binDDRet(fn,wr,T1,T2) =
+  template fn*(x: T1, y: T2): untyped =
+    wr(fn(x[], y[]))
+  #template fn*(x: T1, y: T2): auto =
+  #  var tmp {.noInit.}: wr(evalType(fn(x[],y[]))
+  #  wr(fn(x[], y[]))
+
+#binDDRet(`+`, asColor, Color, Color2)
+binDDRet(`-`, asColor, Color, Color2)
+#binDDRet(`*`, asColor, Color, Color2)
+binDDRet(`/`, asColor, Color, Color2)
+
+setBinop(`+`, add, Color, Color2, asColor(evalType(x[]+y[])))
+setBinop(`*`, mul, Color, Color2, asColor(evalType(x[]*y[])))
+
+template load1*(x: Color): untyped = asColor(load1(x[]))
+template `-`*(x: Color): untyped = asColor(-(x[]))
+template assign*(r: var Color, x: SomeNumber) =
+  assign(r[], x)
+template assign*(r: var Color, x: AsComplex) =
+  assign(r[], x)
+template assign*(r: var Color, x: Color2) =
+  assign(r[], x[])
+template `:=`*(r: var Color, x: SomeNumber) =
+  r[] := x
+template `:=`*(r: var Color, x: AsComplex) =
+  r[] := x
+template `:=`*(r: var Color, x: Color2) =
+  r[] := x[]
+template `+=`*(r: var Color, x: Color2) =
+  r[] += x[]
+template `-=`*(r: var Color, x: Color2) =
+  r[] -= x[]
+template `*=`*(r: var Color, x: SomeNumber) =
+  r[] *= x
+template iadd*(r: var Color, x: SomeNumber) =
+  iadd(r[], x)
+template iadd*(r: var Color, x: AsComplex) =
+  iadd(r[], x)
+template iadd*(r: var Color, x: Color2) =
+  iadd(r[], x[])
+template isub*(r: var Color, x: SomeNumber) =
+  isub(r[], x)
+template isub*(r: var Color, x: Color2) =
+  isub(r[], x[])
+template imul*(r: var Color, x: SomeNumber) =
+  imul(r[], x)
+template imadd*(r: var Color, x: Color2, y: Color3) =
+  imadd(r[], x[], y[])
+template imadd*(r: var Color, x: AsComplex, y: Color3) =
+  imadd(r[], x, y[])
+template imsub*(r: var Color, x: Color2, y: Color3) =
+  imsub(r[], x[], y[])
+template peqOuter*(r: var Color, x: Color2, y: Color3) =
+  peqOuter(r[], x[], y[])
+template meqOuter*(r: var Color, x: Color2, y: Color3) =
+  meqOuter(r[], x[], y[])
+template `+`*(r: Color, x: SomeNumber): untyped =
+  asColor(r[] + x)
+template `+`*(r: SomeNumber, x: Color): untyped =
+  asColor(r + x[])
+template `-`*(r: Color, x: SomeNumber): untyped =
+  asColor(r[] - x)
+template `+`*(r: Color, x: AsComplex): untyped =
+  asColor(r[] + x)
+template `-`*(r: Color, x: AsComplex): untyped =
+  asColor(r[] - x)
+template add*(r: var Color, x: Color2, y: Color3) =
+  add(r[], x[], y[])
+template sub*(r: var Color, x: Color2, y: Color3) =
+  sub(r[], x[], y[])
+template `*`*(x: Color, y: SomeNumber): untyped =
+  asColor(x[] * y)
+#template `*`*(x: SomeNumber, y: Color): auto =
+  #asColor(x * y[])
+template `*`*[X:SomeNumber,Y:Color](x: X, y: Y): auto =
+  #static: echo "SomeNumber * Color"
+  #var tmp {.noInit.}: asColor(type(X)*type(Y)[])
+  var tmp {.noInit.}: asColor(evalType(x*y[]))
+  #static: echo $type(tmp)
+  mul(tmp[], x, y[])
+  tmp
+template `*`*(x: Simd, y: Color2): untyped =
+  asColor(x * y[])
+#template `*`*(x: AsReal, y: Color2): untyped =
+#  asColor(x * y[])
+template `*`*[X:AsReal,Y:Color](x: X, y: Y): auto =
+  #var tmp {.noInit.}: asColor(type(X)*type(Y)[])
+  var tmp {.noInit.}: asColor(evalType(x*y[]))
+  mul(tmp[], x, y[])
+  tmp
+template `*`*(x: AsImag, y: Color2): untyped =
+  asColor(x * y[])
+template `*`*(x: AsComplex, y: Color2): untyped =
+  asColor(x * y[])
+template `*`*(x: Color, y: AsComplex): untyped =
+  asColor(x[] * y)
+#template mul*(x: SomeNumber, y: Color2): untyped =
+#  asColor(`*`(x, y[]))
+template mul*[X:SomeNumber,Y:Color](x: X, y: Y): auto =
+  #var tmp {.noInit.}: asColor(type(X)*type(Y)[])
+  var tmp {.noInit.}: asColor(evalType(x*y[]))
+  mul(tmp[], x, y[])
+  tmp
+template mul*(x: Color, y: Color2): untyped =
+  asColor(`*`(x[], y[]))
+template mul*(r: var Color, x: Color2, y: Color3) =
+  mul(r[], x[], y[])
+template mul*(r: var Color, x: SomeNumber, y: Color3) =
+  mul(r[], x, y[])
+template mul*(r: var Color, x: Color2, y: SomeNumber) =
+  mul(r[], x[], y)
+template mul*(r: var Color, x: AsComplex, y: Color3) =
+  mul(r[], x, y[])
+template mul*(x: AsComplex, y: Color2): untyped =
+  asColor(mul(x, y[]))
+template random*(x: var Color) =
+  gaussian(x[], r)
+template gaussian*(x: Color, r: untyped) =
+  gaussian(x[], r)
+template uniform*(x: var Color, r: var untyped) =
+  uniform(x[], r)
+template z4*(x: var Color, r: var untyped) =
+  z4(x[], r)
+template z2*(x: var Color, r: var untyped) =
+  z2(x[], r)
+template u1*(x: var Color, r: var untyped) =
+  u1(x[], r)
+template projectU*(r: var Color) =
+  projectU(r[])
+template projectU*(r: var Color, x: Color2) =
+  projectU(r[], x[])
+template projectUderiv*(r: var Color, u: Color2, x: Color3, chain: Color4) =
+  projectUderiv(r[], u[], x[], chain[])
+template projectUderiv*(r: var Color, x: Color3, chain: Color4) =
+  projectUderiv(r[], x[], chain[])
+template projectSU*(r: var Color) =
+  projectSU(r[])
+template projectSU*(r: var Color, x: Color2) =
+  projectSU(r[], x[])
+template projectTAH*(r: var Color) =
+  projectTAH(r[])
+template projectTAH*(r: var Color, x: Color2) =
+  projectTAH(r[], x[])
+template checkU*(x: Color):untyped = checkU(x[])
+template checkSU*(x: Color):untyped = checkSU(x[])
+template norm2*(x: Color): untyped = norm2(x[])
+template norm2*(r: var auto, x: Color): untyped = norm2(r, x[])
+template inorm2*(r: var auto, x: Color2) = inorm2(r, x[])
+template dot*(x: Color, y: Color2): untyped =
+  dot(x[], y[])
+template idot*(r: var auto, x: Color2, y: Color3) = idot(r, x[], y[])
+template redot*(x: Color, y: Color2): untyped =
+  redot(x[], y[])
+template trace*(x: Color): untyped = trace(x[])
+template simdSum*(x: Color): untyped = asColor(simdSum(x[]))
+template re*(x: Color): untyped = asColor(re(x[]))
+template im*(x: Color): untyped = asColor(im(x[]))
+template exp*(x: Color): untyped = asColor(exp(x[]))
+template expDeriv*(x: Color, c: Color2): untyped = asColor(expDeriv(x[], c[]))
+template ln*(x: Color): untyped = asColor(ln(x[]))
diff --git a/src/physics/colorTensor.nim b/src/physics/colorTensor.nim
new file mode 100644
index 0000000..e9d8bfa
--- /dev/null
+++ b/src/physics/colorTensor.nim
@@ -0,0 +1,57 @@
+import tensorwrap
+export tensorwrap
+
+type
+  ColorType* = object
+  Color*[T] = TensorObj[ColorType,T]
+  Color2*[T] = Color[T]
+  Color3*[T] = Color[T]
+  Color4*[T] = Color[T]
+
+#[
+template getNc*[K,T](x: SomeTensor[K,T]): int =
+  when K is ColorType:
+    when T is Mat1:
+      x[].nrows
+    elif T is Vec1:
+      x[].len
+    else:
+      static:
+        echo "error: unknown Nc"
+        echo x.repr
+        echo type(x).name
+        qexExit 1
+      0
+  else:
+    getNc(x[])
+]#
+
+template getNc*[K,T](x: typedesc[SomeTensor[K,T]]): int =
+  getNc(x[])
+template getNc*[K,T](x: SomeTensor[K,T]): int =
+  getNc(x[])
+
+template getNc*[T](x: typedesc[SomeTensor[ColorType,T]]): int =
+  when T is Mat1:
+    x[].nrows
+  elif T is Vec1:
+    x[].len
+  else:
+    static:
+      echo "error: unknown Nc"
+      echo x.repr
+      echo type(x).name
+      qexExit 1
+    0
+template getNc*[T](x: SomeTensor[ColorType,T]): int =
+  when T is Mat1:
+    x[].nrows
+  elif T is Vec1:
+    x[].len
+  else:
+    static:
+      echo "error: unknown Nc"
+      echo x.repr
+      echo type(x).name
+      qexExit 1
+    0
diff --git a/src/physics/tensorwrap.nim b/src/physics/tensorwrap.nim
new file mode 100644
index 0000000..d4598e1
--- /dev/null
+++ b/src/physics/tensorwrap.nim
@@ -0,0 +1,347 @@
+import base/basicOps
+import base/wrapperTypes  # for getPtr
+export wrapperTypes
+import maths
+import simd/simdWrap
+getOptimPragmas()
+
+type
+  TensorObj*[K,T] = object
+    obj*: T
+  SomeTensor*[K,T] = TensorObj[K,T]
+  SomeTensor2*[K,T] = TensorObj[K,T]
+  SomeTensor3*[K,T] = TensorObj[K,T]
+  SomeTensor4*[K,T] = TensorObj[K,T]
+
+template tensorObj*[K,T](k: typedesc[K], t: typedesc[T]): typedesc =
+  TensorObj[K,T]
+template tensorObj*[K,T](k: typedesc[K], x: T): auto =
+  TensorObj[K,T](obj: x)
+
+template kind*[K,T](x: typedesc[SomeTensor[K,T]]): typedesc = K
+template kind*[K,T](x: SomeTensor[K,T]): typedesc = K
+
+template `[]`*[K,T](x: typedesc[SomeTensor[K,T]]): typedesc = T
+template `[]`*[K,T](x: TensorObj[K,T]): auto = x.obj
+#proc `[]`*[K,T](x: var TensorObj[K,T]): var T {.alwaysInline.} = x.obj
+
+template isWrapper*[K,T](x: typedesc[SomeTensor[K,T]]): bool = true
+template isWrapper*[K,T](x: SomeTensor[K,T]): bool = true
+template asWrapper*[K,T,Y](x: typedesc[SomeTensor[K,T]], y: typedesc[Y]): typedesc =
+  tensorObj(K, Y)
+template asWrapper*[K,T](x: typedesc[SomeTensor[K,T]], y: auto): auto =
+  tensorObj(K, y)
+template asWrapper*[K,T](x: SomeTensor[K,T], y: auto): auto =
+  tensorObj(K, y)
+
+template eval*[K,T](x: typedesc[TensorObj[K,T]]): typedesc =
+  tensorObj(K,eval(T))
+
+template has*[K,T,Y](x: typedesc[SomeTensor[K,T]], y: typedesc[Y]): bool =
+  mixin has
+  when Y is SomeTensor[K,auto]: true
+  else: has(T[], Y)
+
+template index*[K,T,I](x: typedesc[TensorObj[K,T]], i: typedesc[I]): typedesc =
+  when I is SomeTensor[K,auto]:
+    index(K, I[])
+  elif I.isWrapper:
+    tensorObj(K, index(T, I))
+  else:
+    index(T.type, I.type)
+
+template index*[K,T](x: typedesc[TensorObj[K,T]], i,j: typedesc[int]): typedesc =
+  index(T.type, int, int)
+
+template `[]`*[K,T,I](x: TensorObj[K,T], i: I): auto =
+  when I is SomeTensor[K,auto]:
+    x[][i[]]
+  elif I.isWrapper:
+    #indexed(x, i)
+    var tTensorObjBracket = tensorObj(K, x[][i])
+    tTensorObjBracket
+  else:
+    x[][i]
+
+template `[]=`*[K,T,I](x: SomeTensor[K,T], i: I, y: auto) =
+  mixin `:=`
+  when I is SomeTensor[K,auto]:
+    x[][i[]] = y
+  else:
+    x[][i] = y
+
+template `[]`*[K,T,I](x: TensorObj[K,T], i,j: I): auto =
+  x[][i,j]
+template `[]=`*[K,T,I](x: TensorObj[K,T], i,j: I, y: auto) =
+  x[][i,j] = y
+
+# forward from value to value
+template forwardVV(f: untyped) {.dirty.} =
+  template f*(x: SomeTensor): auto =
+    mixin f
+    f(x[])
+# forward from value to type
+#template forwardVT(f: untyped) {.dirty.} =
+#  template f*[T](x: SomeTensor[T]): auto =
+#    mixin f
+#    f(type T)
+# forward from type to type
+template forwardTT(f: untyped) {.dirty.} =
+  template f*[K,T](x: typedesc[SomeTensor[K,T]]): auto =
+    mixin f
+    f(type T)
+# forward from type to type and wrap
+template forwardTTW(f: untyped) {.dirty.} =
+  template f*[K,T](x: typedesc[TensorObj[K,T]]): auto =
+    mixin f
+    tensorObj(K, f(type T))
+
+forwardVV(len)
+forwardVV(nrows)
+forwardVV(ncols)
+forwardVV(nVectors)
+forwardVV(simdType)
+#forwardVV(simdLength)
+forwardVV(getNs)
+forwardVV(numberType)
+forwardVV(numNumbers)
+
+forwardTT(len)
+forwardTT(nrows)
+forwardTT(ncols)
+forwardTT(nVectors)
+forwardTT(simdType)
+forwardTT(simdLength)
+forwardTT(getNs)
+forwardTT(numberType)
+
+forwardTTW(toSingle)
+forwardTTW(toDouble)
+
+template row*(x: SomeTensor, i: auto): auto =
+  mixin row
+  tensorObj(x.kind, row(x[],i))
+template setRow*[K,T1,T2](r: SomeTensor[K,T1]; x: SomeTensor2[K,T2]; i: auto): auto =
+  setRow(r[], x[], i)
+
+#** Unary assignment functions
+template setUnaryAssignT(f) {.dirty.} =
+  template f*[R:SomeTensor,X:SomeTensor2](r: var R, x: X) =
+    when R.kind is X.kind:
+      f(r[], x[])
+    else:  # assume R has X
+      f(r[], x)
+template setUnaryAssignX(f,X) {.dirty.} =
+  template f*[R:SomeTensor](r: var R, x: X) =
+    f(r[], x)
+
+setUnaryAssignT(assign)
+setUnaryAssignX(assign, SomeNumber)
+setUnaryAssignX(assign, AsComplex)
+
+setUnaryAssignT(`:=`)
+setUnaryAssignX(`:=`, SomeNumber)
+setUnaryAssignX(`:=`, AsComplex)
+
+setUnaryAssignT(`+=`)
+setUnaryAssignX(`+=`, SomeNumber)
+setUnaryAssignX(`+=`, AsComplex)
+
+setUnaryAssignT(`-=`)
+setUnaryAssignX(`-=`, SomeNumber)
+setUnaryAssignX(`-=`, AsComplex)
+
+setUnaryAssignT(`*=`)
+setUnaryAssignX(`*=`, SomeNumber)
+setUnaryAssignX(`*=`, AsComplex)
+
+setUnaryAssignT(iadd)
+setUnaryAssignX(iadd, SomeNumber)
+setUnaryAssignX(iadd, AsComplex)
+
+setUnaryAssignT(isub)
+setUnaryAssignX(isub, SomeNumber)
+setUnaryAssignX(isub, AsComplex)
+
+setUnaryAssignT(imul)
+setUnaryAssignX(imul, SomeNumber)
+setUnaryAssignX(imul, AsComplex)
+
+#** Unary functions
+template setUnopT(f, g) {.dirty.} =
+  template f*(x: typedesc[SomeTensor]): typedesc =
+    tensorObj(x.kind, f(x[]))
+  template f*(x: SomeTensor): auto =
+    tensorObj(x.kind, f(x[]))
+
+setUnopT(load1, load1)
+setUnopT(`-`, neg)
+setUnopT(re, re)
+setUnopT(im, im)
+setUnopT(simdSum, simdSum)
+
+#** Binary assignment functions
+
+template setBinAssignTT(f) {.dirty.} =
+  template f*[R:SomeTensor,X:SomeTensor2,Y:SomeTensor3](r: var R, x: X, y: Y) =
+    when R.kind is X.kind:
+      when R.kind is Y.kind:
+        f(r[], x[], y[])
+      else:
+        f(r[], x[], y)
+    else:
+      when R.kind is Y.kind:
+        f(r[], x, y[])
+      else:
+        f(r[], x, y)
+template setBinAssignTX(f,Y) {.dirty.} =
+  template f*[R:SomeTensor,X:SomeTensor2](r: var R, x: X, y: Y) =
+    when R.kind is X.kind:
+      f(r[], x[], y)
+    else:
+      f(r[], x, y)
+template setBinAssignXT(f,X) {.dirty.} =
+  template f*[R:SomeTensor,Y:SomeTensor3](r: var R, x: X, y: Y) =
+    when R.kind is Y.kind:
+      f(r[], x, y[])
+    else:
+      f(r[], x, y)
+
+setBinAssignTT(add)
+setBinAssignTT(sub)
+
+setBinAssignTT(mul)
+setBinAssignXT(mul, SomeNumber)
+
+setBinAssignTT(imadd)
+setBinAssignXT(imadd, AsComplex)
+
+setBinAssignTT(imsub)
+
+#[
+template mul*(r: var SomeTensor, x: SomeTensor2, y: SomeNumber) =
+  mul(r[], x[], y)
+template mul*(r: var SomeTensor, x: AsComplex, y: SomeTensor3) =
+  mul(r[], x, y[])
+template peqOuter*(r: var SomeTensor, x: SomeTensor2, y: SomeTensor3) =
+  peqOuter(r[], x[], y[])
+template meqOuter*(r: var SomeTensor, x: SomeTensor2, y: SomeTensor3) =
+  meqOuter(r[], x[], y[])
+]#
+
+
+#** Binary functions
+
+template setBinopTT(f, g) {.dirty.} =
+  template f*(x: typedesc[SomeTensor], y: typedesc[SomeTensor2]): typedesc =
+    when x.kind is y.kind:
+      tensorObj(x.kind, f(x[],y[]))
+    elif x.kind.has y.kind:
+      f(x[], y)
+    else:  # assume y.kind has x.kind
+      f(x, y[])
+  template f*(x: SomeTensor, y: SomeTensor2): auto =
+    when x.kind is y.kind:
+      tensorObj(x.kind, f(x[],y[]))
+    elif x.kind.has y.kind:
+      f(x[], y)
+    else:  # assume y.kind has x.kind
+      f(x, y[])
+  #template f*(x: SomeTensor, y: SomeTensor2): auto =
+  #  var r {.noInit.}: f(x.type, y.type)
+  #  g(r, x, y)
+  #  r
+template setBinopTX(f,g,X) {.dirty.} =
+  template f*(x: typedesc[SomeTensor], y: typedesc[X]): typedesc =
+    tensorObj(x.kind, f(x[], y))
+  template f*(x: SomeTensor, y: X): auto =
+    tensorObj(x.kind, f(x[], y))
+  #template f*(x: SomeTensor, y: X): auto =
+  #  var r {.noInit.}: f(x.type, X)
+  #  g(r, x, y)
+  #  r
+template setBinopXT(f,g,X) {.dirty.} =
+  template f*(x: typedesc[X], y: typedesc[SomeTensor]): typedesc =
+    tensorObj(y.kind, f(x, y[]))
+  template f*(x: X, y: SomeTensor): auto =
+    tensorObj(y.kind, f(x, y[]))
+  #template f*(x: X, y: SomeTensor): auto =
+  #  var r {.noInit.}: f(X, y.type)
+  #  g(r, x, y)
+  #  r
+
+setBinopTT(`+`, add)
+setBinopTX(`+`, add, SomeNumber)
+setBinopXT(`+`, add, SomeNumber)
+
+setBinopTT(`-`, sub)
+setBinopTX(`-`, sub, SomeNumber)
+setBinopXT(`-`, sub, SomeNumber)
+
+setBinopTT(`*`, mul)
+setBinopTX(`*`, mul, SomeNumber)
+setBinopXT(`*`, mul, SomeNumber)
+setBinopTX(`*`, mul, Simd)
+setBinopXT(`*`, mul, Simd)
+setBinopXT(`*`, mul, AsReal)
+setBinopTX(`*`, mul, AsReal)
+setBinopXT(`*`, mul, AsImag)
+setBinopTX(`*`, mul, AsImag)
+setBinopXT(`*`, mul, AsComplex)
+setBinopTX(`*`, mul, AsComplex)
+
+#[
+#template mul*(x: SomeNumber, y: SomeTensor2): auto =
+#  asSomeTensor(`*`(x, y[]))
+template mul*[X:SomeNumber,Y:SomeTensor](x: X, y: Y): auto =
+  #var tmp {.noInit.}: asSomeTensor(type(X)*type(Y)[])
+  var tmp {.noInit.}: asSomeTensor(evalType(x*y[]))
+  mul(tmp[], x, y[])
+  tmp
+template mul*(x: SomeTensor, y: SomeTensor2): auto =
+  asSomeTensor(`*`(x[], y[]))
+template mul*(x: AsComplex, y: SomeTensor2): auto =
+  asSomeTensor(mul(x, y[]))
+template random*(x: var SomeTensor) =
+  gaussian(x[], r)
+]#
+setUnaryAssignX(gaussian, auto)
+setUnaryAssignX(uniform, auto)
+setUnaryAssignX(z2, auto)
+setUnaryAssignX(z4, auto)
+setUnaryAssignX(u1, auto)
+
+template projectU*(r: var SomeTensor) =
+  projectU(r[])
+template projectU*(r: var SomeTensor, x: SomeTensor2) =
+  projectU(r[], x[])
+template projectUderiv*(r: var SomeTensor, u: SomeTensor2, x: SomeTensor3, chain: SomeTensor4) =
+  projectUderiv(r[], u[], x[], chain[])
+template projectUderiv*(r: var SomeTensor, x: SomeTensor3, chain: SomeTensor4) =
+  projectUderiv(r[], x[], chain[])
+template projectSU*(r: var SomeTensor) =
+  projectSU(r[])
+template projectSU*(r: var SomeTensor, x: SomeTensor) =
+  projectSU(r[], x[])
+template projectTAH*(r: var SomeTensor) =
+  projectTAH(r[])
+template projectTAH*(r: var SomeTensor, x: SomeTensor2) =
+  projectTAH(r[], x[])
+template checkU*(x: SomeTensor): auto = checkU(x[])
+template checkSU*(x: SomeTensor): auto = checkSU(x[])
+template norm2*(x: SomeTensor): auto = norm2(x[])
+template norm2*(r: var auto, x: SomeTensor) = norm2(r, x[])
+template inorm2*(r: var auto, x: SomeTensor) = inorm2(r, x[])
+template dot*(x: SomeTensor, y: SomeTensor2): auto =
+  dot(x[], y[])
+template idot*(r: var auto, x: SomeTensor2, y: SomeTensor3) =
+  idot(r, x[], y[])
+template redot*(x: SomeTensor, y: SomeTensor2): auto =
+  redot(x[], y[])
+template trace*(x: SomeTensor): auto = trace(x[])
+template exp*(x: SomeTensor): auto =
+  tensorObj(x.kind, exp(x[]))
+template expDeriv*(x: SomeTensor, c: SomeTensor2): auto =
+  tensorObj(x.kind, expDeriv(x[], c[]))
+template ln*(x: SomeTensor): auto =
+  tensorObj(x.kind, ln(x[]))