<commit_before># -*- coding: utf-8 -*-
#
"""Testing script for wlsqm, doubles as a usage example.

-JJ 2016-11-10
"""

from __future__ import division
from __future__ import absolute_import

import time

import numpy as np
import sympy as sy

import scipy.spatial  # cKDTree

import pylab as pl
import mpl_toolkits.mplot3d.axes3d as p3

try:
    import wlsqm
except ImportError:
    print "WLSQM not found; is it installed?"
    from sys import exit
    exit(1)

import sudoku_lhs


# from various scripts, e.g. miniprojects/misc/tworods/main2.py
def axis_marginize(ax, epsx, epsy):
    a = ax.axis()
    w = a[1] - a[0]
    h = a[3] - a[2]
    ax.axis( [ a[0] - w*epsx, a[1] + w*epsx, a[2] - h*epsy, a[3] + h*epsy] )


# from find_neighbors2.py
class SimpleTimer:
    def __init__(self, label="", n=None):
        self.label = label
        self.n     = n      # number of repetitions done inside the "with..." section (for averaging in timing info)

    def __enter__(self):
        self.t0 = time.time()
        return self

    def __exit__(self, errtype, errvalue, traceback):
        dt         = time.time() - self.t0
        identifier = ("%s" % self.label) if len(self.label) else "time taken: "
        avg        = (", avg. %gs per run" % (dt/self.n)) if self.n is not None else ""
        print "%s%gs%s" % (identifier, dt, avg)


# many simultaneous local models, 2D
#
def testmany2d():
    #########################
    # config
    #########################

    ntasks = 8  # OpenMP parallelization

    axislims = [0., 1., 0., 1.]  # [xmin, xmax, ymin, ymax], for plotting
    nvis = 201  # number of visualization points per axis

    expr    = sy.sympify("sin(pi*x) * cos(pi*y)")

    points_per_axis = 100  # for point cloud generation

    r      = 5e-2  # neighborhood radius
    max_nk = 100   # maximum number of neighbor points to accept into each neighborhood (affects memory allocation)

    knowns = 1  # function value is known

    fit_order = 4
    weighting_method = wlsqm.WEIGHT_CENTER

    max_iter = 10  # for iterative fitting method

    reps = 20  # for demonstration of solving multiple times using the same geometry

    #########################
    # the test itself
    #########################

    print
    print "=" * 79
    print "many neighborhoods, 2D case"
    print "=" * 79
    print

    # create a stratified point cloud
    print "generating sudoku sample"
    with SimpleTimer(label=("    done in ")) as s:
        S,m = sudoku_lhs.sample(2, points_per_axis, 1)
        bins_per_axis = m*points_per_axis
        S = S / float(bins_per_axis - 1)  # scale the sample from [0, bins_per_axis-1]**2 to [0, 1]**2
        npoints = len(S)
        print "    %d points" % (npoints)

    # index the point cloud for fast neighbor searching
    print "indexing sample"
    with SimpleTimer(label=("    done in ")) as s:
        tree = scipy.spatial.cKDTree( data=S )

    # If this was an IBVP, we would here get the previous state of the unknown field.
    #
    # In this example, we just sample our function f().
    #
    lambdify_numpy_2d = lambda expr: sy.lambdify(("x","y"), expr, modules="numpy")  # SymPy expr --> lambda(x,y)
    f                 = lambdify_numpy_2d(expr)
    dfdx              = lambdify_numpy_2d(sy.diff(expr, "x"))
    dfdy              = lambdify_numpy_2d(sy.diff(expr, "y"))

    print "evaluating example function"
    with SimpleTimer(label=("    done in ")) as s:
        no = wlsqm.number_of_dofs( dimension=2, order=fit_order )
        fi = np.empty( (npoints,no), dtype=np.float64 )
        fi[:,0] = f( S[:,0], S[:,1] )  # fi[i,0] contains the function value at point S[i,:]

    # find the neighborhoods
    print "generating neighborhoods for each point"
    with SimpleTimer(label=("    done in ")) as s:
        hoods = np.zeros( (npoints,max_nk), dtype=np.int32 )  # neighbor point indices (pointing to rows in S[])
        nk    = np.empty( (npoints,), dtype=np.int32 )        # number of neighbors, i.e. nk[i] is the number of actually used columns in hoods[i,:]
        for i in xrange(npoints):
            I = tree.query_ball_point( S[i], r )  # indices of neighbors of S[i] at distance <= r  (but also including S[i] itself!)
            I = filter( lambda idx: idx != i, I )  # exclude S[i] itself
            if len(I) > max_nk:
                I = I[:max_nk]
            I = np.array( I, dtype=np.int32 )
            nk[i] = len(I)
            hoods[i,:nk[i]] = I

    # DEBUG
    print "number of neighbors min = %g, avg = %g, max = %g" % ( np.min(nk), np.mean(nk), np.max(nk) )
    print "neighbor lists for each problem instance:"
    print hoods
    print "number of neighbors for each problem instance:"
    print nk

    # perform the fitting
    print "fitting %d local surrogate models of order %d, driver mode (fit each model once)" % (npoints, fit_order)
    fit_order_array = fit_order        * np.ones( (npoints,), dtype=np.int32 )
    knowns_array    = knowns           * np.ones( (npoints,), dtype=np.int64 )
    wm_array        = weighting_method * np.ones( (npoints,), dtype=np.int32 )
    with SimpleTimer(label=("    done in ")) as s:
