MATLAB: Finding wilkinson notation parser in backend for fitlm Matlab R2015a

custom loss functiondesign matrixfitlmterms matrixwilkinson notation

Version
MATLAB Production Server/R2015a
Problem
I want to be able to create a regression design matrix from a wilkinson notation string as is done in the backend of fitlm.
From using the debugger it appears that the wilkinson notation string is expanded out into a full wilkinson notation string (time^2 -> time + time^2) in line 856 of MATLAB\MATLAB Production Server\R2015a\toolbox\stats\classreg\@LinearModel\LinearModel.m using LinearModel.createFormula.
After this it seems that in either line 858 or 863 of the same file the regression design matrix (with columns for each term in the regression) is created by calling either assignData or doFit. However I cannot open either of these files, and after grepping the installation directory for identically named files and putting break points in them it appears they're not there, or at least not in the form assignData.* or doFit.* anyway.
Could you tell me where I can find this backend code?
Example
If I'm modelling something where time ~ a*dist + b*dist^2, (a = b = 1 in the code below, plus a small amount of noise) the wilkinson notation would be 'time ~ dist^2'. When using fitlm it would generate a dist^2 column from the dist column. But the matrix with that generated factor is kept internal to the model. I would like to be able to get that matrix.
table = array2table([1,2,3,4;2.01,5.99,12.01,19.99]', 'VariableNames', {'dist','time'})
mdl = fitlm(table,'time ~ dist^2')
The matrix I'm looking for is the one below:
D = array2table([1,2,3,4;1,4,9,16]', 'VariableNames', {'time','time^2'})
This is the matrix that the linear model is actually fit on, so I could use fmincon or fminunc on it with my custom loss function to get the optimised parameters. Obviously with a simple relationship like 'time -> time^2' it's not difficult to write my own parser, but once I have a large complicated wilkinson notation formula with lots of factors which functions of multiple variables to various powers (eg 'x:(y^3):(time^2)') it get's more complicated and it seems like it would be easier to use the internal Matlab one.
Context
The reasoning for this is that I'd like to be able to fit a linear regression model with a different loss function than RMSE. Specifically I'd like to optimise for percentage mean absolute error.
pmae = @(x,y) mean(abs(x - y) ./ y)
I am open to other ways of doing this without wilkison notation parser, but it would be more useful to have the parser as I'll need it elsewhere as well.

Best Answer

After some digging, let's see if this works for you.
The Statistics and Machine Learning Toolbox has a function D = x2fx(X,model) that converts a predictor matrix into a design matrix using variables provided by the linear model object returned by mdl = fitlm(X,y).
Using your demo (thanks for providing it!), 
table = array2table([1,2,3,4;2.01,5.99,12.01,19.99]', 'VariableNames', {'dist','time'})
mdl = fitlm(table,'time ~ dist^2')
pred = [(1:numel(mdl.predict))', mdl.predict]; 
D = x2fx(pred, mdl.Formula.Terms)
% D =
%      1     1     1
%      1     2     4
%      1     3     9
%      1     4    16
mdl.Formula     % = time ~ 1 + dist + dist^2
Related Question