{-# LANGUAGE ConstraintKinds, ScopedTypeVariables, TypeApplications, AllowAmbiguousTypes, NoImplicitPrelude,
    GeneralizedNewtypeDeriving, FlexibleContexts
#-}
module LinRel where

import Numeric.LinearAlgebra
import Prelude hiding ((<>))
import Debug.Trace
type BEnum a = (Enum a, Bounded a) 
enumAll :: (BEnum a) => [a] -- What about Void?
enumAll = [minBound .. maxBound]



-- A Hacked "Void" type to get card @Void = 0
data Void = Void' | Void'' -- SORRY MOM.
instance Enum Void where
  fromEnum Void' =  1
  fromEnum Void'' = 0
  toEnum 1 = Void'
  toEnum 0 = Void''
instance Bounded Void where
  maxBound = Void''
  minBound = Void'

instance (Enum a, Enum b, Bounded a) => Enum (Either a b) where
  fromEnum (Left a) = fromEnum a
  fromEnum (Right b) = fromEnum b + fromEnum (maxBound @a) + 1
  toEnum n | n <= ((fromEnum (maxBound @a)) - (fromEnum (minBound @a))) = Left (toEnum n)
           | otherwise = Right (toEnum (n - fromEnum (maxBound @a))) 


instance (Bounded a, Bounded b) => Bounded (Either a b) where
  maxBound = Right maxBound
  minBound = Left minBound

card :: forall a. (BEnum a) => Int
card = (fromEnum (maxBound @a)) - (fromEnum (minBound @a)) + 1

-- HLinRel holds A x = b constraint
data HLinRel a b = HLinRel (Matrix Double) (Vector Double) deriving Show

-- x = A l + b. Generator constraint. 
data VLinRel a b = VLinRel (Matrix Double) (Vector Double) deriving Show


-- f(x) = xQx + bx
data QuadOp a b = QuadOp (Matrix Double) (Vector Double)
{-
qid :: QuadOp a a
qid = ? -- a = a ? Need Relations too.

qcompose :: QuadOp b c -> QuadOp a b -> QuadOp a c
qcompose (QuadOp q c) (QuadOp q' c') = QuadOp q'' c'' where
        where 
                ca = card @a
                cb = card @b 
                cc = card @c
                a = subMatrix (cb, cb) q
                b = subMatrix (cb,cc) q
                c = subMatrix (cc, cb) q
                d = subMatrix (cc,cc) q
                a' = subMatrix (ca, ca) q'
                b' = subMatrix (ca, cb) q'
                c' = subMatrix (cb, ca  q'
                d' = subMatrix (cb, cb) q'
                m = - (a + d')
                q'' = fromBlocks [[a' + b' <> m </> c'   , b' <> m </> c ],  -- can memoize some of this
                                  [b <> m </> c' , d - ]]
                
                a'' =  a' - b' <> m </> c'
                [v1', v2'] = takesV [ca,cb] c'
                [v1,  v2] = takesV [cb, cc] c
                v3 = (v2' + v1) 
                c'' = vJoin [v1' + m </> c #>  v3, ]

        -}
    {-
    
    Any cost outside the constraint space is irrelevant.
    Q = V m V where m is fullrank.
    c = Vc also

    minimal sets
    -}    
        
        -- break into pieces, form schur complement.

-- if A x = b then x is in the nullspace + a vector b' solves the equation
h2v :: HLinRel a b -> VLinRel a b
h2v (HLinRel a b) = VLinRel a' b' where
        b' = a <\> b -- least squares solution
        a' = nullspace a

-- if x = A l + b, then A' . x = A' A l + A' b = A' b because A' A = 0
v2h :: VLinRel a b -> HLinRel a b
v2h (VLinRel a' b') = HLinRel a b where
        b = a #> b' -- matrix multiply
        a = tr $ nullspace (tr a') -- orthogonal space to range of a.

hid :: forall a. BEnum a => HLinRel a a
hid =  HLinRel (i ||| (- i)) (vzero s) where 
                            s = card @a
                            i = ident s

vzero :: Konst Double d c => d -> c Double
vzero = konst 0

hcompose :: forall a b c. (BEnum a, BEnum b, BEnum c) => HLinRel b c -> HLinRel a b -> HLinRel a c
hcompose (HLinRel m b) (HLinRel m' b') = let a'' = fromBlocks [[       ma',           mb' ,    0       ],
                                                               [         0 ,    mb,        mc          ]] in
                                         let b'' = vjoin [b', b] in 
                                         let (VLinRel q p) = h2v (HLinRel a'' b'') in -- kind of a misuse
                                         let q' = (takeRows ca q)  -- drop rows belonging to @b
                                                       === 
                                                  (dropRows (ca + cb) q) in
                                         let [x,y,z] =  takesV [ca,cb,cc] p in
                                         let p'=  vjoin [x,z] in -- rebuild without rows for @b
                                         v2h (VLinRel q' p') -- reconstruct HLinRel
                                       where 
                                           ca = card @a
                                           cb = card @b 
                                           cc = card @c
                                           sb = size b -- number of constraints in first relation
                                           sb' = size b' -- number of constraints in second relation
                                           ma' = takeColumns ca m'
                                           mb' = dropColumns ca m'
                                           mb = takeColumns cb m
                                           mc = dropColumns cb m

(<<<) :: forall a b c. (BEnum a, BEnum b, BEnum c) => HLinRel b c -> HLinRel a b -> HLinRel a c
(<<<) = hcompose
-- stack the constraints
hmeet :: HLinRel a b -> HLinRel a b -> HLinRel a b
hmeet (HLinRel a b) (HLinRel a' b') = HLinRel (a === a') (vjoin [b,b'])




{- If they don't meet are we still ok? 

I am not sure. Might be weird corner cases?

-}

hjoin :: HLinRel a b -> HLinRel a b -> HLinRel a b
hjoin v w = v2h $ vjoin' (h2v v) (h2v w)

-- hmatrix took vjoin from me :(
-- joining means combining generators and adding a new generator
-- Closed under affine combination l * x1 + (1 - l) * x2 
vjoin' :: VLinRel a b -> VLinRel a b -> VLinRel a b
vjoin' (VLinRel a b) (VLinRel a' b') = VLinRel (a ||| a' ||| (asColumn (b - b'))) b

-- no constraints, everything
-- trivially true
htop :: forall a b. (BEnum a, BEnum b) => HLinRel a b 
htop = HLinRel (vzero (1,ca + cb)) (konst 0 1) where 
                                      ca = card @a
                                      cb = card @b 
{-
hbottom :: forall a b. (BEnum a, BEnum b) => HLinRel a b 
hbottom = HLinRel (vzero (1,ca + cb)) (konst 1 1) where 
                                      ca = card @a
                                      cb = card @b       
                                      -}                                
-- all the constraints! Only the origin.
-- no. it should be the empty set. Impossible to satisfy.
-- 0 x = 1 is impossible
-- not gonna play nice.
{-
hbottom :: forall a b. (BEnum a, BEnum b) => HLinRel a b 
hbottom = HLinRel (ident (ca + cb)) (konst 0 (ca + cb)) where 
                                    ca = card @a
                                    cb = card @b  
  -}                            

hconverse :: forall a b. (BEnum a, BEnum b) => HLinRel a b -> HLinRel b a 
hconverse (HLinRel a b) = HLinRel ( (dropColumns ca a) |||  (takeColumns ca a)) b where 
    ca = card @a
    cb = card @b  

    -- this is numerically unacceptable
-- forall l. A' ( A l + b) == b'
vhsub :: VLinRel a b -> HLinRel a b -> Bool
vhsub (VLinRel a b) (HLinRel a' b') = (naa' <=  1e-10 * (norm_2 a') * (norm_2 a)  ) && ((norm_2 ((a' #> b) - b')) <= 1e-10 * (norm_2 b')  ) where
          naa' = norm_2 (a' <> a)

hsub :: HLinRel a b -> HLinRel a b -> Bool
hsub h1 h2 = vhsub (h2v h1) h2

heq :: HLinRel a b -> HLinRel a b -> Bool
heq a b = (hsub a b) && (hsub b a)


instance Ord (HLinRel a b) where
  (<=) = hsub
  (>=) = flip hsub 

instance Eq (HLinRel a b) where
  (==) = heq


-- I can't do this right?
-- hcomplement :: HLinRel a b -> HLinRel a b
-- hcomplement  

hpar :: HLinRel a b -> HLinRel c d -> HLinRel (Either a c) (Either b d)
hpar (HLinRel mab v) (HLinRel mcd v') = HLinRel (fromBlocks [ [mab, 0], [0 , mcd]]) (vjoin [v, v']) where


hassoc :: forall a b c. (BEnum a, BEnum b, BEnum c) => HLinRel (Either (Either a b) c) (Either a (Either b c))
hassoc = HLinRel m v where HLinRel m v = hid @((Either (Either a b) c))

hassoc' :: forall a b c. (BEnum a, BEnum b, BEnum c) => HLinRel  (Either a (Either b c)) (Either (Either a b) c)
hassoc' = HLinRel m v where HLinRel m v = hid @((Either (Either a b) c))

{-
        -- Void is unit for Either.
-- void has no inhabitants.... This is a bad boy.


hcup :: HLinRel Void (Either a a)
hcap :: HLinRel (Either a a) Void




-}
hleft :: forall a b. (BEnum a, BEnum b) => HLinRel a (Either a b)
hleft = HLinRel ( i ||| (- i) ||| (konst 0 (ca,cb))) (konst 0 ca) where 
    ca = card @a
    cb = card @b  
    i = ident ca

hright :: forall a b. (BEnum a, BEnum b) => HLinRel b (Either a b)
hright = HLinRel ( i ||| (konst 0 (cb,ca)) ||| (- i) ) (konst 0 cb) where 
    ca = card @a
    cb = card @b  
    i = ident cb


hfan ::  forall a b c. BEnum a => HLinRel a b -> HLinRel a c -> HLinRel a (Either b c)
hfan (HLinRel m v) (HLinRel m' v') = HLinRel (fromBlocks [ [ma, mb, 0], [ma', 0, mc']]) (vjoin [v,v']) where
        ca = card @a
        ma = takeColumns ca m 
        mb = dropColumns ca m 
        ma' = takeColumns ca m' 
        mc' = dropColumns ca m' 


hdump :: HLinRel a Void
hdump = HLinRel 0 0
{-
hlabsorb :: HLinRel a b -> HLinRel (Either Void a) b
hlabsorb (HLinRel m v) = (HLinRel m v)
-}
hlabsorb ::forall a. BEnum a => HLinRel (Either Void a) a
hlabsorb = HLinRel m v where (HLinRel m v) = hid @a 

htrans :: HLinRel a (Either b c) -> HLinRel (Either a b) c 
htrans (HLinRel m v) = HLinRel m v

hswap :: forall a b. (BEnum a, BEnum b) => HLinRel (Either a b) (Either b a)
hswap = HLinRel (fromBlocks [[ia ,0,0 ,-ia], [0, ib,-ib,0]]) (konst 0 (ca + cb)) where 
        ca = card @a
        cb = card @b  
        ia = ident ca
        ib = ident cb


hsum :: forall a. BEnum a => HLinRel (Either a a) a
hsum = HLinRel ( i ||| i ||| - i ) (konst 0 ca)  where 
        ca = card @a 
        i= ident ca

hdup :: forall a. BEnum a => HLinRel a (Either a a)
hdup = HLinRel (fromBlocks [[i, -i,0 ], [i, 0, -i]]) (konst 0 (ca + ca))  where 
        ca = card @a 
        i= ident ca

-- hcup :: forall a. BEnum a => HLinRel Void (Either a a)
-- hcup = -- or mainulate hid

-- hcap :: forall a. BEnum a => HLinRel (Either a a) Void

-- smart constructors
hLinRel :: forall a b. (BEnum a, BEnum b) => Matrix Double -> Vector Double -> Maybe (HLinRel a b) 
hLinRel m v | cols m == (ca + cb) && (size v == rows m)  = Just (HLinRel m v)
            |  otherwise = Nothing  where 
                 ca = card @a
                 cb = card @b  

-- a 2d space at every wire or current and voltage.
data IV = I | V deriving (Show, Enum, Bounded, Eq, Ord)


resistor :: Double -> HLinRel IV IV
resistor r = HLinRel ( (2><4)  [ 1,0,-1,   0,
                                 r, 1, 0, -1]) (konst 0 2)  

bridge :: Double -> HLinRel (Either IV IV) (Either IV IV)
bridge r = HLinRel (  (4><8) [ 1,0, 1,  0, -1, 0, -1,  0, -- current conservation
                               0, 1, 0, 0, 0, -1 , 0,  0, --voltage maintained
                               0, 0, 0, 1, 0,  0,  0, -1, -- voltage maintained
                               r, 1, 0,-1, -r,  0,  0, 0  ]) (konst 0 4)  
short = bridge 0



first :: BEnum c => HLinRel a b -> HLinRel (Either a c) (Either b c)
first f = hpar f hid 

second :: BEnum a => HLinRel b c -> HLinRel (Either a b) (Either a c)
second f = hpar hid f

type HLinRel2D u d l r = HLinRel (Either u l) (Either d r)


{-
A stencil  of 2d resistors for tiling


          u
          /
          \
          /
         |
l -/\/\/----/\/\/\-  r
         |
         /
         \
         /
         |
         d


        -}
stencil :: HLinRel2D IV IV IV IV
stencil = (hpar r10 r10) <<< short <<< (hpar r10 r10) where r10 = resistor 10

horicomp :: forall w w' w'' w''' a b c. (BEnum w, BEnum w', BEnum a, BEnum w''', BEnum b, BEnum w'', BEnum c ) => HLinRel2D w' w'' b c -> HLinRel2D w w''' a b -> HLinRel2D (Either w' w) (Either w'' w''') a c
horicomp f g = hcompose f' g' where 
               f' :: HLinRel (Either (Either w' w''') b) (Either (Either w'' w''') c)
               f' = (first hswap) <<< hassoc' <<< (hpar hid f) <<< hassoc <<<  (first hswap) 
               g' :: HLinRel (Either (Either w' w) a) (Either (Either w' w''') b)
               g' = hassoc' <<< (hpar hid g) <<< hassoc


rotate :: (BEnum w, BEnum w', BEnum a, BEnum b) => HLinRel2D w w' a b -> HLinRel2D a b w w'                                      
rotate f = hswap <<< f <<< hswap

vertcomp :: (BEnum w, BEnum w', BEnum a, BEnum d, BEnum b, BEnum w'', BEnum c ) => HLinRel2D w'  w'' c d -> HLinRel2D w w' a b -> HLinRel2D w w'' (Either c a) (Either d b)
vertcomp f g = rotate (horicomp (rotate f)  (rotate g) ) 

{-
traceV :: ->
traceH :: ->


  -}

stencil2 = vertcomp h h where h = stencil `horicomp` stencil



{- Legendre transformations in thermo are for open systems. SOmething to that -}
{- Dependent sources. Well, these are goddamn cheating.
We could do it though
|   |
I   alpha I
|   |

small signal models of transistor, op amps etc


gyrator - would need polynomials from caps and inds
would be kind of nice for boundary models

kron instead of dsum - Quantum fields
Kron relation might be nice for discussing remnant space of topological matter



-}

{-

A wire in a circuit has a potential (questionably) and a current running through it
So our wires should have both of these variables.

-}
newtype VProbe = VProbe () deriving (Enum, Bounded, Show, Eq, Ord)
vprobe :: HLinRel IV VProbe
vprobe = HLinRel ( (2><3)  [1,0,0,
                            0,1,-1]) (konst 0 2)                

vsource :: Double -> HLinRel IV IV
vsource v = HLinRel ( (2><4) [ 1,0,-1,   0,
                               0, 1, 0, -1]) (fromList [0,v])  

isource :: Double -> HLinRel IV IV
isource i = HLinRel ( (2><4) [  1,0, -1,   0 , -- current conservation
                                 1, 0, 0,  0]) (fromList [0,i])  


-- the currents add, but the voltages dup. sum and dup are dual
-- Or should it be |--|   a parallel short?
-- Ad then we could open circuit one of them and absorb the Void
-- to derive this
cmerge :: HLinRel (Either IV IV) IV
cmerge = HLinRel ( (3><4)  [1, 0, 1, 0, -1, 0,
                            0,1,0,0,0 ,  -1  ,
                            0,0,0,1, 0, -1])  (konst 0 3)

open :: HLinRel IV Void
open = HLinRel ( (1><2) [1,0]) (konst 0 1)


cap :: HLinRel  (Either IV IV) Void
cap  = hcompose open cmerge

cup :: HLinRel Void (Either IV IV)
cup = hconverse cap

ground :: HLinRel IV Void
ground = HLinRel ( (1><2) [ 0 , 1 ]) (vzero 1) 

-- resistors in parallel.

ex1 = hcompose (bridge 10) (bridge 10)
ex2 = hcompose (resistor 10) (resistor 30) -- resistors in series.
r20 :: HLinRel IV IV
r20 = resistor 20

divider :: Double -> Double -> HLinRel (Either IV IV) (Either IV IV)
divider r1 r2 = hcompose (bridge r2) (hpar (resistor r1) hid) 

{-
type StateCoState s = Either s s
type ValueD s = s


dynamics :: HLinRel (Either SHOState Control) SHOState
dynamics = HLinRel ((2><5)  [ 1,  dt, 0, -1 ,  0,
                              -dt, 1, dt, 0, -1 ])  (vzero 2)
   where dt = 0.01  

-- labsorb <<< (par initial_cond id)
-}


-- state of an oscillator
data SHOState = X | P deriving (Show, Enum, Bounded, Eq, Ord)
data Control = F  deriving (Show, Enum, Bounded, Eq, Ord)
-- Costate newtype wrapper
newtype Co a = Co a deriving (Show, Enum, Bounded, Eq, Ord)

type M = Matrix Double
dynamics :: forall x u. (BEnum x, BEnum u) =>  Matrix Double -> Matrix Double ->  
  HLinRel (Either x u) x
dynamics a b =  HLinRel (a ||| b ||| -i )  (vzero cx) where
  cx = card @x
  cu = card @u
  i = ident cx

initial_cond :: forall x. BEnum x => Vector Double-> HLinRel Void x
initial_cond x0 =  HLinRel i x0 where
  cx = card @x
  i = ident cx

valueUpdate :: forall x l. (BEnum x, BEnum l) =>  M -> M -> HLinRel (Either x l) l
valueUpdate a q = HLinRel ((tr a) ||| q ||| i)  (vzero cl) where
  cl = card @l
  i = ident cl  

optimal_u :: forall u l. (BEnum u, BEnum l) => 
     M -> M -> HLinRel u l
optimal_u r b = HLinRel (r ||| tr b) (vzero cu) where
  cu = card @u
 
step :: forall x u l. (BEnum x, BEnum u, BEnum l) => M -> M -> M -> M 
        -> HLinRel (Either x l) (Either x l)
step a b r q = 
  f5 <<< f4 <<< hassoc' <<< f3 <<< f2 <<< hassoc <<< f1 where
  f1 :: HLinRel (Either x l) (Either (Either x x) l)
  f1 = first hdup
  f2 :: HLinRel (Either x (Either x l)) (Either x l)
  f2 = second (valueUpdate a q)
  f3 ::  HLinRel (Either x l) (Either x (Either l l))
  f3 = second hdup
  f4 ::  HLinRel (Either (Either x l) l) (Either (Either x u) l)
  f4 = first (second (hconverse (optimal_u r b)))
  f5 ::  HLinRel (Either (Either x u) l) (Either x l)
  f5 = first (dynamics a b)

-- iterate (hcompose (step a b r q)) :: [HLinRel (Either x l) (Either x l)]



 {-
step :: forall x l. (BEnum x, BEnum l) => 
   HLinRel (State + ValueD) (State + ValueD)
step a b q r =                
= dynamics, 
  valueUpdate <<< fst
  second optimal_u
  par (par hid cup) hid
   ::  , optimal_u
-}
{-
cost :: HLinRel (Either (Co SHOState) SHOState Control) (Co SHOState)
cost = HLinRel ((2><5) [ 1,  0, 0,  -1 ,  0,
                        -dt, 1,  dt,  0,  -1 ])  (vzero 2)
-}





{-

Is there a reasonable intepretation of kron?

-}
{-
everything can be definedc inefficiently via v2s and h2v functions

right division

-}

    {-
Call them affine relations

Join and meet aren't union and intersection.
They are the affine closure of union and intersection.




Linear has some niceness.
Homgeonous coordinates usually do.
For clarity and familiaryt I have chosebn not to do it this way
Or maybe I will do it?

par


-}


{-

import numpy as np


def meet(a,b):
    pass
def compose(a,b): # a after b
    assert(a.inN == b.outN)

    combo = np.block([[a.constraints, np.zeros((a.constraints.shape[0] , b.inN) )     ],
                      [np.zeros((b.constraints.shape[0] , a.outN)) ,      b.constraints]])
    print(combo)
    gens = LinRel(a.inN + a.outN, b.outN, gens=combo).gens
    print("gens",gens) 
    gens = np.vstack((gens[:a.outN, :], gens[-b.inN: , :]) )
    print(gens)
    return LinRel(a.outN, b.inN, gens=gens)
def top(outN, inN):
    return LinRel(outN, inN, constraints = np.array([[]]))
def bottom(outN, inN):
    return LinRel(outN, inN, gens = np.array([[]]))

def converse(a):
    return LinRel(a.inN, a.outN, np.hstack((a.constraints[:, a.inN:], a.constraints[:, :a.inN])))
def complement(a):
    return LinRel(a.outN, a.inN, constraints = a.gens.T.conj())

def right_div(a,b):
    pass # return complement( compose(a, complement(b)) ) # something like this
def inclusion(a,b):
    
    s = a.constraints @ b.gens
    np.all(a.constraints @ b.gens <= tol )

    if rcond is None:
        rcond = np.finfo(s.dtype).eps * max(max(a.shape), max(b.shape))
        tol = max(np.amax(a), np.amax(b))
        tol = np.amax(s) * rcond
def fromMat(mat):
    (outN, inN) = mat.shape
    return LinRel(inN, outN,constraints = np.hstack((mat,-np.eye(outN))))
def id(N):
    return fromMat(np.eye(N))
# make 0,1 first index for in/out? Oh, but then in out have to be same size.
# A[0, ...] @ x + A[1, ...] @ y = 0   
# Then I can form the kron of linear relations. 
# store sperate A B matrices? A @in + B @ out
class LinRel():
    def __init__(self, outN, inN, constraints = None, gens = None, rcond=None):
        #assert(inN <= constraints.shape[1])
        self.inN = inN
        self.outN = outN
        if constraints is not None: #baiscally scipy.linalg.null_space
            u, s, vh = np.linalg.svd(constraints, full_matrices=True)
            M, N = u.shape[0], vh.shape[1]
            if rcond is None:
                rcond = np.finfo(s.dtype).eps * max(M, N)
            tol = np.amax(s) * rcond
            num = np.sum(s > tol, dtype=int)
            self.gens = vh[num:,:].T.conj()
            self.constraints = vh[:num,:]
        if gens is not None: #basically scipy.linalg.orth
            u, s, vh = np.linalg.svd(gens, full_matrices=True)
            M, N = u.shape[0], vh.shape[1]
            if rcond is None:
                rcond = np.finfo(s.dtype).eps * max(M, N)
            tol = np.amax(s) * rcond
            num = np.sum(s > tol, dtype=int)
            self.gens = u[:, :num]
            self.constraints = u[:, num:].T.conj()

    def shape(self):
        return (self.outN, self.inN)
    def size(self):
        return self.outN + self.inN
    # operator overloadings
    def __matmul__(a,b):
        return compose(a,b)
    def __invert__(a): # ~
        return complement(a)
    def __or__(a,b): # an argument could be made for + and *
        return join(a,b)
    def __and__(a,b):
        return meet(a,b)
    def __sub__(a,b):
        return  (a) & (-b)
    def __le__(a,b): # Are the others automatic?
        return inclusion(a,b)
    def __str__(self):
        return " Constraints: \n%s, \nGens:\n %s\n" % (str(self.constraints), str(self.gens))




ex = LinRel(1,2, np.array([[3,4,0]]))
e2 = LinRel(2,1, np.array([[3,4,0]]))

assert(np.all(LinRel(1,2, ex.constraints).constraints == ex.constraints) )
print( ex.constraints @ ex.gens)
assert(np.all( np.abs(ex.constraints @ ex.gens) <= 1e-15) )
print(ex @ e2)
print(e2 @ ex)
print(e2 @ id(e2.inN))

'''
Quadratic optimization can be mixed in.
Quad && LinRel

AffineRel = maintain homgenous coord, or always insert 1 -1 keeping homoeg coord
+ discrete? Maintain a couple copies of 


'''


Linear relations
hrep - Ax = 0
or
vrep - y = sum x_i

hrep <-> vrep = row echelon

in and out variables. In and out subspaces.
in : [] - list of indices
out : []

in is a vrep of input space.

in = projection/injection matrix n x d
out = projection/injection matrix. (d-n) x d
in * out = 0. orthogonal

auxiliary variables allowed

compose relations
in1 out1
in2 out2

in = d x n
stack out1 and in2 into matrix. with them equal.

np.hstack(A1, out1 - in2, A2)
in = no.hatck(in, zeores)
in = no.hatck(zeores, out)

drect sum as monoidal product
block([ A, 0  ],
      [ 0, A  ])
in = [in1, in2]
out = [out1, out2]

converse = flip in out

meet = combine the two constraint matrices
join = convert? combine in out?

1-d as unit object
<= is subspace ordering

negation = orthogonalization
division => 


a linear problem is a linear relation of 1-d.

use in1 and out2 as new in/out

fan 
snd(30,10) = project bottom 10 = idetnity matrix stacked.
fst(30, 10) project top 10
id(20)
id(n) = LinRel(np.zeros(0,2*n), [zeros, eye], [eye, zero] ) 
id(n) = LinRel((0,n), eye(n), eye(n)  )
        LinRel [I, -I], [I, 0], [0,I]

"internal" space
class LinRel():
    def __init__(A, in, out):

svd(in * A , smallest )
1 - A*A

Ax + By = 0

vstack adds constraints
hstack adds variables

def idRel(n):
  return LinRel(sparse.lil_matrix((0,n)), n, n)
def mapRel(A):
  (r,c) = A.shape
  newA = sparse.hstack(A, - sparse.eye(r))
  return LinRel(newA, r, c) 

class LinRelV():
class LinRelH():

class LinRel():
  def init(self,A, in, out):
    self.A = A
    self.in = in
    self.out = out
  def compose(self, b):
    assert(self.out == b.in, "Shapes must match")
    
    i = sparse.eye(b.in)
    cons = sparse.hstack([0, i, -i, 0])
    ina = self.A[:,:self.in]
    auxa = self.A[:,self.in:-self.out]
    outa = self.A[:,-self.out:]
    inb = b.A[:,:b.in]
    auxb = b.A[:,b.in:-b.out]
    outb = b.A[:,-b.out:]

    newA = sparse.bmat(  [[ina, auxa, outa, 0, 0,    0],
                          [0  , 0   , i,   -i, 0 ,   0]
                          [0,  0,   0,    inb, auxb,outb]])
    LinRel(newA, self.in, b.out)
  def meet(self,b):
    #hmm. I suppose we acutlaly should split the thing apart again
    assert(self.in == b.in)
    assert(self.out == b.out)
    assert()
    newA = sparse.vstack([self.A, b.A])
    return ()
  def complement(self):
    linalg.svd(self.A)
    return LinRel(get_nonzeroeigs)
  def __negate__(self):
    return self.complement()
  def rdiv(self):
  def transpose():
    self.converse()
  def T():
    self.converse()
  
  # the svd gives you most of what you need?
  def inclusion():
    x = linalg.nullspace(self.A)
    return b.A @ x == 0 #check that every generator is in. Makes sense. Except numrically is trash.
  def __leq__(self):
    self.inclusion(b)

  

  def converse(rel):
      newA = np.hstack( [ rel.A[:,-rel.out: ] , rel.A[:,rel.in:-rel.out], rel.A[:,:rel.in ] ])
      return LinRel(newA, rel.out, rel.in)


compose(A)
linalg.nullspace(A)
range?



(cons1,d1) = self.A.shape
(cons2,d2) = b.A.shape
constrain = np.hstack 
newA = sparse.vstack 


using bigM I can encode the complement of a H-polytope 
but then what?
I do it again I guess?


complementation of relation ->
At least one must be inverted.

polytope inclusion
-> search for point not in B.
Or, do sandardinni encoding

Really It should generate new variables for an instantiation.
Snd should not reuse the same variables every time.
class PolyRel()
  invars = []
  constraints = [] # store >= values, not full constraints?
  outvars = []
  def __init__():
    all fresh variables
  def compose(self,b):
    PolyRel(self.constraints + self., invars = outvars )
  def complement():
    zs = []
    for c in constraints:
      z, constraints = reify(c)
      # z = cvx.Variable(1, boolean=True) # actually make same shape as c
      # c += c + M * z
    sum(zs) >= 1 # one or more constraints is disatisfied.
  def rdiv():

yeah. We should use a dsl, compile it, then encode it.
data Rel a b where
  Compose ::
  Complement ::
  Converse ::


relu

l1 >= 0
l2 >= 0
l1 <= M * z
l2 <= M * (1 - z)
x = lambda1 - lambda2
y = lambda1

Maybe insetado f subspaces, we should be thinking ellipses and svd.



-}