LibLip Maple toolbox 

School of Engineering and IT, Deakin University, Australia 

Copyright Gleb Beliakov, 2006 

Australia 

gbeliako@gmail.com 

Installation instructions 

To use this worksheet you need to have MapleLibLip.dll library in the curent directory or on the path 

Introduction 

This worksheet illustrates how LibLip package can be used with Maple (under MS Windows).   

LibLip is a programming library for interpolation and smoothing of multivariate scattered data with Lipschitz functions. 

 

Please see the user manual supplied with this worksheet. 

This worksheet is provided as a template for user's own calculations. It shows how to declare the data and  

supply it to LibLip methods. 

It should be used in conjunction with the supplied dll (MapleLibLip.dll), and the user manual (liblipmaple.pdf). 

The typical calling sequence is as follows: 

1. execute define_external procedures in the next subsection 

2. Define your data 

3. call the relevant  interpolaiton routines with correct arguments
 

Please change the LIBLIPPATH line below to the location of this package on your computer! 

Initialization: Declarations of LibLip procedures found in the dll 

>	restart; with(plots):

 

Execute the commands below to declare the external subroutines from LibLip library. 

>

DLLPATH:="mapleliblip.dll":    # if you need to provide full path, use, e.g.,DLLPATH:="c:/work/maple/mapleliblip.dll": LipIntValue := define_external('MWRAP_LipIntValue','MAPLE', LIB=DLLPATH): LipIntValueAuto := define_external('MWRAP_LipIntValueAuto','MAPLE', LIB=DLLPATH): LipIntValueCons := define_external('MWRAP_LipIntValueCons','MAPLE', LIB=DLLPATH): LipIntValueConsLeftRegion := define_external('MWRAP_LipIntValueConsLeftRegion','MAPLE', LIB=DLLPATH): LipIntValueConsRightRegion := define_external('MWRAP_LipIntValueConsRightRegion','MAPLE', LIB=DLLPATH): LipIntValueLocal := define_external('MWRAP_LipIntValueLocal','MAPLE', LIB=DLLPATH): LipIntValueLocalCons := define_external('MWRAP_LipIntValueLocalCons','MAPLE', LIB=DLLPATH): LipIntValueLocalConsLeftRegion := define_external('MWRAP_LipIntValueLocalConsLeftRegion','MAPLE', LIB=DLLPATH): LipIntValueLocalConsRightRegion := define_external('MWRAP_LipIntValueLocalConsRightRegion','MAPLE', LIB=DLLPATH): LipIntComputeLipschitz := define_external('MWRAP_LipIntComputeLipschitz','MAPLE', LIB=DLLPATH): LipIntComputeLocalLipschitz := define_external('MWRAP_LipIntComputeLocalLipschitz','MAPLE', LIB=DLLPATH): LipIntComputeLipschitzCV := define_external('MWRAP_LipIntComputeLipschitzCV','MAPLE', LIB=DLLPATH): LipIntComputeLipschitzSplit := define_external('MWRAP_LipIntComputeLipschitzSplit','MAPLE', LIB=DLLPATH): LipIntSmoothLipschitz := define_external('MWRAP_LipIntSmoothLipschitz','MAPLE', LIB=DLLPATH): LipIntGetLipConst := define_external('MWRAP_LipIntGetLipConst','MAPLE', LIB=DLLPATH): LipIntGetScaling := define_external('MWRAP_LipIntGetScaling','MAPLE', LIB=DLLPATH): LipIntComputeScaling := define_external('MWRAP_LipIntComputeScaling','MAPLE', LIB=DLLPATH): LipIntVerifyMonotonicity := define_external('MWRAP_LipIntVerifyMonotonicity','MAPLE', LIB=DLLPATH): LipIntVerifyMonotonicityLeftRegion := define_external('MWRAP_LipIntVerifyMonotonicityLeftRegion','MAPLE', LIB=DLLPATH): LipIntVerifyMonotonicityRightRegion := define_external('MWRAP_LipIntVerifyMonotonicityRightRegion','MAPLE', LIB=DLLPATH): LipIntInfValue := define_external('MWRAP_LipIntInfValue','MAPLE', LIB=DLLPATH): LipIntInfValueAuto := define_external('MWRAP_LipIntInfValueAuto','MAPLE', LIB=DLLPATH): LipIntInfValueCons := define_external('MWRAP_LipIntInfValueCons','MAPLE', LIB=DLLPATH): LipIntInfValueConsLeftRegion := define_external('MWRAP_LipIntInfValueConsLeftRegion','MAPLE', LIB=DLLPATH): LipIntInfValueConsRightRegion := define_external('MWRAP_LipIntInfValueConsRightRegion','MAPLE', LIB=DLLPATH): LipIntInfValueLocal := define_external('MWRAP_LipIntInfValueLocal','MAPLE', LIB=DLLPATH): LipIntInfValueLocalCons := define_external('MWRAP_LipIntInfValueLocalCons','MAPLE', LIB=DLLPATH): LipIntInfValueLocalConsLeftRegion := define_external('MWRAP_LipIntInfValueLocalConsLeftRegion','MAPLE', LIB=DLLPATH): LipIntInfValueLocalConsRightRegion := define_external('MWRAP_LipIntInfValueLocalConsRightRegion','MAPLE', LIB=DLLPATH): LipIntInfComputeLipschitz := define_external('MWRAP_LipIntInfComputeLipschitz','MAPLE', LIB=DLLPATH): LipIntInfComputeLocalLipschitz := define_external('MWRAP_LipIntInfComputeLocalLipschitz','MAPLE', LIB=DLLPATH): LipIntInfComputeLipschitzCV := define_external('MWRAP_LipIntInfComputeLipschitzCV','MAPLE', LIB=DLLPATH): LipIntInfComputeLipschitzSplit := define_external('MWRAP_LipIntInfComputeLipschitzSplit','MAPLE', LIB=DLLPATH): LipIntInfSmoothLipschitz := define_external('MWRAP_LipIntInfSmoothLipschitz','MAPLE', LIB=DLLPATH): LipIntInfGetLipConst := define_external('MWRAP_LipIntInfGetLipConst','MAPLE', LIB=DLLPATH): LipIntInfGetScaling := define_external('MWRAP_LipIntInfGetScaling','MAPLE', LIB=DLLPATH): LipIntInfComputeScaling := define_external('MWRAP_LipIntInfComputeScaling','MAPLE', LIB=DLLPATH): LipIntInfVerifyMonotonicity := define_external('MWRAP_LipIntInfVerifyMonotonicity','MAPLE', LIB=DLLPATH): LipIntInfVerifyMonotonicityLeftRegion := define_external('MWRAP_LipIntInfVerifyMonotonicityLeftRegion','MAPLE', LIB=DLLPATH): LipIntInfVerifyMonotonicityRightRegion := define_external('MWRAP_LipIntInfVerifyMonotonicityRightRegion','MAPLE', LIB=DLLPATH): LipIntSetBounds := define_external('MWRAP_SetBoundaries','MAPLE', LIB=DLLPATH): LipIntClearBounds := define_external('MWRAP_ClearBoundaries','MAPLE', LIB=DLLPATH): STCBuildLipInterpolant := define_external('MWRAP_STCBuildLipInterpolant','MAPLE', LIB=DLLPATH): STCBuildLipInterpolantExplicit := define_external('MWRAP_STCBuildLipInterpolantExplicit','MAPLE', LIB=DLLPATH): STCValue := define_external('MWRAP_STCValue','MAPLE', LIB=DLLPATH): STCValueExplicit := define_external('MWRAP_STCValueExplicit','MAPLE', LIB=DLLPATH): STCFreeMemory := define_external('MWRAP_STCFreeMemory','MAPLE', LIB=DLLPATH):

 

