RuleOfCramer/main.py

'''
Created on 06.08.2016

@author: wn
'''

import time
import random


startTime = 0
lastTime = 0

def log(msg, reset=False):
    global startTime, lastTime
    currentTime = time.time()
    if reset:
        startTime = currentTime
    if lastTime == 0:
        lastTime = currentTime
    print("%f %f: %s" % (currentTime - startTime, currentTime - lastTime, msg))
    lastTime = currentTime

def s_permutations(seq):
    items = [[]]
    for j in seq:
        new_items = []
        for i, item in enumerate(items):
            if i % 2:
                # step up
                new_items += [item[:i] + [j] + item[i:]
                              for i in range(len(item) + 1)]
            else:
                # step down
                new_items += [item[:i] + [j] + item[i:]
                              for i in range(len(item), -1, -1)]
        items = new_items
 
    return [(tuple(item), -1 if i % 2 else 1)
            for i, item in enumerate(items)]
    
cachedIndexJ = {}

def det(m):
    log("Calculating determinant started")

    l = len(m)
    indexListI = range(1, l+1)
    if l in cachedIndexJ:
        indexListJ = cachedIndexJ[l]
        log("J for %i taken from cache" % l)
    else:
        indexListJ = s_permutations(indexListI)
        cachedIndexJ[l] = indexListJ
        cnt = len(indexListJ)
        #print(indexListJ)
        log("J for %i generated, %i entries" % (l, cnt))

    summ = 0
    steps = 0
    for indexJ in indexListJ:
        sign = indexJ[1]
        prod = sign
        for indexI in indexListI:
            prod *= m[indexI-1][indexJ[0][indexI-1]-1]
            steps += 1
        summ += prod
    log("Calculating determinant done, %i steps" % steps)

    return summ



def solve(g):
    coeffMatrix = g[:-1]
    coeffDet = det(coeffMatrix)

    order = len(g)-1
    unknowns = []
    for i in range(order):
        varMatrix = []
        for j in range(order):
            if (i == j):
                varMatrix.append(g[order])
            else:
                varMatrix.append(g[j])
        varDet = det(varMatrix)
        var = varDet / coeffDet
        unknowns.append(var)
    
    return unknowns


def generate(order, lowerBound, upperBound):
    log("Generating system started")
    matrix = []
    for i in range(order+1):
        vector = []
        for j in range(order):
            vector.append(random.uniform(lowerBound, upperBound))
        matrix.append(vector)
    log("Generating system done")
    return matrix


def verify(m, v):
    coeffs = m[:-1]
    results = m[-1]
    epsilons = []
    good = True
    MAX_EPSILON = 1E-10
    worstEpsilon = 0.0
    
    for i in range(len(results)):
        res = 0.0
        for j in range(len(results)):
            res += coeffs[j][i] * v[j]
        epsilon = res - results[i]
        if epsilon > MAX_EPSILON:
            good = False
        if abs(epsilon) > worstEpsilon:
            worstEpsilon = abs(epsilon)
        epsilons.append(res - results[i])
        
    return (good, worstEpsilon, epsilons)
    


#g = [[1.0,2.0,10.0], [5.0,3.0,11.0], [7.0,4.0,3.0], [8.0,9.0,13.0]]
#g = [[599.1290806728277, 80.40433203454268, 615.5120772027384], [187.33506822036273, 663.7490700388757, 677.5997857592452], [956.929845461785, 456.8619830134234, 419.21973493934263], [61.824473951393344, 342.8555419523176, 216.32654519677496]]
#g = [[8.019558599768107, 903.9449140455304, 919.612242767639, 753.9105365307668, 59.01030128588536, 920.4023338915356, 202.58884034378954, 431.09289053573997, 339.17491130065537, 313.11544472585905], [317.26651051567023, 727.3387587224795, 152.77343043145507, 408.44599685908923, 155.79032672640025, 425.6401391270683, 864.8657207590146, 298.9063131849389, 430.3786466384899, 425.29127250155807], [865.8564877305874, 669.481784333669, 538.9130960586408, 961.8480073483083, 753.9136870076827, 912.3571914611318, 929.0546428258709, 873.5276061317886, 459.59802608296684, 455.8093236762648], [318.8613025219013, 497.149024070492, 343.8577826382859, 780.6069003672162, 319.67906522418554, 628.3446507753014, 37.526503978800974, 490.48198671795006, 527.0266127960343, 180.8362568347207], [156.8237212311948, 510.1641000612821, 229.2078194439, 637.1024703370609, 972.0604419040393, 571.3909858582223, 184.34543232118418, 479.07607501617986, 382.12667148701297, 150.05251231416116], [297.8478733948842, 981.2792974599201, 271.8121041641507, 990.1956266273802, 645.839991521914, 288.9135934264093, 845.7496599336642, 19.88633912497506, 97.46153065555741, 343.1102890796861], [384.5208626880663, 145.38972551894247, 560.1624430301523, 224.07537517241872, 704.7921876231894, 781.2588369275445, 752.0786768851573, 77.91470589306748, 851.5354388193074, 146.64240712820654], [252.4899031289375, 928.4755757806711, 176.3416092863044, 914.652591937962, 379.5647442648579, 418.6432906860023, 397.44937187401274, 943.8352970973906, 495.8158468453109, 746.6463673774376], [330.80903398935936, 560.2262982840332, 881.0811315973745, 181.77380227358952, 315.50808931102347, 158.5072721277919, 947.7031957044982, 929.3796342482069, 131.55810107335486, 982.8660461523164], [310.9354667285096, 770.115136026887, 334.49254019770905, 236.96386990210726, 918.8431629974137, 893.3379619077041, 431.3668259114833, 285.69624238120986, 64.53926288027478, 799.2392946390154], [20.70259163693389, 306.46447894527165, 420.64284211501723, 37.622373302959474, 336.4417266492545, 247.6600784425259, 215.21800533449152, 426.10683135701044, 783.5055424440518, 817.507304763315]]


for i in range(2, 10+1):
    log("Starting ...", reset=True)
    g = generate(i, 0.0, 1000.0)

    v = solve(g)
    log("Solved")
    print(i)
    print(g)
    print(v)
    
    e = verify(g, v)
    print(e)
    print