#!/usr/bin/env python

import numpy

def solve(point_list):
    """
    This function solves the linear equation system involved in the n
    dimensional linear extrapolation of a vector field to an arbitrary point.

    f(x) = x * A + b

    with:
        A - The "slope" of the affine function in an n x n matrix.
        b - The "offset" value for the n dimensional zero vector.

    The function takes a list of n+1 point-value tuples (x, f(x)) and returns
    the matrix A and the vector b. In case anything goes wrong, the function
    returns the tuple (None, None).

    These can then be used to compute directly any value in the linear
    vector field.
    """
    # Some helpers.
    dimensions = len(point_list[0][0])
    unknowns = dimensions ** 2 + dimensions
    number_points = len(point_list[0])

    # Bail out if we do not have enough data.
    if number_points < unknowns:
        print 'For a %d dimensional problem I need at least %d data points.' \
              % (dimensions, unknowns)
        print 'Only %d data points were given.' % number_points
        return None, None

    # Ensure we are working with a NumPy array.
    point_list = numpy.asarray(point_list)

    # For the solver we are stating the problem as
    # C * x = d
    # with the problem_matrix C and the problem_vector d

    # We're going to feed our linear problem into these arrays.
    # This one is the matrix C.
    problem_matrix = numpy.zeros([unknowns, unknowns])
    # This one is the vector d.
    problem_vector = numpy.zeros([unknowns])

    # Populate data matrix C and vector d.
    x_values, y_values = point_list[0], point_list[1]
    for i in range(dimensions):
        x_i, y_i = x_values[:, i], y_values[:, i]
        for j in range(dimensions):
            y_j = y_values[:, j]
            row = dimensions * i + j
            problem_vector[row] = (x_i * y_j).sum()
            problem_matrix[row, dimensions ** 2 + j] = x_i.sum()
            problem_matrix[dimensions ** 2 + j, dimensions * i + j] = x_i.sum()
            for k in range(dimensions):
                x_k = x_values[:, k]
                problem_matrix[row, dimensions * k + j] = (x_i * x_k).sum()
        row = dimensions ** 2 + i
        problem_vector[row] = y_i.sum()
        problem_matrix[row, dimensions ** 2 + i] = number_points


    matrix_A, vector_b = None, None
    try:
        result_vector = numpy.linalg.solve(problem_matrix, problem_vector)

        # Check whether we really did get the right answer.
        # This is advised by the NumPy doc string.
        if numpy.linalg.norm(numpy.dot(problem_matrix, result_vector)
                             - problem_vector) < 1e-6:
            # We're good, so hack up the result into the matrix and vector.
            matrix_A = result_vector[:dimensions ** 2]
            matrix_A.shape = (dimensions, dimensions)
            vector_b = result_vector[dimensions ** 2:]
        else:
            print "For whatever reason our linear equations didn't solve."
            print numpy.linalg.norm(numpy.dot(result_vector, problem_matrix)
                                    - problem_vector)
    except numpy.linalg.linalg.LinAlgError:
        print "Things didn't work out as expected, eh."

    return matrix_A, vector_b


def main():
    """
    We're testing the 3D case of a mapping

        y = x * A + b

    We are setting up a mystery matrix (A) and vector (b), creating a solid
    number of random samples with noise from that, and then we're finding
    out a good approximation for A and b by solving the regression problem.
    The approximation is compared to the original mystery values.
    """
    # We're testing the 3D case. And using 100 data points.
    dimensions = 3
    data_points = 100

    # Let's make a mystery matrix and a mystery vector with elements < 10.0.
    mystery_matrix = numpy.random.uniform(low=0.0, high=10.0,
                                          size=dimensions ** 2)
    mystery_matrix.shape = (dimensions, dimensions)
    mystery_vector = numpy.random.uniform(low=0.0, high=10.0,
                                          size=dimensions)

    # Now let's make 100 sample points ...
    x_values = numpy.random.uniform(low=0.0, high=10.0,
                                    size=dimensions * data_points)
    x_values.shape = (data_points, dimensions)
    
    # ... and calculate their corresponding values ...
    y_values = numpy.dot(x_values, mystery_matrix) + mystery_vector

    # ... and add a bit of noise.
    noise_scale = 1.5
    noise = numpy.random.normal(loc=0.0, scale=noise_scale,
                                size=dimensions * data_points)
    noise.shape = (data_points, dimensions)
    y_values += noise

    # Solve the n-D linear regression problem.
    estimated_A, estimated_b = solve([x_values, y_values])

    if estimated_A is not None:
        # And test it on a nice point somewhere
        dest_point = numpy.random.uniform(low=0.0, high=10.0,
                                          size=dimensions)
        expected_value = numpy.dot(dest_point, mystery_matrix) + mystery_vector
        estimated_value = numpy.dot(dest_point, estimated_A) + estimated_b
        distance = numpy.linalg.norm(expected_value - estimated_value)
        
        print 'Point at:\n\t%s' % dest_point.round(3)
        print 'Expected:\n\t%s' % expected_value.round(3)
        print 'Estimated:\n\t%s' % estimated_value.round(3)
        print 'Distance:\n\t%s' % distance.round(3)
    else:
        print 'Sorry, problem could not be solved.'


if __name__ == '__main__':
    main()
<br />
<b>Fatal error</b>:  Call to a member function getRouter() on a non-object in <b>/var/www/ppapers/ojs.pythonpapers.org/htdocs/lib/pkp/classes/template/PKPTemplateManager.inc.php</b> on line <b>64</b><br />