Section 1: Generation of the data and its interpolation 

We need to declare two arrays: abscissae XD and function values YD. We will generate the data using some test function (you may change it) 

Declaration of the function f(x) that generates the data 

Declare f as a procedure.  

>	f := proc( x, y)  # user's density function        exp(-(y-x^2)^2 - (x^2+y^2)/2);    end proc;

 

 

>	# test the density function:  evaluate f(x), plot it `f(0.2,0.6)=`; f(0.2,0.6); plot3d(f(x,y),x=-2..2,y=-1..2, axes=boxed,orientation=[25,50]);

 

 

 

 

Generating the data (random abscissae, uniformly distributed in the domain) 

Boundaries, and also the dimension of x and declaration of various vectors.  

>

Lo:=Vector(1..2,datatype=float[8],[-2,-1]):   Up:=Vector(1..2,datatype=float[8],[2,2]): dim:=2: Ndata:=600: Xr:=[stats[random, uniform]((Ndata)*dim)]: # uniformly distributed numbers XD:=Matrix(1..Ndata, 1..dim, datatype=float[8], order=C_order): # will hold the data abscissae YD:=Vector(1..Ndata,datatype=float[8], order=C_order): # will hold the function values X:=Vector(1..dim,datatype=float[8]): scale:=(x,a,b)->x*(b-a)+a: # just a scaling function for i from 1 to Ndata do for j from 1 to dim do   X[j]:=scale(Xr[i*j],Lo[j],Up[j]): od: YD[i]:= f(X[1],X[2]):        # record function values for j from 1 to dim do   XD[i,j]:=X[j]:             # record data abscisae od: od:

 

>	pointplot({seq([XD[k,1],XD[k,2]],k=1..Ndata)}); # this is how the data are distributed

 

 

>	pointplot3d({seq([XD[k,1],XD[k,2],YD[k]],k=1..Ndata)},axes=BOXED,orientation=[25,50]);

 

 

Executing LibLip methods to approximate function values 

Interpolating the data: choose any random vector X and calculate the value of f 

>

LC:=10: lprint(`approximation       \t \t exact value     \t error`); for i from 1 to 10 do  Xtest:=[stats[random, uniform](dim)]:   # uniformly distributed numbers  for j from 1 to dim do   X[j]:=scale(Xtest[j],Lo[j],Up[j]):  od:  fint:=LipIntValue(dim,Ndata,X,XD,YD,LC):  fexact:=f(X[1],X[2]):  lprint(fint,`    \t`,fexact,`    \t`,abs(fint-fexact)); od:

 

`approximation        exact value      error` .517837017099999964, `     `, .5140309449, `     `, 0.38060722e-2 0.793549996800035907e-3, `     `, 0.2822637558e-3, `     `, 0.5112862410e-3 0.858336488551714182e-7, `     `, 0.1188342648e-7, `     `, 0.7395022238e-7 0.452840224099999578e-1, `     `, 0.3301228124e-1, `     `, 0.1227174117e-1 0.155828984200023780e-3, `     `, 0.1166064469e-3, `     `, 0.392225373e-4 0.901900033200000173e-1, `     `, 0.7685018221e-1, `     `, 0.1333982111e-1 .123249067172205573, `     `, .1188084203, `     `, 0.44406469e-2 .163830715800000026, `     `, .1463113837, `     `, 0.175193321e-1 .880383238700000126, `     `, .9070838534, `     `, 0.267006147e-1 0.403472440399998966e-2, `     `, 0.6573609521e-2, `     `, 0.2538885117e-2

 

Compute the Lipschitz constant from the data 

>	LipIntComputeLipschitz(dim, Ndata, XD,YD); LC := LipIntGetLipConst();

 

 

Plot the interpolant together with the data. We define an auxiliary function for this 

>	f_aux:=(t1,t2)->LipIntValue(dim,Ndata, convert([t1,t2],Vector,datatype=float[8]), XD,YD,LC): f_aux(0.4,0.2);  # test its usage

 

 

>

R:=20.0: # plot resolutions: determines how many points to plot sdata :=[seq([ seq([i/R,j/R,f_aux(i/R,j/R)], i=R*Lo[1]..R*Up[1])], j=R*Lo[2]..R*Up[2])]: S1:=surfdata(sdata,axes=normal,grid=[100,100],style=HIDDEN, axes=BOXED,orientation=[25,50]): S2:=pointplot3d({seq([XD[k,1],XD[k,2],YD[k]],k=1..Ndata)},axes=BOXED,color=BLACK): display({S1,S2});

 

 

Now use a different method of locally Lipschitz functions 

>	LipIntComputeLocalLipschitz(dim,Ndata,XD,YD); f_aux1:=(t1,t2)->LipIntValueLocal(dim,Ndata, convert([t1,t2],Vector,datatype=float[8]), XD,YD): f_aux1(0.1,0.2);  #test the auxiliary function

 

 

>

