Commit 79b51595 authored by gonciarz's avatar gonciarz
Browse files

Latest benchmark in R added

parent 09cc1583
%%
dimension = 2;
decomposition = 2;
%% Input parameters
m = dimension;
decomp = decomposition;
scale = 1/decomp;
cube = decomp^m;
disp("numOfCubes=" + cube);
range = [-1; 1];
mask = 2^dimension - 1;
for
\ No newline at end of file
## To install 'chebpol' uncomment following line
# install.packages("chebpol")
library(chebpol)
library(R.matlab)
# clear all data in workspace
rm(list = ls())
## ---------- read points saved in matlab
# d <- readMat('/Users/gonciarz/ws/repo/polyapprox/Newton Interpol/finalResults/benchInterpolation_deg_2-1-150_dim_2_lp_2_decomp_2_rndPts_100_testPoints.mat')
# d <- readMat('/Users/gonciarz/ws/repo/polyapprox/Newton Interpol/finalResults/benchInterpolation_deg_2-1-150_dim_3_lp_2_decomp_2_rndPts_100_testPoints.mat')
# d <- readMat('/Users/gonciarz/ws/repo/polyapprox/Newton Interpol/finalResults/benchInterpolation_deg_2-1-50_dim_4_lp_2_decomp_2_rndPts_100_testPoints.mat')
d <- readMat('/Users/gonciarz/ws/repo/polyapprox/Newton Interpol/finalResultsRunge1decomp1/benchInterpolation_deg_2-1-150_dim_3_lp_2_decomp_1_rndPts_100_testPoints.mat')
summary(d)
numOfDataPts=d$numOfDataPts;
degreeRange=d$degreeRange;
randPts=d$randPts;
# degreeRange = degreeRange[1:10] # limit range of degrees if needed
## tested function - Runge
f <- function(x) 1/(1+1*sum(x^2))
benchFunction <- function(f, fi, currPts) {
maxError = 0;
for(col in 1:ncol(currPts)) {
x = currPts[,col]
expectedVal = f(x)
interpolVal = fi(x)
ae = abs(expectedVal - interpolVal)
if (ae > maxError) {maxError = ae}
}
return(maxError);
}
## go through all
errorsPerDegree = c()
for (degree in degreeRange) {
currPts = randPts[[degree]][[1]];
m=dim(currPts)[1]; ## first value is a number of elements per data point = dimension
currNumOfDataPts = numOfDataPts[[degree]][[1]];
N=numOfDataPts[[degree]][[1]];
NinOneDir = ceiling(N^(1/m))
Nall = NinOneDir^m
## ------------------- prepare grid/knots needed later to generate interpolant
gridOneDim = seq(-1, 1, length.out=NinOneDir);
grid=list()
dims=c()
i = 1;
while(i <= m) {
grid[[i]] <- gridOneDim
dims=c(dims, NinOneDir + 0)
i = i+1
}
knots=t(expand.grid(grid))
cat(sprintf("Degree=%d Number of points: PIP=%d, rounded in one dim=%d in all dims=%d\n", degree, N, NinOneDir, Nall))
## ----------- generate interpolant ------------------------
# method='multilin'; fm <- ipol(f, grid=grid, method=method)
# method='fh'; fh <- ipol(f, grid=grid, method=method);
# method='hstalker'; fi <- ipol(f, grid=grid, method=method)
# method='chebyshev'; fc <- ipol(f, dims=dims, method=method)
# method='uniform'; fi <- ipol(f, dims=dims, method=method)
# method='polyharmonic'; fi <- ipol(f, knots=knots, k=2, method=method);
methodNames=c()
method='hstalker'; fhs <- ipol(f, grid=grid, method=method)
methodNames=c(methodNames, method)
maxError <- benchFunction(f, fhs, currPts);
errorsPerDegree = c(errorsPerDegree, maxError)
method='uniform'; fu <- ipol(f, dims=dims, method=method)
methodNames=c(methodNames, method)
maxError <- benchFunction(f, fu, currPts);
errorsPerDegree = c(errorsPerDegree, maxError)
method='multilin'; fm <- ipol(f, grid=grid, method=method)
methodNames=c(methodNames, method)
maxError <- benchFunction(f, fm, currPts);
errorsPerDegree = c(errorsPerDegree, maxError)
method='fh'; fh <- ipol(f, grid=grid, method=method);
methodNames=c(methodNames, 'Floater-Hormann')
maxError <- benchFunction(f, fh, currPts);
errorsPerDegree = c(errorsPerDegree, maxError)
method='chebyshev'; fc <- ipol(f, dims=dims, method=method)
methodNames=c(methodNames, 'Chebyshev')
maxError <- benchFunction(f, fc, currPts);
errorsPerDegree = c(errorsPerDegree, maxError)
# knots <- matrix(runif(m*N), m)
method='simplexlinear'; fsl <- ipol(f, knots=knots, method=method)
methodNames=c(methodNames, method)
maxError <- benchFunction(f, fsl, currPts);
errorsPerDegree = c(errorsPerDegree, maxError);
}
numPlots=6;
errorsPerDegree = matrix(errorsPerDegree, nrow = numPlots)
for (i in 1:numPlots) {
if (i==1) plot(degreeRange, errorsPerDegree[i,], ylab='error', xlab='degree', log='y', type='b', col='red', xlim=range(degreeRange), ylim=range(errorsPerDegree))
else lines(degreeRange, errorsPerDegree[i,], type='b', col='red')
}
d <- writeMat('/Users/gonciarz/ws/repo/polyapprox/Newton Interpol/finalResultsRunge1decomp1/benchInterpolation_deg_2-1-150_dim_3_lp_2_decomp_1_rndPts_100_R.mat', RDATA=errorsPerDegree, degreeRange=degreeRange, methodNames=methodNames)
......@@ -8,7 +8,7 @@ degreeRange = degreeStart:degreeStep:degreeStop;
m = dimension;
p = lpDegree;
S = numOfRandom;
B =round(numOfRandom * 0.3); % num of normalize to reach boundary
B =round(numOfRandom * 0.0); % num of normalize to reach boundary
decomp = decomposition;
scale = 1/decomp;
cube = decomp^m;
......@@ -313,8 +313,7 @@ end
Kn=degreeRange;
figure(1)
figure %(3)
clf;
hold on
plot(Kn, DATA(1,:), '-.ro');
......@@ -324,8 +323,7 @@ plot(Kn, DATA(9,:), '-.co');
legend({'PIP', 'Chebyshev', 'cubic spline', '5th order spline'});
set(gca, 'YScale', 'log')
figure(2);
figure%(4);
clf;
hold on;
plot(DATA(2,:),DATA(1,:), '-.ro');
......@@ -334,4 +332,37 @@ plot(DATA(4,:),DATA(6,:), '-.go');
plot(DATA(4,:),DATA(9,:), '-.co');
legend({'PIP', 'Chebyshev', 'cubic spline', '5th order spline'});
set(gca, 'YScale', 'log')
% ----------------- plots for chebfun1 R vs Matlab
Kn=degreeRange;
figure;
clf;
hold on
plot(Kn, DATA(1,:), '-.ro');
plot(Kn, RDATA, '-.bo');
legend({'PIP Matlab', 'Chebyshev R'});
set(gca, 'YScale', 'log')
figure
clf;
Kn=degreeRange';
hold on
plot(Kn, DATA(1,:)', '-.ro');
plot(repmat(Kn', [size(RDATA, 1) 1])', RDATA', '-.o');
methodNames=['PIP'; methodNames]
legend(methodNames);
set(gca, 'YScale', 'log')
end
degreeStart = 2;
degreeStep = 1;
degreeStop = 38;
dimension = 5;
lpDegree = 2;
decomp = 1;
numOfRandom = 100;
outputFileString = "_benchInterpolation_deg_" + degreeStart+"-"+degreeStep+"-"+degreeStop+"_dim_"+dimension+"_lp_"+lpDegree+"_decomp_"+decomp+"_rndPts_"+numOfRandom;
% Data in DATA starts with degreeStart for index 1
degreeRange = degreeStart:degreeStep:degreeStop;
DATA=DATA(:,1:length(degreeRange));
save(outputFileString, 'DATA', 'degreeRange');
% data in randPts* is always from degree 1 so if degreeStart is > 1 then
% first entries are empty
degreeRange = degreeStart:degreeStep:degreeStop;
randPts=randPts(1:degreeStop);
randPtsNotShifted=randPtsNotShifted(1:degreeStop);
numOfDataPts=numOfDataPts(1:degreeStop);
save(outputFileString + '_testPoints', 'randPts', 'randPtsNotShifted', 'degreeRange', 'numOfDataPts');
......@@ -2,7 +2,7 @@ function result=interpolatedFunction(varargin)
% Currently implemented: Runge
s = 0;
for k = 1:nargin
s = s + (varargin{k}).^2;
s = s + (0.0 + varargin{k}).^2;
end
result = 1 ./ (1 + 25 * s);
end
x = linspace(-1, 1, 100);
y = interpolatedFunction(x);
S=3;
B=1;
decomp = 4;
m = 2;
cube = decomp^m;
hold on;
% figure;
% clf;
plot(x, y);
% random points per 'cube'
Y=cell(cube, 1);
for c = 1:cube
X = 2*rand(m,S)-1;
for i = 1:B
X(:,i) = X(:,i) / norm(X(:,i), Inf);
end
Y{c} = X;
end
% all random points in one array
XA = cat(2, Y{:});
%--------------------------- complexity PIP -------------------------------
x = 1:1:400;
figure(1);
clf;
hold on;
%dim5
y=[0.068595,0.080618,0.164914,0.429148,0.861793,1.109939,1.955869,3.120737,4.610426,7.375446,9.024691,13.720329,18.581056,26.546650,35.485266,45.240309,57.670418,72.810503,94.934136,113.919853,149.565036,177.916349,218.346931,262.904574,309.189193,375.669441,444.467968,540.464604,625.657736,730.426194,831.922605,951.727794,1124.472836,1275.465440,1491.960815,1684.027319,1973.791268,2189.458625,2444.460587,2752.667929,3097.610491,3487.305102,3889.743916,4383.066208];
%dim3
y=[0.042923,0.010981,0.038157,0.039698,0.064943,0.065572,0.069849,0.093965,0.166853,0.159949,0.229381,0.284271,0.352861,0.361768,0.444579,0.536229,0.489327,0.505584,0.581555,0.664864,0.755573,0.871308,0.961611,1.175756,1.609597,1.595194,1.880962,1.762606,1.799570,2.232791,2.403479,2.501581,2.735326,2.764621,3.396517,3.252242,4.274000,4.163604,4.169254,4.522867,4.988126,6.031072,5.662595,6.845605,6.484866,8.257599,7.291528,8.332987,8.924471,9.921514,9.181155,10.498393,10.857187,10.790682,12.188421,12.118465,13.290955,14.366728,13.771533,14.656861,16.310239,16.616867,17.642616,17.698830,19.253415,21.074192,22.251436,22.296505,23.498756,24.315911,25.196793,26.940557,27.086491,28.192804,29.357527,30.584914,31.764962,33.101138,36.809997,35.670681,37.234316,38.847478,41.539658,42.307584,42.937759,45.948805,47.558506,46.885071,48.939754,50.556941,51.952392,53.387160,55.414307,57.348534,59.710202,61.544933,62.842772,65.222995,70.733295,71.299272,73.135282,75.718759,78.980019,81.578860,82.581590,87.194608,89.014305,89.450071,91.557574,94.560760,98.174638,100.289032,102.324618,103.502937,108.524388,111.914311,115.808172,118.977738,121.040690,125.326070,126.686040,127.828721,131.107967,137.804792,141.288137,142.678981,148.191089,150.995706,157.531562,161.301085,163.976770,168.405755,174.477934,176.183050,182.581646,180.780023,191.658143,194.585915,197.123733,199.172245,202.587527,208.993202,218.611585,224.391471,225.131756,239.026079,241.510131,247.398488,252.294571];
%dim4
y = [0.036567,0.035758,0.070747,0.090545,0.204911,0.251867,0.385055,0.546604,0.752298,1.004189,1.189549,1.576459,2.047989,2.638625,3.280919,4.072443,5.018704,6.115649,7.394873,8.868044,10.522573,12.338405,14.544691,16.845152,19.646209,22.996866,26.433046,30.537880,34.526292,39.962571,46.314780,49.902327,56.368123,62.187966,68.871624,76.984102,85.035701,95.838810,105.561667,117.162753,126.931202,142.664241,153.367717,167.070117,180.846971,196.602749,213.630991,233.630647,254.265501,271.558199,297.292263,321.382869,343.816157,371.637104,402.888305,432.030409,460.237491,492.678402,527.265299,564.172925,605.831012,649.773613,690.193999,733.061375,782.256328,835.960005,885.642042,942.796434,997.808570,1049.808735,1124.219801,1184.622175,1254.549762,1341.800759,1416.828128,1489.658453,1574.207647,1666.052247,1760.246439];
%dim4 prep
y = [0.023330,0.011125,0.004075,0.006301,0.019437,0.016700,0.024261,0.038164,0.053137,0.074170,0.088373,0.119216,0.187748,0.247137,0.289953,0.350036,0.447462,0.562917,0.719119,0.875249,1.090068,1.325238,1.621264,1.991195,2.379992,2.865156,3.460546,4.128418,4.981556,5.992856,7.205303,8.937905,10.179605,12.050135,14.290907,17.005148,19.968695,23.588693,27.700260,32.250183,37.574921,43.384922,50.646408,58.680900,66.629561,76.660273,87.457001,101.760076,114.329306,129.183745,148.387160,166.760960,191.576186,210.641163,238.893234,273.041613,303.106786,341.666070,380.895036,423.077304,472.675391,523.977158,573.333400,648.064024,714.259692,789.020118,872.684348,966.418254,1049.104581,1161.230738,1264.581963,1378.165200,1516.151536,1673.737331,1822.768997,1978.143036,2163.801008,2351.156403,2563.257184];
y = 1 .* exp(x);
plot(x, y);
y = 0.1 .* exp(x);
plot(x, y);
y = (x.^14) .* exp(x);
plot(x, y);
y = x.^5;
plot(x, y);
y = 1 .* exp(2.*x);
plot(x, y);
y = 1 .* exp(0.5.*x);
plot(x, y);
set(gca, 'YScale', 'log')
y = y';
x = 2:(numel(y) + 1); x=x';
set(gca, 'XScale', 'log')
g = fittype('b * x^c');
f0 = fit(x, y, g);%, 'StartPoint', [2, 6, 2]);%, 'StartPoint', [ [ones(size(x)), -exp(-x) ] \ y; 1]);
xx = linspace(1,numel(y)+2,50);
figure(4);
plot(x,y,'o-',xx,f0(xx),'r-');
% ------------------------------------- chebyshev diff MATLAB vs R -------
figure(1)
load('R.mat')
clf;
CFR1= chebfun(@interpolatedFunction, degree, 'chebkind', 1);
CFR2= chebfun(@interpolatedFunction, degree, 'chebkind', 2);
x = linspace(-1, 1, 100);
y1 = CFR1(x);
y2 = CFR2(x);
plot(x, yf, 'c')
hold on
plot(x, yi, 'b-.', 'LineWidth', 2)
plot(x, y1, 'r-.', 'LineWidth', 1.5)
plot(x, y2, 'g', 'LineWidth', 2)
classdef A
properties
end
% set(gca, 'YScale', 'log')
end
pts1 = chebpts(degree, [-1,1], 1)
pts2 = chebpts(degree, [-1,1], 2)
plot(pts1, CFR1(pts1), 'mo', 'LineWidth', 3);
plot(pts2, CFR2(pts2), 'k*', 'LineWidth', 3);
classdef B
legend({'Runge', 'R interp', 'Matlab kind 1', 'Matlab kind 2', 'chebpts kind 1', 'chebpts kind 2'});
end
args <- commandArgs(trailingOnly = TRUE)
degreeStart = as.numeric(args[1])
degreeStep = as.numeric(args[2])
degreeStop = as.numeric(args[3])
dimension = as.numeric(args[4])
lpDegree = as.numeric(args[5])
decomposition = as.numeric(args[6])
numOfRandom = as.numeric(args[7])
baseName = sprintf("benchInterpolation_deg_%d-%d-%d_dim_%d_lp_%d_decomp_%d_rndPts_%d", degreeStart, degreeStep, degreeStop, dimension, lpDegree, decomposition, numOfRandom)
cat(baseName)
#benchInterpolation_deg_2-1-150_dim_2_lp_2_decomp_2_rndPts_100_testPoints.mat
#benchInterpolation_deg_2-1-150_dim_2_lp_2_decomp_2_rndPts_100
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment