File: TasmanianAddons.py

package info (click to toggle)
tasmanian 8.2-2
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 4,852 kB
  • sloc: cpp: 34,523; python: 7,039; f90: 5,080; makefile: 224; sh: 64; ansic: 8
file content (442 lines) | stat: -rw-r--r-- 23,721 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
##############################################################################################################################################################################
# Copyright (c) 2017, Miroslav Stoyanov
#
# This file is part of
# Toolkit for Adaptive Stochastic Modeling And Non-Intrusive ApproximatioN: TASMANIAN
#
# Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions
#    and the following disclaimer in the documentation and/or other materials provided with the distribution.
#
# 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse
#    or promote products derived from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
# INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
# OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
# OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# UT-BATTELLE, LLC AND THE UNITED STATES GOVERNMENT MAKE NO REPRESENTATIONS AND DISCLAIM ALL WARRANTIES, BOTH EXPRESSED AND IMPLIED.
# THERE ARE NO EXPRESS OR IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, OR THAT THE USE OF THE SOFTWARE WILL NOT INFRINGE ANY PATENT,
# COPYRIGHT, TRADEMARK, OR OTHER PROPRIETARY RIGHTS, OR THAT THE SOFTWARE WILL ACCOMPLISH THE INTENDED RESULTS OR THAT THE SOFTWARE OR ITS USE WILL NOT RESULT IN INJURY OR DAMAGE.
# THE USER ASSUMES RESPONSIBILITY FOR ALL LIABILITIES, PENALTIES, FINES, CLAIMS, CAUSES OF ACTION, AND COSTS AND EXPENSES, CAUSED BY, RESULTING FROM OR ARISING OUT OF,
# IN WHOLE OR IN PART THE USE, STORAGE OR DISPOSAL OF THE SOFTWARE.
##############################################################################################################################################################################

from ctypes import c_int, c_double, c_char_p, c_void_p, POINTER, CFUNCTYPE, CDLL, RTLD_GLOBAL
import numpy as np
import sys

import TasmanianConfig
import TasmanianSG
TasmanianInputError = TasmanianConfig.TasmanianInputError

pLibCTSG = CDLL(TasmanianConfig.__path_libcaddons__, mode = RTLD_GLOBAL)

type_1Dfunc = CFUNCTYPE(c_double, c_double)
type_lpnmodel = CFUNCTYPE(None, c_int, POINTER(c_double), c_int, POINTER(c_double), c_int, POINTER(c_int))
type_scsmodel = CFUNCTYPE(None, c_int, c_int, POINTER(c_double), c_int, POINTER(c_double), c_int, POINTER(c_int))
type_icsmodel = CFUNCTYPE(None, c_int, c_int, POINTER(c_double), c_int, c_int, POINTER(c_double), c_int, POINTER(c_int))

pLibCTSG.tsgLoadNeededValues.argtypes = [c_int, type_lpnmodel, c_void_p, c_int, POINTER(c_int)]

pLibCTSG.tsgConstructSurrogateNoIGSurplus.argtypes = [type_scsmodel, c_int, c_int, c_int, c_void_p, c_double, c_char_p,
                                                      c_int, POINTER(c_int), c_char_p, POINTER(c_int)]
pLibCTSG.tsgConstructSurrogateNoIGAniso.argtypes = [type_scsmodel, c_int, c_int, c_int, c_void_p, c_char_p, c_int, POINTER(c_int), c_char_p, POINTER(c_int)]
pLibCTSG.tsgConstructSurrogateNoIGAnisoFixed.argtypes = [type_scsmodel, c_int, c_int, c_int, c_void_p, c_char_p,
                                                         POINTER(c_int), POINTER(c_int), c_char_p, POINTER(c_int)]

pLibCTSG.tsgConstructSurrogateWiIGSurplus.argtypes = [type_icsmodel, c_int, c_int, c_int, c_void_p, c_double, c_char_p,
                                                      c_int, POINTER(c_int), c_char_p, POINTER(c_int)]