R:=20.0: #plot resolution sdata :=[seq([ seq([i/R,j/R,f_aux1(i/R,j/R)], i=R*Lo[1]..R*Up[1])], j=R*Lo[2]..R*Up[2])]: S1:=surfdata(sdata,axes=normal,grid=[100,100],style=HIDDEN, axes=BOXED,orientation=[25,50]): S2:=pointplot3d({seq([XD[k,1],XD[k,2],YD[k]],k=1..Ndata)},axes=BOXED,color=BLACK): display({S1,S2});

 

 

The next section gives an example of smoothing noisy data 

Smoothing noisy data 

We shall use the same data with some added uniformly distributed noise 

Add noise to the data 

>	Noiselevel:=0.05: Noise:=[stats[random, uniform[-1,1]](Ndata)]: # uniformly distributed numbers

 

>	for i from 1 to Ndata do  YD[i]:=YD[i]+Noiselevel*Noise[i]: od:

 

Let us plot the interpolant of the noisy data, to see that it is really inadequate 

>

LipIntComputeLipschitz(dim, Ndata, XD,YD); LC := LipIntGetLipConst(); R:=20.0: #plot resolution sdata :=[seq([ seq([i/R,j/R,f_aux(i/R,j/R)], i=R*Lo[1]..R*Up[1])], j=R*Lo[2]..R*Up[2])]: S1:=surfdata(sdata,axes=normal,grid=[100,100],style=HIDDEN, axes=BOXED,orientation=[25,50]): S2:=pointplot3d({seq([XD[k,1],XD[k,2],YD[k]],k=1..Ndata)},axes=BOXED,color=BLACK): display({S1,S2});

 

 

 

Smoothing using a desired Lipschitz constant 

>	LC:=1.5:  # this is a desired value TD:=Vector(1..Ndata,datatype=float[8], order=C_order): # will hold the smoothened values

 

>	LipIntSmoothLipschitz(dim, Ndata,XD,YD,TD,LC,0,0,0);

 

Plot interpolant of the smoothened data 

Notice that we use TD (the array that contains smoothened data) rather than YD 

>	f_aux2:=(t1,t2)->LipIntValue(dim,Ndata, convert([t1,t2],Vector,datatype=float[8]), XD,TD,LC): f_aux2(0.4,0.2);  #test it

 

 

>

R:=20.0: #plot resolution sdata :=[seq([ seq([i/R,j/R,f_aux2(i/R,j/R)], i=R*Lo[1]..R*Up[1])], j=R*Lo[2]..R*Up[2])]: S1:=surfdata(sdata,axes=normal,grid=[100,100],style=HIDDEN, axes=BOXED,orientation=[25,50]): S2:=pointplot3d({seq([XD[k,1],XD[k,2],YD[k]],k=1..Ndata)},axes=BOXED,color=BLACK): display({S1,S2});

 

 

Next we determine the Lipschitz constant automatically by using sampe splitting 

Determining the Lipschitz constant 

>	LipIntComputeLipschitzSplit(dim,Ndata,XD,YD,TD,0.2,0); #LipIntComputeLipschitzCV(dim,40,XD,YD,TD,0); `Estimated Lipschitz constant from noisy data`;LC := LipIntGetLipConst();

 

 

 

The estimation routine has already computed the Lipschitz constant, and has smoothened all data (it returned array TD) 

>

f_aux2:=(t1,t2)->LipIntValue(dim,Ndata, convert([t1,t2],Vector,datatype=float[8]), XD,TD,LC): R:=20.0: #plot resolution sdata :=[seq([ seq([i/R,j/R,f_aux2(i/R,j/R)], i=R*Lo[1]..R*Up[1])], j=R*Lo[2]..R*Up[2])]: S1:=surfdata(sdata,axes=normal,grid=[100,100],style=HIDDEN, axes=BOXED,orientation=[25,50]): S2:=pointplot3d({seq([XD[k,1],XD[k,2],YD[k]],k=1..Ndata)},axes=BOXED,color=BLACK): display({S1,S2});

 

 

Using very fast interpolation method with preprocessing, with much more data 

>

Ndata:=20000: Xr:=[stats[random, uniform]((Ndata)*dim)]: # uniformly distributed numbers XD:=Matrix(1..Ndata, 1..2, datatype=float[8], order=C_order): # will hold the data abscissae YD:=Vector(1..Ndata,datatype=float[8], order=C_order): # will hold the function values X:=Vector(1..2,datatype=float[8]): for i from 1 to Ndata do for j from 1 to dim do   X[j]:=scale(Xr[i*j],Lo[j],Up[j]): od: YD[i]:= f(X[1],X[2]):        # record function values for j from 1 to dim do   XD[i,j]:=X[j]:             # record data abscisae od: od:

 