#        max_iterations_taken = wlsqm.fit_2D_many( xk=S[hoods], fk=fi[hoods,0], nk=nk,
#                                                   xi=S, fi=fi,
#                                                   sens=None, do_sens=False,
#                                                   order=fit_order_array, knowns=knowns_array, weighting_method=wm_array,
#                                                   debug=False )
        max_iterations_taken = wlsqm.fit_2D_many_parallel( xk=S[hoods], fk=fi[hoods,0], nk=nk,
                                                            xi=S, fi=fi,
                                                            sens=None, do_sens=False,
                                                            order=fit_order_array, knowns=knowns_array, weighting_method=wm_array,
                                                            ntasks=ntasks, debug=False )
#        max_iterations_taken = wlsqm.fit_2D_iterative_many( xk=S[hoods], fk=fi[hoods,0], nk=nk,
#                                                             xi=S, fi=fi,
#                                                             sens=None, do_sens=False,
#                                                             order=fit_order_array, knowns=knowns_array, weighting_method=wm_array,
#                                                             max_iter=max_iter, debug=False )
#        max_iterations_taken = wlsqm.fit_2D_iterative_many_parallel( xk=S[hoods], fk=fi[hoods,0], nk=nk,
#                                                            xi=S, fi=fi,
#                                                            sens=None, do_sens=False,
#                                                            order=fit_order_array, knowns=knowns_array, weighting_method=wm_array,
#                                                            max_iter=max_iter, ntasks=ntasks, debug=False )

    # Expert mode: allows solving multiple times (with new fk data) in the same geometry, performing the prepare step only once.
    #
    # This is especially good for a large number of repetitions with ALGO_BASIC, where a large majority of the computational cost comes from the prepare step.
    #
    # The total advantage is slightly smaller for a small number of repetitions with ALGO_ITERATIVE,
    # since the iterative mode already uses this strategy internally (also when invoked in driver mode).
    #
    print "fitting %d local surrogate models of order %d, expert mode" % (npoints, fit_order)
    print "    init"
    with SimpleTimer(label=("        done in ")) as s:
        solver = wlsqm.ExpertSolver( dimension=2, nk=nk, order=fit_order_array, knowns=knowns_array, weighting_method=wm_array, algorithm=wlsqm.ALGO_BASIC, do_sens=False, max_iter=max_iter, ntasks=ntasks, debug=False )
    print "    prepare"
    with SimpleTimer(label=("        done in ")) as s:
        solver.prepare( xi=S, xk=S[hoods] )
    print "    fit (each model %d times)" % (reps)
    with SimpleTimer(label=("        %d reps done in " % reps), n=reps) as s:
        for k in xrange(reps):
            solver.solve( fk=fi[hoods,0], fi=fi, sens=None )

    # DEBUG
    print "max_iterations_taken: %d" % (max_iterations_taken)
    # see that we got the derivatives at each point
    if fit_order > 0:  # no derivatives if piecewise constant fit
        print dfdx( S[:,0], S[:,1] ) - fi[:,1]
        print dfdy( S[:,0], S[:,1] ) - fi[:,2]

    #########################
    # plotting
    #########################

    xx   = np.linspace(axislims[0], axislims[1], nvis)
    yy   = np.linspace(axislims[2], axislims[3], nvis)<commit_msg>wlsqm_example.py: max_iterations_taken is the number of *corrective* iterations (not including the initial fitting), word this more clearly<commit_after># -*- coding: utf-8 -*-