pLibCTSG.tsgConstructSurrogateWiIGAniso.argtypes = [type_icsmodel, c_int, c_int, c_int, c_void_p, c_char_p, c_int, POINTER(c_int), c_char_p, POINTER(c_int)]
pLibCTSG.tsgConstructSurrogateWiIGAnisoFixed.argtypes = [type_icsmodel, c_int, c_int, c_int, c_void_p, c_char_p,
                                                         POINTER(c_int), POINTER(c_int), c_char_p, POINTER(c_int)]

pLibCTSG.tsgLoadUnstructuredDataL2.argtypes = [POINTER(c_double), c_int, POINTER(c_double), c_double, c_void_p]

pLibCTSG.tsgCreateExoticQuadratureFromGrid.argtypes = [c_void_p, c_int, c_double, c_void_p, c_char_p, c_int]
pLibCTSG.tsgCreateExoticQuadratureFromFunction.argtypes = [c_void_p, c_int, c_double, type_1Dfunc, c_int, c_char_p, c_int]

def tsgLnpModelWrapper(oUserModel, iSizeX, pX, iSizeY, pY, iThreadID, pErrInfo):
    '''
    DO NOT CALL DIRECTLY
    This is callback from C++, see TasGrid::loadNeededValues()

    Creates an interface between a user callable object and
    the Tasmanian callback from C++.
    The callback passes in two raw-arrays with given sizes
    (that correspond to grid inputs and outputs)
    and the id of the running thread.
    The raw-array is wrapped in a numpy array structure and given
    to the user, the result is written to the output array.

    oUserModel: is callable object, e.g., function or lambda

    iSizeX: number of entries in x, equal to getNumDimensions()
    pX:     raw-array corresponding to the model inputs
    iSizeY: number of entries in y, equal to getNumOutputs()
    pY:     raw-array corresponding to the model outputs

    iThreadID: the id of the running thread
    '''
    pErrInfo[0] = 1
    aX = np.ctypeslib.as_array(pX, (iSizeX,))
    aY = np.ctypeslib.as_array(pY, (iSizeY,))
    aResult = oUserModel(aX, iThreadID)
    if aResult.shape != (iSizeY,):
        if TasmanianConfig.enableVerboseErrors:
            print("ERROR: incorrect model output dimensions, should be (iNumOutputs,)")
        return
    aY[0:iSizeY] = aResult[0:iSizeY]
    pErrInfo[0] = 0


def loadNeededValues(callableModel, grid, iNumThreads = 1):
    '''
    Wrapper to TasGrid::loadNeededValues(), non-overwrite version.

    If the grid has needed points, the callableModel will be called
    for each grid point (i.e., model input) and the resulting values
    will be loaded in the grid.

    callableModel: is callable object, e.g., function or lambda
        The object must accept two inputs and give one output:
        aY = callableModel(aX, iThreadID)

        aX: is a one dimensional numpy.ndarray with size equal
            to the number of model inputs

        iThreadID: is the ID of the thread executing the model,
            always between 0 and iNumThreads -1
            Two simultaneous calls to callableModel() will always
            have different ids.

        Return: aY must be a one dimensional numpy.ndarray with
            size equal to the number of model outputs

        Note: if iNumThreads > 1, then callableModel() must be thread-safe.

    grid: must be an instance of Tasmanian.SparseGrid()
        model values will be loaded in the grid

    iNumThreads: integer, if greater than 1 the model will be called
        in parallel from multiple threads.
        See TasGrid::loadNeededValues().

    '''
    iOverwrite = 0 # do not overwrite
    pErrorCode = (c_int * 1)()
    pLibCTSG.tsgLoadNeededValues(iOverwrite,
                                 type_lpnmodel(lambda nx, x, ny, y, tid, err : tsgLnpModelWrapper(callableModel, nx, x, ny, y, tid, err)),
                                 grid.pGrid, iNumThreads, pErrorCode)
    if pErrorCode[0] != 0:
        raise TasmanianInputError("loadNeededValues", "An error occurred during the call to Tasmanian.")

def loadNeededPoints(callableModel, grid, iNumThreads = 1):
    '''
    Alias to loadNeededValues()
    '''
    loadNeededValues(callableModel, grid, iNumThreads)