>	STCFreeMemory(); # just in case we are repeating STCBuildLipInterpolant(dim,Ndata,XD,YD); # nonzero result indicates an error

 

 

Observe the speed with which this method generates the surface. If you change STCValue to STCValueExplicit in the libe below, the calculation time will be much higher 

>	f_aux3:=(t1,t2)->STCValue(convert([t1,t2],Vector,datatype=float[8])):

 

>

R:=20.0: # plot resolution sdata :=[seq([ seq([i/R,j/R,f_aux3(i/R,j/R)], i=R*Lo[1]..R*Up[1])], j=R*Lo[2]..R*Up[2])]: S1:=surfdata(sdata,axes=normal,grid=[100,100],style=HIDDEN, axes=BOXED,orientation=[25,50]): #S2:=pointplot3d({seq([XD[k,1],XD[k,2],YD[k]],k=1..Ndata)},axes=BOXED,color=BLACK): display({S1});

 

>

 

 

Section 2: Mootone interpolation 

We illustrate monotone interpolation on an example of a monotone increasing function in x1 and decreasing in x2 

Generate data by using a monotone function  

>	g := proc( x, y)  # user's test function        min(0.5,x*(1-y)/(2-x+y-x*y));    end proc:

 

>	plot3d(g(x,y),x=0..1,y=0..1, axes=boxed,orientation=[-140,60]);

 

 

The commands below just generate the data to interpolate 

>

Lo:=Vector(1..2,datatype=float[8],[0,0]):   Up:=Vector(1..2,datatype=float[8],[1,1]): Cons:=Vector(1..2,datatype=integer[4],[1,-1]): dim:=2: Ndata:=100: Xr:=[stats[random, uniform]((Ndata)*dim)]: # uniformly distributed numbers Noise:=[stats[random, uniform[-1,1]](Ndata)]: # uniformly distributed numbers XD:=Matrix(1..Ndata, 1..2, datatype=float[8], order=C_order): # will hold the data abscissae YD:=Vector(1..Ndata,datatype=float[8], order=C_order): # will hold the function values X:=Vector(1..2,datatype=float[8]): for i from 1 to Ndata do for j from 1 to dim do    X[j]:=Xr[i*j]:  od: YD[i]:= g(X[1],X[2])+ 0.01*Noise[i]:        # record function values for j from 1 to dim do   XD[i,j]:=X[j]:             # record data abscisae od: od:

 

>	pointplot3d({seq([XD[k,1],XD[k,2],YD[k]],k=1..Ndata)},axes=BOXED,orientation=[-140,60]);

 

 

Interpolation and approximation 

First calculate the Lipschitz constant, then plot the interpolant with no monotonicity constraints 

>	LipIntComputeLipschitz(dim, Ndata, XD,YD); LC := LipIntGetLipConst();

 

 

>

f_aux:=(t1,t2)->LipIntValue(dim,Ndata, convert([t1,t2],Vector,datatype=float[8]), XD,YD,LC): R:=20.0: #plot resolution sdata :=[seq([ seq([i/R,j/R,f_aux(i/R,j/R)], i=R*Lo[1]..R*Up[1])], j=R*Lo[2]..R*Up[2])]: S1:=surfdata(sdata,axes=normal,grid=[100,100],style=HIDDEN, axes=BOXED,orientation=[-140,60]): S2:=pointplot3d({seq([XD[k,1],XD[k,2],YD[k]],k=1..Ndata)},axes=BOXED,color=BLACK): display({S1,S2});

 

 

The interpolant above is not necessarily monotone, but we can impose monotonicity as shown below (using LipIntValueCons command) 

>

f_aux:=(t1,t2)->LipIntValueCons(dim,Ndata, Cons, convert([t1,t2],Vector,datatype=float[8]), XD,YD,LC): sdata :=[seq([ seq([i/R,j/R,f_aux(i/R,j/R)], i=R*Lo[1]..R*Up[1])], j=R*Lo[2]..R*Up[2])]: S1:=surfdata(sdata,axes=normal,grid=[100,100],style=HIDDEN, axes=BOXED,orientation=[-140,60]): S2:=pointplot3d({seq([XD[k,1],XD[k,2],YD[k]],k=1..Ndata)},axes=BOXED,color=BLACK): display({S1,S2});

 

 

