from numpy import *
import bisect

"""
# Create an ixj matrix filled with 0s
>>> A = [3,4,8,4]
>>> [[[I(i,j) for j in range(len(A))]] for i in range(len(A))]
[[[0, 0, 0, 0]], [[0, 0, 0, 0]], [[0, 0, 0, 0]], [[0, 0, 0, 0,]]]
"""
def I(i,j) : return 0
# Create an ixj matrix that "expands" list A into 1s
# For example, if A is [5,6,2,3,4],
# Row 1 will have 5 1's in positions 0,1,2,3,4
# Row 2 will have 6 1's in positions 5,6,7,8,9,10
# Row 3 will have 2 1's in positions 11,12
# Row 4 will have 3 1's in positions 13,14,15
# Row 5 will have 4 1's in positions 16,17,18,19
"""
>>> A = [3,4,8,4]
>>> [[[II(i,j,[sum(A[0:x]) for x in range(len(A)+1)]) for j in range(sum(A))]] for i in range(len(A))]
[[[1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1]]]
"""
def II(i,j,A) :
  return 1 if bisect.bisect(A,j) - 1 == i else 0

"""
# Create an ixj matrix similar to III, except rotated 90 degrees & the expansion
# occurs only at one point, which is the midpoint of a given region.
# If the region has an even number of elements, the midpoint is split in 2.
# For example, if A is [5,6,2,3,4],
# Column 1 will have 1 at position 2
# Column 2 will have 0.5 at position 7 and 0.5 at position 8
# Column 3 will have 0.5 at position 11 and 0.5 at position 12
# Column 4 will have 1 at position 14
# Column 5 will have 0.5 at position 17 and 0.5 at position 18
>>> A = [3,4,8,4]
>>> [[[III(i,j,sum(A),len(A),[sum(A[0:x]) for x in range(len(A)+1)]) for j in range(sum(A))]] for i in range(len(A))]
[[[0, 0, 0, 0]], [[1, 0, 0, 0]], [[0, 0, 0, 0]], [[0, 0, 0, 0]], [[0, 0.5, 0, 0]], [[0, 0.5, 0, 0]], [[0, 0, 0, 0]], [[0, 0, 0, 0]], [[0, 0, 0, 0]], [[0, 0, 0, 0]], [[0, 0, 0.5, 0]], [[0, 0, 0.5, 0]], [[0, 0, 0, 0]], [[0, 0, 0, 0]], [[0, 0, 0, 0]], [[0, 0, 0, 0]], [[0, 0, 0, 0.5]], [[0, 0, 0, 0.5]], [[0, 0, 0, 0]]]
"""
def III(i,j,cardA,cardB,integralA) :
  if (integralA[j+1]- integralA[j]) % 2 :
    if (i == int((integralA[j+1]- integralA[j])/2. + integralA[j])) : return 1
    return 0
  if (i == int((integralA[j+1] - integralA[j])/2. + integralA[j])) : return 1/2.
  if (i == int((integralA[j+1] - integralA[j])/2. + integralA[j]) - 1) : return 1/2.
  return 0

"""
# Creates 1x1 rows of 0,0,0,...,0,0.5,-1,0.5,0,...0,0,0
# Where the position of 0.5,-1,0.5 shifts by 1 to the right and wraps around the end if necessary
>>> A = [3,4,8,4]
>>> [[[IV(i,j,19,4) for j in range(sum(A))]] for i in range(sum(A))]
[[[-1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5]], [[0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5]], [[0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1]]]
"""
def IV(i,j,cardA,cardB) :
  if j == i : return -1
  if (i < 1) & (j > cardA + i - 2) :   
    j = j - cardA
  if (i >= cardA - 1) & (j <= i - cardA + 1) : 
    j = j + cardA
  if abs(i - j) > 1 : return 0
  return 1/2.

"""
# Creates matrix with properties:
# [ I(len(A)xlen(A))   II(len(A)xsum(A))]
# [ III(sum(A)xlen(A)) IV(sum(A)xsum(A))]
>>> A = [3,4,8,4]
>>> [[[M(i,j,sum(A),len(A),[sum(A[0:x]) for x in range(len(A)+1)]) for j in range(sum(A)+len(A))]] for i in range(sum(A)+len(A))]
[[[0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1]], [[0, 0, 0, 0, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5]], [[1, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0.5, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0.5, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0, 0]], [[0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0, 0]], [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0, 0]], [[0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5, 0]], [[0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1, 0.5]], [[0, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, -1]]]
"""
def M(i,j,cardA,cardB,A) :
  if (i < cardB) & (j < cardB) : return I(i,j)
  if (i < cardB) : return II(i,j-cardB,A)
  if (j < cardB) : return III(i-cardB,j,cardA,cardB,A)
  return IV(i-cardB,j-cardB,cardA,cardB) 

"""
# Creates a vector that fills the first len(B) values with the values of B and 0 for everything else
>>> B = [1.3,3.3,5.5,4.4]
>>> [V(i,B,sum(A)) for i in range(sum(A)+len(A))]
[1.3, 3.2999999999999998, 5.5, 4.4000000000000004, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
"""
def V(i,B,cardA) :
  if i < len(B) : return B[i]
  return 0

"""
# Solves the equation M*X = V
>>> A = [3,4,8,4]
>>> B = B = [1.3,3.3,5.5,4.4]
>>> solve_sys(A,B)
array([-0.25131579,  0.14270243, -0.11447368,  0.22308704,  0.56787449,
        0.26578947,  0.46633603,  0.66688259,  0.86742915,  0.92527328,
        0.84041498,  0.75555668,  0.67069838,  0.58584008,  0.50098178,
        0.53059717,  0.67468623,  0.8187753 ,  0.96286437,  1.10695344,
        1.25104251,  1.17204453,  0.86995951])
"""
def solve_sys(A, B) :
  integralA = [sum(A[0:x]) for x in range(len(A)+1)]
  cardA = integralA[-1]
  cardB = len(A)
  m = array([[M(i,j,cardA,cardB,integralA) for j in range(cardA+cardB)] for i in range(cardA+cardB)])
  v = array([V(x,B,cardA) for x in range(cardA+cardB)])
  # REPLACE ME!!!
  result = linalg.solve(m,v)
  return result


# Some tests
BIG = [[[10,10,20,10,20,20,20,30],
        [1.0,1.5,2.5,3.5,3.3,3.3,1.4,3.3]],
       [[20,10,20,30,20,20,20,25],
        [1.0,1.5,2.5,3.5,3.3,3.3,1.4,3.3]],
       [[20,10,20,30,20,20,20],
        [1.0,1.5,2.5,3.5,3.3,3.3,1.4]],
       [[10,50,20,45],
        [1.0,4.0,2.5,3.5]],
       [[20,20,20,20,25],
        [1.0,1.5,1.5,1.0,1.3]],
       [[20,20,40,20,10,15],
        [2.0, 1.8,3.3,1.8,1.6,1.5]],
       [[20,20,20,20,20,20],
        [1.0,1.0,4.0,1.0,1.0,0.4]],
       [[20,20,40,20,10,15],
        [1.0,2.8,2.1,4.0,1.6,1.0]]]

for A,B in BIG :
  a = solve_sys(A,B)
  print "RESULT"
  print a
  integralA = [sum(A[0:x]) for x in range(len(A)+1)]
  print "SHOULD BE", B
  print "IS",
  for x in range(len(A)) :
    print sum(a[integralA[x]+len(B):integralA[x+1]+len(B)]),
  print " "