def reloadLoadedValues(callableModel, grid, iNumThreads = 1):
    '''
    Wrapper to TasGrid::loadNeededPoints(), overwrite version.

    Clears any pending refinement (i.e., needed points) and overwrites the model
    values associated with the existing loaded points.

    The inputs are identical to Tasmanian.loadNeededPoints().

    '''
    iOverwrite = 1 # do overwrite
    pErrorCode = (c_int * 1)()
    pLibCTSG.tsgLoadNeededValues(iOverwrite,
                                 type_lpnmodel(lambda nx, x, ny, y, tid, err : tsgLnpModelWrapper(callableModel, nx, x, ny, y, tid, err)),
                                 grid.pGrid, iNumThreads, pErrorCode)
    if pErrorCode[0] != 0:
        raise TasmanianInputError("reloadLoadedPoints()", "An error occurred during the call to Tasmanian.")

def reloadLoadedPoints(callableModel, grid, iNumThreads = 1):
    '''
    Alias to reloadLoadedValues()
    '''
    reloadLoadedValues(callableModel, grid, iNumThreads)


###############################################################################
################### Construct Surrogate #######################################
###############################################################################

def tsgScsModelWrapper(oUserModel, iNumSamples, iNumDims, pX, iNumOuts, pY, iThreadID, pErrInfo):
    '''
    DO NOT CALL DIRECTLY
    This is callback from C++, see TasGrid::constructSurrogate()

    Handles the case of batch models:
    oUserModel: user defined model that takes a two dimensional
        array of inputs and returns a two dimensional array
        of outputs, and a thread ID

    iNumSamples: number of samples in the batch
    iNumDims:    number of model inputs per sample
    iNumOuts:    number of model outputs per sample

    pX and pY are 2D arrays of c_doubles with size iNumSamples
        times iNumDims and iNumOuts respectively

    iThreadID: the id of the running thread

    '''
    pErrInfo[0] = 1
    aX = np.ctypeslib.as_array(pX, (iNumSamples,iNumDims))
    aY = np.ctypeslib.as_array(pY, (iNumSamples,iNumOuts))
    aResult = oUserModel(aX, iThreadID)
    if aResult.shape != (iNumSamples, iNumOuts):
        if TasmanianConfig.enableVerboseErrors:
            print("ERROR: incorrect model output dimensions, should be (iNumSamples, iNumOutputs)")
        return
    aY[0:iNumSamples, 0:iNumOuts] = aResult[0:iNumSamples, 0:iNumOuts]
    pErrInfo[0] = 0

def tsgIcsModelWrapper(oUserModel, iNumSamples, iNumDims, pX, iHasGuess, iNumOuts, pY, iThreadID, pErrInfo):
    '''
    DO NOT CALL DIRECTLY
    This is callback from C++, see TasGrid::constructSurrogate()

    See tsgScsModelWrapper(), the only difference is that
    the user model oUserModel() takes two arrays,
    one with the inputs and one with the initial guess.
    The initial guess could be empty.

    iHasGuess is a boolean that determines whether an
    initial guess has been loaded in pY.
    '''
    pErrInfo[0] = 1
    aX = np.ctypeslib.as_array(pX, (iNumSamples,iNumDims))
    aY = np.ctypeslib.as_array(pY, (iNumSamples,iNumOuts))
    if (iHasGuess == 0): # no guess
        aResult = oUserModel(aX, np.empty([0,0], np.float64), iThreadID)
    else:
        aResult = oUserModel(aX, aY, iThreadID)
    if aResult.shape != (iNumSamples, iNumOuts):
        if TasmanianConfig.enableVerboseErrors:
            print("ERROR: incorrect model output dimensions, should be (iNumSamples, iNumOutputs)")
        return
    aY[0:iNumSamples, 0:iNumOuts] = aResult[0:iNumSamples, 0:iNumOuts]
    pErrInfo[0] = 0

def constructAnisotropicSurrogate(callableModel, iMaxPoints, iMaxParallel, iMaxSamplesPerCall, grid,
                                  sDepthType, liAnisotropicWeightsOrOutput,
                                  liLevelLimits = [], bUseInitialGuess = False,
                                  sCheckpointFilename = ""):
    '''
    Construct a surrogate model to the callableModel using
    anisotropic refinement until the iMaxPoints is reached.

    See the documentation for TasGrid::constructSurrogate()
    This is wrapper around the anisotropic refinement variant.

    callableModel is a function (or lambda) that returns a two dimensional
        numpy.ndarray of outputs for the model,
        If bUseInitialGuess is true, the model has to take two inputs,
        otherwise the model uses one.
        The first input is a two dimensional numpy.ndarray of inputs
        similar to TasmanianSparseGrid.evaluateBatch().
        The second input corresponds to the initial guess, the size will
        either match the expected output or will be empty, if no
        guess can be computed.

    iMaxPoints: is a positive integer indicating the maximum number of points
        that the grid will have.
    iMaxParallel: is a positive integer indicating the number of simultaneous
        calls to the user model, i.e., the number of threads.
    iMaxSamplesPerCall: maximum number of samples that will be given to
        a single call to the user model.
    grid: must be an instance of Tasmanian.SparseGrid() with either
        global, sequence or Fourier grid.
    sDepthType: the type used for refinement, see TasGrid::constructSurrogate()

    liAnisotropicWeightsOrOutput: is either an output to use to determine
        the model anisotropy or a list/typle/ndarray of user selected
        anisotropic weights.
    liLevelLimits: same as in all other refinement calls, the refinement will
        never add points below the given level in the diven direction
        even if the budget has not been reached yet.

    sCheckpointFilename: filename to use to checkpoint the algorithm so that
        construction can proceed from a saved point in case of a crash.
    '''
    iNumDims = grid.getNumDimensions()
    pLevelLimits = None
    if (len(liLevelLimits) > 0):
        if (len(liLevelLimits) != iNumDims):
            raise TasmanianInputError("liLevelLimits", "ERROR: invalid number of level limits, must be equal to the grid dimension")
        pLevelLimits = (c_int*iNumDims)()
        for iI in range(iNumDims):
            pLevelLimits[iI] = liLevelLimits[iI]

    if (sys.version_info.major == 3):
        sDepthTypeCtypes = bytes(sDepthType, encoding='utf8')
        if (sCheckpointFilename):
            sCheckpointFilename = bytes(sCheckpointFilename, encoding='utf8')
    else:
        sDepthTypeCtypes = sDepthType

    pCPFname = None
    if (sCheckpointFilename):
        pCPFname = c_char_p(sCheckpointFilename)

    pErrorCode = (c_int * 1)()

    if (((sys.version_info.major == 3) and isinstance(liAnisotropicWeightsOrOutput, int))
            or ((sys.version_info.major == 2) and isinstance(liAnisotropicWeightsOrOutput, (int, long)))):
        # will call the algorithm to dynamically estimate the weights
        iOutput = liAnisotropicWeightsOrOutput

        if (bUseInitialGuess):
            pLibCTSG.tsgConstructSurrogateWiIGAniso(
                type_icsmodel(lambda nx, nd, x, f, ny, y, tid, err : tsgIcsModelWrapper(callableModel, nx, nd, x, f, ny, y, tid, err)),
                iMaxPoints, iMaxParallel, iMaxSamplesPerCall, grid.pGrid,
                c_char_p(sDepthTypeCtypes), iOutput, pLevelLimits, pCPFname, pErrorCode)
        else:
            pLibCTSG.tsgConstructSurrogateNoIGAniso(
                type_scsmodel(lambda nx, nd, x, ny, y, tid, err : tsgScsModelWrapper(callableModel, nx, nd, x, ny, y, tid, err)),
                iMaxPoints, iMaxParallel, iMaxSamplesPerCall, grid.pGrid,
                c_char_p(sDepthTypeCtypes), iOutput, pLevelLimits, pCPFname, pErrorCode)
    else:
        # weights are set by the user
        pAnisoWeights = None
        if (len(liAnisotropicWeightsOrOutput) > 0):
            if (sDepthType in TasmanianSG.lsTsgCurvedTypes):
                iNumWeights = 2*grid.getNumDimensions()
            else:
                iNumWeights = grid.getNumDimensions()
            if (len(liAnisotropicWeightsOrOutput) != iNumWeights):
                raise TasmanianInputError("liAnisotropicWeightsOrOutput", "ERROR: wrong number of liAnisotropicWeightsOrOutput, sType '{0:s}' needs {1:1d} weights but len(liAnisotropicWeightsOrOutput) == {2:1d}".format(sDepthType, iNumWeights, len(liAnisotropicWeightsOrOutput)))
            else:
                aAWeights = np.array([liAnisotropicWeightsOrOutput[i] for i in range(iNumWeights)], np.int32)
                pAnisoWeights = np.ctypeslib.as_ctypes(aAWeights)

        if (bUseInitialGuess):
            pLibCTSG.tsgConstructSurrogateWiIGAnisoFixed(
                type_icsmodel(lambda nx, nd, x, f, ny, y, tid, err : tsgIcsModelWrapper(callableModel, nx, nd, x, f, ny, y, tid, err)),
                iMaxPoints, iMaxParallel, iMaxSamplesPerCall, grid.pGrid,
                c_char_p(sDepthTypeCtypes), pAnisoWeights, pLevelLimits, pCPFname, pErrorCode)
        else:
            pLibCTSG.tsgConstructSurrogateNoIGAnisoFixed(
                type_scsmodel(lambda nx, nd, x, ny, y, tid, err : tsgScsModelWrapper(callableModel, nx, nd, x, ny, y, tid, err)),
                iMaxPoints, iMaxParallel, iMaxSamplesPerCall, grid.pGrid,
                c_char_p(sDepthTypeCtypes), pAnisoWeights, pLevelLimits, pCPFname, pErrorCode)

    if pErrorCode[0] != 0:
        raise TasmanianInputError("constructSurplusSurrogate", "An error occurred during the call to Tasmanian.")