Next we impose monotonicity only in part of the domain, namely for x1,x2 <=0.5 (left region, shaded) 

>	Region:=Vector(1..2,datatype=float[8],[0.5,0.5]): Cons[1]:=+1; Cons[2]:=-1;       # try to change these to see what happens LipIntComputeLocalLipschitz(dim,Ndata,XD,YD); f_aux1:=(t1,t2)->LipIntValueLocalConsLeftRegion(dim,Ndata, Cons, convert([t1,t2],Vector,datatype=float[8]), XD,YD,Region):

 

 

 

>

sdata :=[seq([ seq([i/R,j/R,f_aux1(i/R,j/R)], i=R*Lo[1]..R*Up[1])], j=R*Lo[2]..R*Up[2])]: S1:=surfdata(sdata,axes=normal,grid=[100,100],style=HIDDEN, axes=BOXED,orientation=[-140,60]): S2:=pointplot3d({seq([XD[k,1],XD[k,2],YD[k]],k=1..Ndata)},axes=BOXED,color=BLACK): S3:=plot3d(0,x=0..Region[1],y=0..Region[2], color=BLACK): display({S1,S2,S3});

 

 

Section 3: Reading data from a file and its interpolation 

The following commands just read some data from a text file 

>

fp := fopen(`aluminiumdata.txt`,READ,TEXT): dimensions := fscanf(fp, `%d %d`): dim:= dimensions[1]: Ndata:=dimensions[2]: XD:=Matrix(1..Ndata, 1..dim,datatype=float[8], order=C_order): YD:=Vector(1..Ndata,datatype=float[8], order=C_order): # will hold the function values X:=Vector(1..dim,datatype=float[8]): for i from 1 to Ndata do  for j from 1 to dim do dat := fscanf(fp, `%f `): XD[i,j]:= dat[1]:  od:  dat := fscanf(fp, `%f`):  YD[i]:= dat[1]: od: fclose(fp): Lo:=Vector(1..2,datatype=float[8],[-2.3,-0.1]):   Up:=Vector(1..2,datatype=float[8],[0,1.15]):

 

The next command calculates the local Lipschitz constants 

>	LipIntComputeLocalLipschitz(dim, Ndata, XD,YD);

 

>

f_aux:=(t1,t2)->LipIntValueLocal(dim,Ndata, convert([t1,t2],Vector,datatype=float[8]), XD,YD): R:=30.0: #plot resolution sdata :=[seq([ seq([i/R,j/R,f_aux(i/R,j/R)], i=R*Lo[1]..R*Up[1])], j=R*Lo[2]..R*Up[2])]: S1:=surfdata(sdata,axes=normal,grid=[100,100],style=HIDDEN, axes=BOXED,orientation=[-140,60]): S2:=pointplot3d({seq([XD[k,1],XD[k,2],YD[k]],k=1..Ndata)},axes=BOXED,color=BLACK): display({S1,S2});

 

 

It looks like the data are monotone, and the interpolant is not. Next we impose monotonicity constraints, and the result  

>

Cons:=Vector(1..2,datatype=integer[4],[1,1]): f_aux1:=(t1,t2)->LipIntValueLocalCons(dim,Ndata, Cons, convert([t1,t2],Vector,datatype=float[8]), XD,YD): R:=20.0: #plot resolution sdata :=[seq([ seq([i/R,j/R,f_aux1(i/R,j/R)], i=R*Lo[1]..R*Up[1])], j=R*Lo[2]..R*Up[2])]: S1:=surfdata(sdata,axes=normal,grid=[100,100],style=HIDDEN, axes=BOXED,orientation=[-140,60]): S2:=pointplot3d({seq([XD[k,1],XD[k,2],YD[k]],k=1..Ndata)},axes=BOXED,color=BLACK): display({S1,S2});

 

 

Section 4: Another example, now using bounds on function values 

>	g := proc( x, y)  # user's test function        min(0.8,x*y/(2-x-y+x*y));    end proc;

 

>	plot3d(g(x,y),x=0..1,y=0..1, axes=boxed,orientation=[-80,70],grid=[50,50]);

 

 

 

Next we define the upper and lower bounds on function values. Note that the last (dim+1)-st argument is the parameter containing the Lipschitz constant 