#
"""Testing script for wlsqm, doubles as a usage example.

-JJ 2016-11-10
"""

from __future__ import division
from __future__ import absolute_import

import time

import numpy as np
import sympy as sy

import scipy.spatial  # cKDTree

import pylab as pl
import mpl_toolkits.mplot3d.axes3d as p3

try:
    import wlsqm
except ImportError:
    print "WLSQM not found; is it installed?"
    from sys import exit
    exit(1)

import sudoku_lhs


# from various scripts, e.g. miniprojects/misc/tworods/main2.py
def axis_marginize(ax, epsx, epsy):
    a = ax.axis()
    w = a[1] - a[0]
    h = a[3] - a[2]
    ax.axis( [ a[0] - w*epsx, a[1] + w*epsx, a[2] - h*epsy, a[3] + h*epsy] )


# from find_neighbors2.py
class SimpleTimer:
    def __init__(self, label="", n=None):
        self.label = label
        self.n     = n      # number of repetitions done inside the "with..." section (for averaging in timing info)

    def __enter__(self):
        self.t0 = time.time()
        return self

    def __exit__(self, errtype, errvalue, traceback):
        dt         = time.time() - self.t0
        identifier = ("%s" % self.label) if len(self.label) else "time taken: "
        avg        = (", avg. %gs per run" % (dt/self.n)) if self.n is not None else ""
        print "%s%gs%s" % (identifier, dt, avg)


# many simultaneous local models, 2D
#
def testmany2d():
    #########################
    # config
    #########################

    ntasks = 8  # OpenMP parallelization

    axislims = [0., 1., 0., 1.]  # [xmin, xmax, ymin, ymax], for plotting
    nvis = 201  # number of visualization points per axis

    expr    = sy.sympify("sin(pi*x) * cos(pi*y)")

    points_per_axis = 100  # for point cloud generation

    r      = 5e-2  # neighborhood radius
    max_nk = 100   # maximum number of neighbor points to accept into each neighborhood (affects memory allocation)

    knowns = 1  # function value is known

    fit_order = 4
    weighting_method = wlsqm.WEIGHT_CENTER

    max_iter = 10  # for iterative fitting method

    reps = 20  # for demonstration of solving multiple times using the same geometry

    #########################
    # the test itself
    #########################

    print
    print "=" * 79
    print "many neighborhoods, 2D case"
    print "=" * 79
    print

    # create a stratified point cloud
    print "generating sudoku sample"
    with SimpleTimer(label=("    done in ")) as s:
        S,m = sudoku_lhs.sample(2, points_per_axis, 1)
        bins_per_axis = m*points_per_axis
        S = S / float(bins_per_axis - 1)  # scale the sample from [0, bins_per_axis-1]**2 to [0, 1]**2
        npoints = len(S)
        print "    %d points" % (npoints)

    # index the point cloud for fast neighbor searching
    print "indexing sample"
    with SimpleTimer(label=("    done in ")) as s:
        tree = scipy.spatial.cKDTree( data=S )

    # If this was an IBVP, we would here get the previous state of the unknown field.
    #
    # In this example, we just sample our function f().
    #
    lambdify_numpy_2d = lambda expr: sy.lambdify(("x","y"), expr, modules="numpy")  # SymPy expr --> lambda(x,y)
    f                 = lambdify_numpy_2d(expr)
    dfdx              = lambdify_numpy_2d(sy.diff(expr, "x"))
    dfdy              = lambdify_numpy_2d(sy.diff(expr, "y"))

    print "evaluating example function"
    with SimpleTimer(label=("    done in ")) as s:
        no = wlsqm.number_of_dofs( dimension=2, order=fit_order )
        fi = np.empty( (npoints,no), dtype=np.float64 )
        fi[:,0] = f( S[:,0], S[:,1] )  # fi[i,0] contains the function value at point S[i,:]

    # find the neighborhoods
    print "generating neighborhoods for each point"
    with SimpleTimer(label=("    done in ")) as s:
        hoods = np.zeros( (npoints,max_nk), dtype=np.int32 )  # neighbor point indices (pointing to rows in S[])
        nk    = np.empty( (npoints,), dtype=np.int32 )        # number of neighbors, i.e. nk[i] is the number of actually used columns in hoods[i,:]
        for i in xrange(npoints):
            I = tree.query_ball_point( S[i], r )  # indices of neighbors of S[i] at distance <= r  (but also including S[i] itself!)
            I = filter( lambda idx: idx != i, I )  # exclude S[i] itself
            if len(I) > max_nk:
                I = I[:max_nk]
            I = np.array( I, dtype=np.int32 )
            nk[i] = len(I)
            hoods[i,:nk[i]] = I

    # DEBUG
    print "number of neighbors min = %g, avg = %g, max = %g" % ( np.min(nk), np.mean(nk), np.max(nk) )
    print "neighbor lists for each problem instance:"
    print hoods
    print "number of neighbors for each problem instance:"
    print nk

    # perform the fitting
    print "fitting %d local surrogate models of order %d, driver mode (fit each model once)" % (npoints, fit_order)
    fit_order_array = fit_order        * np.ones( (npoints,), dtype=np.int32 )
    knowns_array    = knowns           * np.ones( (npoints,), dtype=np.int64 )
    wm_array        = weighting_method * np.ones( (npoints,), dtype=np.int32 )
    with SimpleTimer(label=("    done in ")) as s:
#        max_iterations_taken = wlsqm.fit_2D_many( xk=S[hoods], fk=fi[hoods,0], nk=nk,
#                                                   xi=S, fi=fi,
#                                                   sens=None, do_sens=False,
#                                                   order=fit_order_array, knowns=knowns_array, weighting_method=wm_array,
#                                                   debug=False )
        max_iterations_taken = wlsqm.fit_2D_many_parallel( xk=S[hoods], fk=fi[hoods,0], nk=nk,
                                                            xi=S, fi=fi,
                                                            sens=None, do_sens=False,
                                                            order=fit_order_array, knowns=knowns_array, weighting_method=wm_array,
                                                            ntasks=ntasks, debug=False )
#        max_iterations_taken = wlsqm.fit_2D_iterative_many( xk=S[hoods], fk=fi[hoods,0], nk=nk,
#                                                             xi=S, fi=fi,
#                                                             sens=None, do_sens=False,
#                                                             order=fit_order_array, knowns=knowns_array, weighting_method=wm_array,
#                                                             max_iter=max_iter, debug=False )
#        max_iterations_taken = wlsqm.fit_2D_iterative_many_parallel( xk=S[hoods], fk=fi[hoods,0], nk=nk,
#                                                            xi=S, fi=fi,
#                                                            sens=None, do_sens=False,
#                                                            order=fit_order_array, knowns=knowns_array, weighting_method=wm_array,
#                                                            max_iter=max_iter, ntasks=ntasks, debug=False )

    # Expert mode: allows solving multiple times (with new fk data) in the same geometry, performing the prepare step only once.
    #
    # This is especially good for a large number of repetitions with ALGO_BASIC, where a large majority of the computational cost comes from the prepare step.
    #
    # The total advantage is slightly smaller for a small number of repetitions with ALGO_ITERATIVE,
    # since the iterative mode already uses this strategy internally (also when invoked in driver mode).
    #
    print "fitting %d local surrogate models of order %d, expert mode" % (npoints, fit_order)
    print "    init"
    with SimpleTimer(label=("        done in ")) as s:
        solver = wlsqm.ExpertSolver( dimension=2, nk=nk, order=fit_order_array, knowns=knowns_array, weighting_method=wm_array, algorithm=wlsqm.ALGO_BASIC, do_sens=False, max_iter=max_iter, ntasks=ntasks, debug=False )
    print "    prepare"
    with SimpleTimer(label=("        done in ")) as s:
        solver.prepare( xi=S, xk=S[hoods] )
    print "    fit (each model %d times)" % (reps)
    with SimpleTimer(label=("        %d reps done in " % reps), n=reps) as s:
        for k in xrange(reps):
            solver.solve( fk=fi[hoods,0], fi=fi, sens=None )

    # DEBUG
    print "max corrective iterations taken: %d" % (max_iterations_taken)
    # see that we got the derivatives at each point
    if fit_order > 0:  # no derivatives if piecewise constant fit
        print dfdx( S[:,0], S[:,1] ) - fi[:,1]
        print dfdy( S[:,0], S[:,1] ) - fi[:,2]

    #########################
    # plotting
    #########################

    xx   = np.linspace(axislims[0], axislims[1], nvis)
    yy   = np.linspace(axislims[2], axislims[3], nvis)
{'fraction_non_alphanumeric': 0.1215154085383093, 'fraction_numerical': 0.0121006502685892, 'mean_word_length': 2.039353840866128, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 64, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7283616', 'n_tokens_mistral': 5753, 'n_tokens_neox': 5448, 'n_words': 1596}