Adding solveMatrix.py from master

pythonScripts/solveMatrix.py

@@ -0,0 +1,786 @@
+import numpy
+
+def integrate(t,y):
+#integrate function y (N,n) or (N,) or a scalar tabulated at t (N,) values
+
+#compose time and value shift arrays
+   t0=[x for x in t]
+   t1=[x for x in t]
+#t0.insert(0,0)
+   t0=t0[:-1]
+   t1=t1[1:]
+   t0=numpy.array(t0)
+   t1=numpy.array(t1)
+   try:
+      n=y.shape[1]
+      y0=numpy.matrix.copy(y)
+      y1=numpy.matrix.copy(y)
+   except IndexError:
+      try:
+         n=len(y)
+         y0=numpy.outer(numpy.ones(t.shape),y)
+         y1=numpy.outer(numpy.ones(t.shape),y)
+      except TypeError:
+         n=1
+         y0=y*numpy.outer(numpy.ones(t.shape),numpy.ones(n))
+         y1=y*numpy.outer(numpy.ones(t.shape),numpy.ones(n))
+
+   y0=numpy.delete(y0,-1,axis=0)
+   y1=numpy.delete(y1,0,axis=0)
+   
+
+#integral (trapezoid) updates
+#(N,n)
+   dt=numpy.outer(t1-t0,numpy.ones(n))
+   #print('{} {}'.format(dt.shape,(y0+y1).shape))
+   d=0.5*dt*(y1+y0)
+
+#use cumsum to compute integrals
+#(N-1,n)
+   return numpy.cumsum(d,axis=0)
+
+
+def integrateFull(t,Y,axis=1):
+    #calculate integrals of Y along axis, assume value i of each vector along axis to be calculated at ti
+
+#the formula for trapezoid is 
+#2*I=y[n]*t[n]-y[0]*t[0]
+#+[y[0]*t[1]+y[1]*t[2]+...+y[n-2]*t[n-1]]+
+#-[y[1]*t[0]+y[2]*t[1]+...+y[n-1]*t[n-2]]=
+#y*(t1-t0)
+#where
+#t0=[t[0],t[0],...,t[n-3],t[n-2]]
+#t1=[t[1],t[2],...,t[n-1],t[n-1]]
+    t0=numpy.concatenate((t[0:1],t[:-1]))
+    t1=numpy.concatenate((t[1:],t[-1:]))
+    #this works for arbitrary dimension of Y, so some Einstein magic has to be applied
+    v='abcdef'
+    n=len(Y.shape)
+    einsumStr=v[:n]+','+v[axis]+'->'+v[:n]
+    print('einsum [{}] {} {}'.format(einsumStr,Y.shape,t0.shape))
+    T=numpy.einsum(einsumStr,numpy.ones(Y.shape),t1-t0)
+    return 0.5*numpy.sum(Y*T,axis)
+
+
+def skipRedundant(model):
+   redundant=[r for r in model.scaled]
+   redundant.append('total')
+   lut=model.lut
+   v_sorted=numpy.sort([lut[x] for x in redundant])[::-1]
+#n=len(Q['lut'])
+#sort by value
+   sMap={x:lut[x] for x in lut if x not in redundant}
+#sMap=dict(sorted(redundantMap.items(),key=lambda item:item[1]))
+   updatedMap={list(sMap.keys())[i]:i for i in range(len(sMap))}
+#print(sMap)
+#print(updatedMap)
+#return v_sorted
+   return updatedMap
+
+
+
+def removeScaled(model,M,mode='matrix'):
+   originalMap=model.lut
+   updatedMap=skipRedundant(model)
+   n=len(updatedMap)
+#print(skip)
+   if mode=='matrix':
+      A=numpy.zeros((n,n))
+      for c in updatedMap:
+         i=updatedMap[c]
+         i0=originalMap[c]
+         for c1 in updatedMap:
+            j=updatedMap[c1]
+            j0=originalMap[c1]
+            A[i,j]=M[i0,j0]
+      return A
+   if mode=='columnsOnly':
+      A=numpy.zeros((M.shape[0],n))
+      for c in updatedMap:
+         i=updatedMap[c]
+         i0=originalMap[c]
+         A[:,i]=M[:,i0]
+      return A
+   if mode=='vector':
+      v=numpy.zeros(n)
+      for c in updatedMap:
+         i=updatedMap[c]
+         i0=originalMap[c]
+         v[i]=M[i0]
+      return v
+   return numpy.zeros(0)
+
+def getSE(s1):
+   s2=numpy.multiply(s1,s1)
+   return numpy.sqrt(numpy.dot(s2,numpy.ones(s1.shape[2])))
+
+def removeImVector(v,eps):
+   for i in range(len(v)):
+      if numpy.abs(numpy.imag(v[i]))<eps*numpy.abs(numpy.real(v[i])):
+         v[i]=numpy.real(v[i])
+   return v
+
+def removeImMatrix(A,eps):
+   f=A.ravel()  
+   f=removeImVector(f,eps)
+   return f.reshape(A.shape)
+ 
+def solveMatrix(model,tmax,nt=201,t0=0,y0=numpy.array([]), sIn=numpy.array([]),method='full'):
+   if method=='full':
+      obj=_solveMatrix(model,tmax,nt,t0,y0,sIn)
+   if method=='separately':
+      obj=_solveMatrixSeparately(model,tmax,nt,t0,y0,sIn)
+
+   return obj['t'],obj['sol'],obj['se'],obj['s1']
+
+def _solveMatrix(model,tmax,nt=201,t0=0,y0=numpy.array([]), sIn=numpy.array([])):
+   t=numpy.linspace(0,tmax-t0,nt)
+
+#solve system Q at time points t for initial values y0
+#SS is gradient of system nxn matrix M with respect to m parameters p 
+   SS=model.fSS(0)
+#large (M+1)*N by (M+1)*N matrix
+   m=len(model.lutSE)
+   A=removeScaled(model,model.M(0))
+#adjusted n where only non-terminal(scaled) compartments are counted
+   n=A.shape[0]
+   SM=numpy.zeros(((m+1)*n,(m+1)*n))
+   for k in range(m):
+#A is on diagonal
+      SM[k*n:(k+1)*n,k*n:(k+1)*n]=A
+#on the far right, gradient matrix SS sits
+      SM[k*n:(k+1)*n,-n:]=removeScaled(model,SS[k])
+#for original solutions, matrix A is used
+      SM[-n:,-n:]=A
+#matrix SM is now of full rank and can be diagonalized
+#print('{}/{}'.format(numpy.linalg.matrix_rank(SM),(m+1)*n))
+
+#RHS
+#SU has shape mxn
+#we have to remove scaled columns
+#Add plain rhs
+   SU=removeScaled(model,model.Su(0),mode='columnsOnly').flatten()
+   u=removeScaled(model,model.u(0),mode='vector')
+   print('{} {}'.format(SU.shape,u.shape))
+   SU=numpy.concatenate((SU,u))
+   print(SU.shape)
+#diagonalize matrix
+   lam,P=numpy.linalg.eig(SM)
+#no effect down to eps=1e-2
+   #eps=1e-4
+   #print(f'Using eps={eps}')
+   #lam=removeImVector(lam,eps)
+   #P=removeImMatrix(P,eps)
+   P1=numpy.linalg.inv(P)
+   #P1=removeImMatrix(P1,eps)
+   D=numpy.diag(lam)
+
+#verify
+#print(numpy.max(P @ D @ P1-M))
+
+#shift input to diagonal system
+   v=P1 @ SU
+   #also no effect
+   #maxV=numpy.max(numpy.abs(v))
+   #v=v*(numpy.abs(v)>eps*maxV)
+   #v=removeImVector(v,eps)
+
+   #print(f'v={v}')
+   #print(f'lam={lam}')
+
+#set initial values if missing
+   if not y0.size:
+      y0p=numpy.zeros(n)
+   else:
+#should move this to diagonal space
+      y0p=removeScaled(model,y0,mode='vector')
+
+   if not sIn.size:
+      s0=numpy.zeros(n*m)
+   else:
+#sIn is n by m matrix, so transpose before reshaping
+      s0=numpy.transpose(sIn)
+      s0=removeScaled(model,s0,mode='columnsOnly')
+      s0=s0.flatten()
+   print('s0 {} y0p {}'.format(s0.shape,y0p.shape))
+   S0=numpy.concatenate((s0,y0p))
+   S0= P1 @ S0   
+
+#present results as a (n*(m+1) x N array)
+#n number of variables/compartments
+#N number of time points
+   W=numpy.outer(lam,t)
+   tE=numpy.ones(t.shape)
+   V0=numpy.outer(S0,tE)
+   V=numpy.outer(v/lam,tE)
+   y1=(V0+V)*numpy.exp(W)-V
+
+#convert back to true system
+   y=P @ y1
+
+#reinsert scaled into array
+#start with y of shape (m+1)*n by N where N is number of time-points and n are only non-scaled compartments
+#should end up with y of shape (m+1)*n by N, where n is now number of ALL compartments
+   originalMap=model.lut
+   updatedMap=skipRedundant(model)
+   n0=len(originalMap)
+   n=len(updatedMap)
+   N=y.shape[1]
+   sOut=numpy.zeros((N,n0,m+1))
+   for k in range(m+1):
+      for c in updatedMap:
+         i=updatedMap[c]
+         i0=originalMap[c]
+         sOut[:,i0,k]=y[k*n+i,:]
+
+#print('Shape: {}, n={} N={} m={}'.format(yout.shape,n0,N,m))
+#equivalent of Sout (N,n,m)
+#return numpy.transpose(yout).reshape(N,n0,m+1)
+   sOut=calculateScaled(model,t,sOut,sIn,y0)
+   sol=sOut[:,:,-1]
+   s1=sOut[:,:,:-1]
+   se=getSE(s1)
+   return {'t':t+t0,'sol':sol,'se':se,'s1':s1,'P':P,'D':D,'SM':SM}
+
+
+
+
+def calculateScaled(model,t,sOut,sIn=numpy.zeros(0),y0=numpy.zeros(0)):
+
+#update sOut w/ calculated scaled values
+#sIn,yIn are initial values
+#sIn is of shape (n,m)
+#yIn is (n,)
+#sOut is (N,n,m+1)
+   lutSE=model.lutSE
+   lut=model.lut
+
+#if reading a model, keep the read lut
+#lutSE=Q['lutSE']
+#lut=Q['lut']
+
+   m=len(lutSE)
+   n=len(lut)
+   N=len(t)
+   
+   if not sIn.size:
+      sIn=numpy.zeros((n,m))
+
+   if not y0.size:
+      y0=numpy.zeros(n)
+
+   #add column for initial values
+   sIn=numpy.c_[sIn,y0]
+
+
+#reshape to n*(m+1) by N
+   y=numpy.zeros((n*(m+1),N))
+
+   for k in range(m+1):
+      for i in range(n):
+         y[n*k+i,:]=sOut[:,i,k]
+
+#full version of SM
+
+   SS=model.fSS(0)
+   A=model.M(0)
+   SM=numpy.zeros(((m+1)*n,(m+1)*n))
+   for k in range(m):
+#A is on diagonal
+      SM[k*n:(k+1)*n,k*n:(k+1)*n]=A
+#on the far right, gradient matrix SS sits
+      SM[k*n:(k+1)*n,-n:]=model.SS[k]
+
+   SM[-n:,-n:]=A
+
+#full version of RHS
+   SU=model.Su(0).flatten()
+   u=model.u(0)
+#n*(m+1) by N
+   SU=numpy.outer(numpy.concatenate((SU,u)),numpy.ones(t.shape))
+
+#integral, n*(m+1) by N
+   fI=integrate(t,y)
+   fU=integrate(t,SU)
+
+#apply couplings
+   fI= SM @ fI
+
+#update values; scale compartments to total exposure
+   iT=lut['total']
+   fTotal=y0[iT]+fU[m*n+iT,:]
+   for k in range(m+1):
+      for c in model.scaled:
+         i=lut[c]
+#sOut[1:,i,k]=fI[k*n+i,:]+fU[k*n+i,:]
+         sOut[1:,i,k]=(sIn[i,k]*y0[iT]+fI[k*n+i,:]+fU[k*n+i,:])/fTotal
+         sOut[0,i,k]=sIn[i,k]
+#set total for solution only
+   sOut[1:,iT,-1]=fTotal
+   sOut[0,iT,-1]=y0[iT]
+   return sOut
+
+def _solveMatrixSeparately(model,tmax,nt=201,t0=0,y0=numpy.array([]), sIn=numpy.array([])):
+#sIn (n,m)
+#y0 (n,)
+   t=numpy.linspace(0,tmax-t0,nt)
+   n=len(model.lut)
+   m=len(model.lutSE)
+   N=len(t)
+   
+   if not y0.size:
+      y0=numpy.zeros(n)
+
+   iT=model.lut['total']
+   u0=y0[iT]
+   solver=solveUnitObj(model,t,u0)
+
+
+   if not sIn.size:
+      sIn=numpy.zeros([n,m])
+   
+  
+   sOut=numpy.zeros([N,n,m+1])
+#y is (N,n)
+   #y=solveUnit(model,t,y0,u0,model.u(0))   
+   y=solver.solve(y0,model.u(0))   
+   sOut[:,:,-1]=y
+#(N,n,1)
+   y1=numpy.expand_dims(y,axis=2)
+
+
+#SS (m,n,n)
+   SS=model.fSS(0)
+#SU (m,n)
+   SU=model.Su(0)  #m by n
+#u (n0,)
+   u=model.u(0) 
+   for p in model.lutSE:
+      k=model.lutSE[p]
+#SS[k] is (n,n)
+#y is (N,n)
+#rhs is (N,n)=(n,n) @ (N,n,1)
+      rhs=(SS[k] @ y1).squeeze()
+#yIn is (n,)
+      yIn=sIn[:,k]
+#return in (N,n)
+#      sOut[:,:,k]=solveUnit(model,t,yIn,u0,u,rhs)
+      #sOut[:,:,k]=solver.solve(sIn[:,k],SU[k],rhs)
+      sOut[:,:,k]=solver.solveS(SS[k],y1,y0,u=u,z0=sIn[:,k],r=SU[k])
+
+
+#      break
+      
+   sol=sOut[:,:,-1]
+   s1=sOut[:,:,:-1]
+   se=getSE(s1)
+   return {'t':t+t0,'sol':sol,'se':se,'s1':s1}
+
+
+class solveUnitObj:
+
+   def __init__(self,model,t,u0):
+
+      self.model=model
+      self.t=t
+      self.u0=u0
+#some dimensions
+      self.originalMap=model.lut
+      self.updatedMap=skipRedundant(model)
+      self.n0=len(self.originalMap)
+      self.n=len(self.updatedMap)
+ 
+      self.N=len(t)
+      
+      self.A=removeScaled(model,model.M(0))
+
+#diagonalize matrix
+      self.lam,self.P=numpy.linalg.eig(self.A)
+      self.P1=numpy.linalg.inv(self.P)
+      D=numpy.diag(self.lam)
+      self.D1=numpy.diag(1/self.lam)
+
+
+#next are all (N,n,n)
+      self.tE=numpy.ones(t.shape)
+      E=numpy.einsum('a,bc',self.tE,numpy.eye(self.n))
+      QT=numpy.einsum('a,bc',t,numpy.eye(self.n))
+      QTL=numpy.einsum('a,bc',t,D)
+      self.W=E*numpy.exp(QTL)
+      invL=numpy.einsum('a,bc',self.tE,self.D1)
+      self.Q=(self.W-E)*invL
+
+#we have exp(lam*(t-t')) as Green's function
+#construct matrix W1 (N,N,n,n) where the last part is diag(lam)
+#W[i,j,k,l]=delta(k,l)*exp(lam_k*(t_i-t'_j))
+#delta(i,j) is Kronecker's delta
+      T=numpy.outer(t,self.tE)-numpy.outer(self.tE,t)
+      T1=numpy.einsum('ab,cd',T,numpy.eye(self.n))
+      F=numpy.einsum('ab,cd',T>=0,numpy.eye(self.n))
+      L1=numpy.einsum('ab,cd',numpy.ones(T.shape),D)
+      self.W1=numpy.zeros(F.shape,dtype='complex')
+      self.W1[F>0]=numpy.exp(L1[F>0]*T1[F>0])
+
+#a matrix with lam_i in rows LR(i,j)=lam(i)
+#(N,n,n)
+      LR=numpy.einsum('a,b,c',self.tE,self.lam,numpy.ones(self.n))
+      LRT=numpy.einsum('a,b,c',t,self.lam,numpy.ones(self.n))
+#a matrix with lam_i in columns LC(i,j)=lam(j)
+#(N,n,n)
+      LC=numpy.einsum('a,b,c',self.tE,numpy.ones(self.n),self.lam)
+      LCT=numpy.einsum('a,b,c',t,numpy.ones(self.n),self.lam)
+
+#a matrix with exp(lam_i) in rows DR(i,j)=exp(lam(i))
+#(N,n,n)
+      DR=numpy.exp(LRT)
+#a matrix with exp(lam_j) in columns DC(i,j)=exp(lam(j))
+#(N,n,n)
+      DC=numpy.exp(LCT)
+
+#a diagonal matrix with t*exp(D*t) on diagonal
+#(N,n,n)
+      self.H=numpy.zeros(self.W.shape,dtype='complex')
+      self.H[E>0]=QT[E>0]*numpy.exp(QTL[E>0])
+#off diagonal is exp(lam(j))-exp(lam(i))/(lam(j)-lam(i))
+      self.H[E==0]=(DC[E==0]-DR[E==0])/(LC[E==0]-LR[E==0])
+      
+      
+#
+
+
+   def solve(self,y0=numpy.zeros(0),u=numpy.zeros(0),rhs=numpy.zeros(0)):
+
+#solve system Q at time points t for initial values y0
+#u is the time independent nohomogeneous part
+#rhs is tabulated non-homogeneous part (N,n)
+#u0 is initial accumulated xenobiote
+
+
+#use Green function theorem, since impulse response is simple
+#yGREEN=exp(lam*(t-t0)
+   
+#matrix solution works for reduced space where all
+#non-trivial components are reduced
+
+#typical variables are denoted with trailing p (for prime)
+#and have a smaller size (fewer compartments) than
+#full variables
+
+
+#get the matrix and remove trivial components
+#remove trivial components from non-homogeneous parts
+      if not u.size:
+         u=numpy.zeros(self.n0)
+
+      up=removeScaled(self.model,u,mode='vector')
+
+#remove trivial components from initial values
+      if not y0.size:
+         y0=numpy.zeros(u.shape)
+      y0p=removeScaled(self.model,y0,mode='vector')
+
+#overload w/ rhs
+#RHS has a time component
+#typically of shape n by N
+      if rhs.size:
+#to N by n shape
+#      rhs=rhs.transpose()
+#remove trivial components
+         rhsp=removeScaled(self.model,rhs,mode='columnsOnly')
+#to (N,n,1)
+         rhsp=numpy.expand_dims(rhsp,axis=2) 
+      else:
+         rhsp=numpy.zeros(0)
+   #print(numpy.linalg.matrix_rank(A))
+
+
+
+#(n,n) @ (N,n,n) @ (n,n) @ (n,)
+#y1 is (N,n)
+#sum works for (N,n)+(n,)
+   #print(up)
+   #print(P1 @ up)
+   #print( (W-E)*invL @ P1 @ up)
+      y= self.P @ self.W @ self.P1 @ y0p + self.P @ self.Q @ self.P1 @ up
+
+      if rhsp.size:
+
+#time dependent component
+#(a,b,c,e)=(a,b,c,d) @ (b,d,e)
+#(N,N,n,1)=(n,n) @ (N,N,n,n) @ (n,n) @ (N,n,1)
+#Y1 is (N,N,n,1)
+         Y1=numpy.real(self.P @ self.W1 @ self.P1 @ rhsp)
+#to apply Green's theorem, integral along axis=1 should be performed
+#(N,n,1)
+      
+         Y2=integrateFull(self.t,Y1,axis=1)
+#to remove last dimension, squeeze at the end
+#(N,n), same as y1
+         y+=Y2.squeeze()
+     
+#back to original system, add trivial compartments
+      yOut=numpy.zeros((self.N,self.n0))
+   
+      for c in self.updatedMap:
+         i=self.updatedMap[c]
+         i0=self.originalMap[c]
+         yOut[:,i0]=y[:,i]
+  
+      #return yOut
+#here we need full sized y0,u,rhs, so make sure they are not overwritten by reduced versions
+      return getScaled(self.model,self.t,yOut,self.u0,y0,u,rhs)
+
+   def solveS(self,S,y,y0=numpy.zeros(0),u=numpy.zeros(0),\
+      z0=numpy.zeros(0),r=numpy.zeros(0)):
+
+#solve gradient system (S) at time points t 
+#y (N,n,1) is solution of the system 
+#y0 (n,) initial values of solution  
+#u (n,)is the non-homogenous (but constant) part of equation
+#z0 (n,)initial values of gradients
+#r (n,) is the time independent gradient of nonhomogenous part
+
+#use Green function theorem, since impulse response is simple
+#yGREEN=exp(lam*(t-t0) and integrate it for exponentail solutions 
+   
+#matrix solution works for reduced space where all
+#non-trivial components are reduced
+
+#typical variables are denoted with trailing p (for prime)
+#and have a smaller size (fewer compartments) than
+#full variables
+
+
+#get the matrix and remove trivial components
+
+#remove trivial components from initial values
+      if not y0.size:
+         y0=numpy.zeros(self.n0)
+      y0p=removeScaled(self.model,y0,mode='vector')
+
+#remove trivial components from non-homogenuous part of base equation
+      if not u.size:
+         u=numpy.zeros(y0.shape)
+      up=removeScaled(self.model,u,mode='vector')
+
+#remove trivial components from gradient of non-homogenous part
+      if not r.size:
+         r=numpy.zeros(self.n0)
+      rp=removeScaled(self.model,r,mode='vector')
+
+#remove trivial components from initial values of gradient
+      if not z0.size:
+         z0=numpy.zeros(r.shape)
+      z0p=removeScaled(self.model,z0,mode='vector')
+
+
+      Sp=removeScaled(self.model,S)
+#Spp is nearly diagonal in (n,n)
+      Spp=self.P1 @ Sp @ self.P
+#converted to time space for H multiplication
+#(N,n,n)
+      HS= numpy.einsum('a,bc',self.tE,Spp) * self.H
+
+#z=P @ W @ P1 @ z0 + P @ Q @ P1 r + 
+# + P @ (S * H) @ P1 @ y0 + P @ (S * H) @ D1 @ P1 @ u +P @ S @ Q @ D1 @ P1 @ u
+#(n,n) @ (N,n,n) @ (n,n) @ (n,)
+#y1 is (N,n)
+#sum works for (N,n)+(n,)
+
+      z1 = self.P @ self.W @ self.P1 @ z0p 
+      z2 = self.P @ self.Q @ self.P1 @ rp
+
+#intermediate variable for speed-up
+#(n,n) @ (n,n) @(n,)
+      _f=self.D1 @ self.P1 @ up
+
+      #combine above for speed-up
+#(n,n) @ (N,n,n) @ ( (n,n) @ (n,) + (n,) )
+      z3 = self.P @ HS @ (self.P1 @ y0p + _f)
+#(n,n) @ (n,n) @ (N,n,n) @ (n,)
+      z4 = self.P @ self.Q @ Spp @ _f
+     
+      z=z1+z2+(z3-z4)
+
+#back to original system, add trivial compartments
+      sOut=numpy.zeros((self.N,self.n0))
+   
+      for c in self.updatedMap:
+         i=self.updatedMap[c]
+         i0=self.originalMap[c]
+         sOut[:,i0]=z[:,i]
+  
+      #return yOut
+#here we need full sized y0,u,rhs, so make sure they are not overwritten by reduced versions
+      #y = self.P @ self. W @ self.P1 @ y + self.P @ self.Q @ self.P1 @ u
+      #(N,n0) = (n0,n0) @ (N,n0,1)
+      rhs = (S @ y).squeeze()
+      return getScaled(self.model,self.t,sOut,self.u0,z0,r,rhs)
+
+
+ 
+
+def solveUnit(model,t,y0=numpy.zeros(0),u0=0,u=numpy.zeros(0),rhs=numpy.zeros(0)):
+#solve system Q at time points t for initial values y0
+#u is the time independent nohomogeneous part
+#rhs is tabulated non-homogeneous part (N,n)
+#u0 is initial accumulated xenobiote
+
+#use Green function theorem, since impulse response is simple
+#yGREEN=exp(lam*(t-t0)
+   
+#some dimensions
+   originalMap=model.lut
+   updatedMap=skipRedundant(model)
+   n0=len(originalMap)
+   n=len(updatedMap)
+   N=len(t)
+
+#matrix solution works for reduced space where all
+#non-trivial components are reduced
+
+#typical variables are denoted with trailing p (for prime)
+#and have a smaller size (fewer compartments) than
+#full variables
+
+#get the matrix and remove trivial components
+   A=removeScaled(model,model.M(0))
+#remove trivial components from non-homogeneous parts
+   if not u.size:
+      u=numpy.zeros(n0)
+
+   up=removeScaled(model,u,mode='vector')
+
+#remove trivial components from initial values
+   if not y0.size:
+      y0=numpy.zeros(u.shape)
+   y0p=removeScaled(model,y0,mode='vector')
+
+#overload w/ rhs
+#RHS has a time component
+#typically of shape n by N
+   if rhs.size:
+#to N by n shape
+#      rhs=rhs.transpose()
+#remove trivial components
+      rhsp=removeScaled(model,rhs,mode='columnsOnly')
+#to (N,n,1)
+      rhsp=numpy.expand_dims(rhsp,axis=2) 
+   else:
+      rhsp=numpy.zeros(0)
+   #print(numpy.linalg.matrix_rank(A))
+
+#diagonalize matrix
+   lam,P=numpy.linalg.eig(A)
+   P1=numpy.linalg.inv(P)
+   D=numpy.diag(lam)
+
+#next are all (N,n,n)
+   W=numpy.exp(numpy.einsum('a,bc',t,numpy.diag(lam)))
+   tE=numpy.ones(t.shape)
+   invL=numpy.einsum('a,bc',tE,numpy.diag(1/lam))
+   E=numpy.ones(W.shape)
+
+#(n,n) @ (N,n,n) @ (n,n) @ (n,)
+#y1 is (N,n)
+#sum works for (N,n)+(n,)
+   #print(up)
+   #print(P1 @ up)
+   #print( (W-E)*invL @ P1 @ up)
+   y= y0p + P @ ((W-E)*invL) @ P1 @ up
+
+   if rhsp.size:
+
+#time dependent component
+#t-t'
+#we have exp(lam*(t-t')) as Green's function
+#construct matrix (N,N,n,n) where the last part is diag(lam)
+#W[i,j,k,l]=delta(k,l)*exp(lam_k*(t_i-t'_j))
+#delta(i,j) is Kronecker's delta
+      T=numpy.outer(t,tE)-numpy.outer(tE,t)
+      T1=numpy.einsum('ab,cd',T,numpy.eye(n))
+      F=numpy.einsum('ab,cd',T>=0,numpy.eye(n))
+      L1=numpy.einsum('ab,cd',numpy.ones(T.shape),numpy.diag(lam))
+      W1=numpy.zeros(F.shape,dtype='complex')
+      W1[F>0]=numpy.exp(L1[F>0]*T1[F>0])
+#(a,b,c,e)=(a,b,c,d) @ (b,d,e)
+#(N,N,n,1)=(n,n) @ (N,N,n,n) @ (n,n) @ (N,n,1)
+#Y1 is (N,N,n,1)
+      Y1=P @ W1 @ P1 @ rhsp
+#to apply Green's theorem, integral along axis=1 should be performed
+#(N,n,1)
+      
+      Y2=integrateFull(t,Y1,axis=1)
+#to remove last dimension, squeeze at the end
+#(N,n), same as y1
+      y+=Y2.squeeze()
+     
+#back to original system, add trivial compartments
+   yOut=numpy.zeros((N,n0))
+   
+   for c in updatedMap:
+      i=updatedMap[c]
+      i0=originalMap[c]
+      yOut[:,i0]=y[:,i]
+   
+#here we need full sized y0,u,rhs, so make sure they are not overwritten by reduced versions
+   return getScaled(model,t,yOut,u0,y0,u,rhs)
+
+def getTotal(model,t,u0=0):
+   u=model.u(0)
+   iT=model.lut['total']
+   U=u[iT]
+#*numpy.ones(t.shape)
+   fU=numpy.zeros(t.shape)
+   fU[1:]=u0+integrate(t,U).squeeze()
+   fU[0]=u0
+   return fU
+
+
+def getScaled(model,t,y,u0=0,y0=numpy.zeros(0),u=numpy.zeros(0),rhs=numpy.zeros(0)):
+#update y with scaled containers, which were excluded from direct calculation to make M of full rank
+#u has shape (n,)
+#rhs has shape (N,n)
+#y is (N,n)
+   M=model.M(0)
+
+#integrate input
+#(N-1,n)
+   fI=integrate(t,y)
+
+#assemble non-homogenous part of equation
+   if not rhs.size:
+      rhs=numpy.zeros(y.shape[1])
+   if u.size:
+#rhs is either (n,) or (N,n) so adding (n,) should work
+      rhs+=u
+
+#(N-1,n) even if rhs is (n,)
+   fU=integrate(t,rhs)
+
+#apply coupling
+#(N-1,n)
+   fI=(M @ numpy.expand_dims(fI,axis=2)).squeeze()
+
+#find initial values
+#(n,)
+   if not y0.size:
+      y0=numpy.zeros(u.shape)
+
+#get starting value for fTotal
+#(N-1,)
+   fTotal=getTotal(model,t,u0)
+
+#update solutions with scaled compartments
+   for c in model.scaled:
+      i=model.lut[c]
+      y[1:,i]=(u0*y0[i]+fI[:,i]+fU[:,i])/fTotal[1:]
+      y[0,i]=y0[i]
+
+   iT=model.lut['total']
+   y[:,iT]=fTotal
+ 
+#(N,n)
+   
+   return y
+