def constructSurplusSurrogate(callableModel, iMaxPoints, iMaxParallel, iMaxSamplesPerCall, grid,
                              fTolerance, sRefinementType, iOutput = -1,
                              liLevelLimits = [], bUseInitialGuess = False,
                              sCheckpointFilename = ""):
    '''
    Construct a surrogate model to the callableModel using surplus refinement
    until either the iMaxPoints or the tolerance are reached.

    See Tasmanian.constructAnisotropicSurrogate() for all matchin inputs,
    except the grid has to be local polynomial and the refinement proceeds
    until the budget is exhausted or the fTolerance is reached.
    The sRefinementType is the same as in the call to local surplus refinement,
    same with the iOutput.
    '''
    iNumDims = grid.getNumDimensions()
    pLevelLimits = None
    if (len(liLevelLimits) > 0):
        if (len(liLevelLimits) != iNumDims):
            raise TasmanianInputError("liLevelLimits", "ERROR: invalid number of level limits, must be equal to the grid dimension")
        pLevelLimits = (c_int*iNumDims)()
        for iI in range(iNumDims):
            pLevelLimits[iI] = liLevelLimits[iI]

    pCPFname = bytes(sCheckpointFilename, encoding='utf8') if sCheckpointFilename else None

    pErrorCode = (c_int * 1)()

    if (bUseInitialGuess):
        pLibCTSG.tsgConstructSurrogateWiIGSurplus(
            type_icsmodel(lambda nx, nd, x, f, ny, y, tid, err : tsgIcsModelWrapper(callableModel, nx, nd, x, f, ny, y, tid, err)),
            iMaxPoints, iMaxParallel, iMaxSamplesPerCall, grid.pGrid,
            fTolerance, bytes(sRefinementType, encoding='utf8'), iOutput, pLevelLimits, pCPFname, pErrorCode)
    else:
        pLibCTSG.tsgConstructSurrogateNoIGSurplus(
            type_scsmodel(lambda nx, nd, x, ny, y, tid, err : tsgScsModelWrapper(callableModel, nx, nd, x, ny, y, tid, err)),
            iMaxPoints, iMaxParallel, iMaxSamplesPerCall, grid.pGrid,
            fTolerance, bytes(sRefinementType, encoding='utf8'), iOutput, pLevelLimits, pCPFname, pErrorCode)

    if pErrorCode[0] != 0:
        raise TasmanianInputError("constructSurplusSurrogate", "An error occurred during the call to Tasmanian.")


def loadUnstructuredDataL2(points, model_data, tolerance, grid):
    '''
    Wrapper around TasGrid::loadUnstructuredDataL2(), see the C++ reference.

    points is 2D ndarray with shape(num_data, grid.getNumDimensions())
    model_data is 2D ndarray with shape(num_data, grid.getNumOutputs())
    grid is an instance of TasmanianSparseGrid
    '''
    if len(points.shape) != 2 or points.shape[1] != grid.getNumDimensions():
        raise TasmanianInputError("points", "ERROR: points must be a 2D numpy.ndarray with points.shape[1] == grid.getNumDimensions()")
    if len(model_data.shape) != 2 or model_data.shape[1] != grid.getNumOutputs():
        raise TasmanianInputError("model_data", "ERROR: model_data must be a 2D numpy.ndarray with model_data.shape[1] == grid.getNumOutputs()")
    if points.shape[0] != model_data.shape[0]:
        raise TasmanianInputError("model_data", "ERROR: mismatch between shape[0] of points and model_data")
    if not hasattr(grid, "TasmanianSparseGridObject"):
        raise TasmanianInputError("grid", "ERROR: grid must be an instance of TasmanianSparseGrid")

    num_data = points.shape[0]
    pLibCTSG.tsgLoadUnstructuredDataL2(np.ctypeslib.as_ctypes(points.reshape((points.size,))), num_data,
                                       np.ctypeslib.as_ctypes(model_data.reshape((model_data.size,))), tolerance, grid.pGrid)

def createExoticQuadratureFromGrid(level, shift, ref_grid, description, is_symmetric = False):
    '''
    Calls TasGrid::getExoticQuadrature() from a one dimensional interpolant/surrogate of the weight function, and output a
    python CustomTabulated object.
    See the C++ reference for more information.

    level:        positive integer representing the level of the exotic quadrature grid.
    shift:        double where [weight_function(x) + shift] is nonegative for every x in the domain of the weight function.
    ref_grid:     Python TasmanianSparseGrid object that represents the one dimensional surrogate/interpolant of the weight function.
    description:  string describing the Exotic quadrature instance.
    is_symmetric: (optional) boolean that should be set to True if the weight function is symmetric.

    output:       a Python CustomTabulated object.
    '''
    if not hasattr(ref_grid, "TasmanianSparseGridObject"):
        raise TasmanianInputError("ref_grid", "ERROR: ref_grid must be an instance of TasmanianSparseGrid")
    ct = TasmanianSG.CustomTabulated()
    pLibCTSG.tsgCreateExoticQuadratureFromGrid(c_void_p(ct.pCustomTabulated), c_int(level), c_double(shift), c_void_p(ref_grid.pGrid),
                                               bytes(description, encoding='utf8'), c_int(is_symmetric))
    return ct

def createExoticQuadratureFromFunction(level, shift, weight_fn, nref, description, is_symmetric = False):
    '''
    Calls TasGrid::getExoticQuadrature() from a function lambda representing the weight function, and output a Python
    CustomTabulated object.
    See the C++ reference for more information.

    level:        positive integer representing the level of the exotic quadrature grid.
    shift:        double where [weight_function(x) + shift] is nonegative for every x in the domain of the weight function.
    weight_fn:    Python lambda function representing the weight function.
    nref:         positive integer representing the number of points used to generate the weight function surrogate/interpolant.
    description:  string describing the Exotic quadrature instance.
    is_symmetric: (optional) boolean that should be set to True if the weight function is symmetric.

    output:       a Python CustomTabulated object.
    '''
    ct = TasmanianSG.CustomTabulated()
    pLibCTSG.tsgCreateExoticQuadratureFromFunction(c_void_p(ct.pCustomTabulated), c_int(level), c_double(shift), type_1Dfunc(weight_fn),
                                                   c_int(nref), bytes(description, encoding='utf8'), c_int(is_symmetric))
    return ct