>	low := proc( x, y,p)   max(0,1-p*((1-x)+(1-y)));  end proc; up:=proc(x,y,p) min(x,y); end proc;

 

>	PL1:=plot3d({low(x,y,1)},x=0..1,y=0..1, axes=boxed,orientation=[-60,80],grid=[40,40],style=WIREFRAME,color=BLUE): PL2:=plot3d({up(x,y,1)},x=0..1,y=0..1, axes=boxed,orientation=[-60,80],grid=[20,20],style=WIREFRAME,color=RED): display({PL1,PL2});

 

 

 

 

In the above plot, the graph of the function must lie between the red and blue surfaces. Next generate the data with noise 

>

Lo:=Vector(1..2,datatype=float[8],[0,0]):   Up:=Vector(1..2,datatype=float[8],[1,1]): Cons:=Vector(1..2,datatype=integer[4],[1,1]): dim:=2: Ndata:=100: Xr:=[stats[random, uniform]((Ndata)*dim)]: # uniformly distributed numbers Noise:=[stats[random, uniform[-1,1]](Ndata)]: # uniformly distributed numbers XD:=Matrix(1..Ndata, 1..2, datatype=float[8], order=C_order): # will hold the data abscissae YD:=Vector(1..Ndata,datatype=float[8], order=C_order): # will hold the function values TD:=Vector(1..Ndata,datatype=float[8], order=C_order): # will hold the smoothened values X:=Vector(1..2,datatype=float[8]): for i from 1 to Ndata do for j from 1 to dim do    X[j]:=Xr[i*j]:  od: YD[i]:= g(X[1],X[2])+ 0.1*Noise[i]:        # record function values for j from 1 to dim do   XD[i,j]:=X[j]:             # record data abscisae od: od:

 

Now we smooth the data, and ensure they also satisfy low(x)<y<up(x) 

>	LC:=1.0: LipIntSetBounds(up,low); LipIntSmoothLipschitz(dim, Ndata,XD,YD,TD,LC,0,0,0); f_aux2:=(t1,t2)->LipIntValue(dim,Ndata, convert([t1,t2],Vector,datatype=float[8]), XD,TD,LC):

 

Plot the surface and the original (red croses) and smoothened data (black dots). It lies between the upper and lower bounds as required 

>

R:=20.0: #plot resolution sdata :=[seq([ seq([i/R,j/R,f_aux2(i/R,j/R)], i=R*Lo[1]..R*Up[1])], j=R*Lo[2]..R*Up[2])]: S1:=surfdata(sdata,axes=normal,grid=[100,100],style=HIDDEN, axes=BOXED,orientation=[-10,65]): S2:=pointplot3d({seq([XD[k,1],XD[k,2],TD[k]],k=1..Ndata)},axes=BOXED,color=BLACK): S3:=pointplot3d({seq([XD[k,1],XD[k,2],YD[k]],k=1..Ndata)},axes=BOXED,color=RED,symbol=CROSS): display({S1,S2,S3});

 

 

The bounds also work when estimating the Lipschitz constant from noisy data using sampl splitting or cross-validation 

>	LipIntComputeLipschitzSplit(dim,Ndata,XD,YD,TD,0.4,0); `Estimated Lipschitz constant from noisy data`;LC := LipIntGetLipConst();;

 

 

 

To clear the bounds use command below. You can plot the graph of the interpolant again to see that it does not satisfy the bounds anymore 

>	LipIntClearBounds();

 

>

 

References 

The methods implemented in LibLip are based on the results outlined in: 

 

Beliakov, G., A review of applications of the Cutting Angle methods, in Continuous Optimization, A. Rubinov and V. Jeyakumar, Editors. 2005, Springer: New York. p. 209-248.
 

Beliakov, G., Interpolation of Lipschitz functions. J. of Comp. and Applied Mathematics, 2006. 196: p. 20-44.

Beliakov, G., Monotonicity preserving approximation of multivariate scattered data. BIT, 2005. 45: p. 653-677. 

 

Please cite these papers when you need to reference LibLip 

Please consult http://www.deakin.edu.au/~gleb for updates. 

We appreciate your comments about your experiences using LibLip and this Maple toolbox. Please send your feedback to 

gbeliako@gmail.com 

This project was carried out with the assistance of Ilya Khriapin, which is greatefully acknowledged. 

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.