close all
clearvars
uqlab

%2.2 RANDOM FIELDS: 4 random field-material properties and corrosion growth
% Random field for material properties(UTS & SMYS) is modelled as a Gaussian or lognormal random process
% A) SMYS
%Specify the type of input
RFInput1.Type = 'RandomField';
%Specify the type of RF
RFInput1.RFType = 'Gaussian' ;
%Specify the correlation family
RFInput1.Corr.Family = 'Gaussian' ;
%Specify the discretization domain (Duration of 30 years, expressed in months 360)
RFInput1.Domain= [0, 360]';
%Specify the number of points for building the mesh
RFInput1.SPD = 361 ;
%Specify the correlation length ( 5years = 60 months)
RFInput1.Corr.Length = 60;
%Specify the expansion order
RFInput1.ExpOrder = 7;
%Specify the mean and standard deviation of the random field
RFInput1.Mean = 624.48 ;
RFInput1.Std = 6.24 ;
%Create an INPUT object based on the specified random field options
myRFInput1 = uq_createInput(RFInput1);
%Print a report of the created INPUT object
%uq_print(myRFInput1);
%uq_display(myRFInput1);

% Lognormal for SMYS
%RFInput1.Type = 'RandomField';
%Specify the type of RF
%RFInput1.RFType = 'Lognormal' ;
%Specify Karhunen-Loève as discretization scheme
%RFInput1.DiscScheme = 'KL' ;
%Specify the correlation family
%RFInput1.Corr.Family = 'Gaussian' ;
%Specify the correlation length
%RFInput1.Corr.Length = 60 ;
%Specify the mesh when the random field is sampled
%RFInput1.Mesh = linspace(0,360,361)' ;
%Specify the mean and standard deviation of the random field
%RFInput1.Mean = 6.436 ;
%RFInput1.Std = 0.053 ;
%Create an INPUT object based on the specified random field options
%myRFInput1 = uq_createInput(RFInput1);
%Print a report of the created INPUT object
%uq_print(myRFInput1);
%uq_display(myRFInput1);

% B) UTS
%Specify the type of input
RFInput2.Type = 'RandomField';
%Specify the type of RF
RFInput2.RFType = 'Gaussian' ;
%Specify the correlation family
RFInput2.Corr.Family = 'Gaussian' ;
%Specify the discretization domain (Duration of 30 years, expressed in 360 months)
RFInput2.Domain= [0, 360]';
%Specify the number of points for building the mesh
RFInput2.SPD = 361 ;
%Specify the correlation length (5 YEARS = 60 months)
RFInput2.Corr.Length = 60;
%Specify the expansion order
RFInput2.ExpOrder = 7 ;
%Specify the mean and standard deviation of the random field
RFInput2.Mean = 732.00 ;
RFInput2.Std = 7.32;
%Create an INPUT object based on the specified random field options
myRFInput2 = uq_createInput(RFInput2);
%Print a report of the created INPUT object
%uq_print(myRFInput2);
%uq_display(myRFInput2);

% Lognormal for UTS
%RFInput2.Type = 'RandomField';
%Specify the type of RF
%RFInput2.RFType = 'Lognormal' ;
%Specify Karhunen-Loève as discretization scheme
%RFInput2.DiscScheme = 'KL' ;
%Specify the correlation family
%RFInput2.Corr.Family = 'Gaussian' ;
%Specify the correlation length
%RFInput2.Corr.Length = 60 ;
%Specify the mesh when the random field is sampled
%RFInput2.Mesh = linspace(0,360,361)' ;
%Specify the mean and standard deviation of the random field
%RFInput2.Mean = 6.595;
%RFInput2.Std = 0.039 ;
%Create an INPUT object based on the specified random field options
%myRFInput2 = uq_createInput(RFInput2);
%Print a report of the created INPUT object
%uq_print(myRFInput2);
%uq_display(myRFInput2);

% Random field for corrosion growth in the defect depth aand length is modelled as stocahstic process
% C) CGR in defect depth
%BD2 dc
%Specify the type of input
RFInput3.Type = 'RandomField' ;
%Specify the type of RF
RFInput3.RFType = 'Gamma' ;
%Specify the type of correlation
RFInput3.Corr.Type = 'Separable' ;
%Specify the correlation family
RFInput3.Corr.Family = 'Gaussian' ;
%Specify the correlation length
RFInput3.Corr.Length = 60 ;
%Specify the expansion order
RFInput3.ExpOrder = 10 ;
%Specify the mesh for the sampling of the random field
RFInput3.Mesh = linspace(0,360,361)';
%Specify the mean and standard deviation of the random field
RFInput3.Mean = 4.9464 ;
RFInput3.Std = 0.1063;
%Specify the conditional data
RFInput3.RFData.X = [0; 60; 120; 180; 240; 300; 360] ;
RFInput3.RFData.Y = [4.8; 4.85; 4.90; 4.95; 4.99; 5.04; 5.10] ;
%Create an INPUT object based on the specified random field options
myRFInput3 = uq_createInput(RFInput3);
%Print a report of the created INPUT object
%uq_print(myRFInput3);
%uq_display(myRFInput3);

% D) CGR in defect length
% Lc
%Specify the type of input
RFInput4.Type = 'RandomField' ;
%Specify the type of RF
RFInput4.RFType = 'Gamma' ;
%Specify the type of correlation
RFInput4.Corr.Type = 'Separable' ;
%Specify the correlation family
RFInput4.Corr.Family = 'Gaussian' ;
%Specify the correlation length
RFInput4.Corr.Length = 60 ;
%Specify the expansion order
RFInput4.ExpOrder = 10 ;
%Specify the mesh for the sampling of the random field
RFInput4.Mesh = linspace(0,360,361)';
%Specify the mean and standard deviation of the random field
RFInput4.Mean =210.92 ;
RFInput4.Std = 0.664;
%Specify the conditional data
RFInput4.RFData.X = [0; 60; 120; 180; 240; 300; 360] ;
RFInput4.RFData.Y = [210.00; 210.30; 210.60; 210.91; 211.22; 211.53; 211.85] ;
%Create an INPUT object based on the specified random field options
myRFInput4 = uq_createInput(RFInput4);
%Print a report of the created INPUT object
%uq_print(myRFInput4);
%uq_display(myRFInput4);

%3 - COMPUTATIONAL MODEL
%Insert developed PCE surrogate model .mat file
% Specify X, Xval, Y and Yval
X =  [624.4800  732.0000    4.7900  458.7600    7.9800  210.0000   30.0000   21.0800
  698.7191  649.0031    6.2734  459.3833    8.1948  240.9830   25.4821   24.2885
  524.8011  810.7650    6.7539  459.5569    8.1941  177.8882   26.0114   18.9514
  702.2099  805.1040    7.1454  459.6793    7.9401  188.4695   25.5121   23.7205
  725.4150  722.3977    3.0112  458.1524    8.1910  234.9569   34.5662   24.3022
  523.0512  653.6285    6.2418  459.4885    8.1267  179.2886   25.5272   24.2865
  541.9590  799.9709    2.5195  457.9773    7.7834  240.0637   25.9824   19.6003
  531.6219  652.9428    3.6635  457.9379    8.1763  239.1703   25.5745   23.7179
  548.0212  815.8323    3.6014  459.5384    8.1621  240.0106   25.8958   20.2834
  543.2446  658.6519    5.2072  457.9081    7.7767  182.8670   34.4758   18.5220
  533.0712  819.9537    3.7935  459.5414    8.1516  239.6965   34.3078   21.9106
  588.3707  805.2872    2.5699  457.9754    7.8369  191.6357   34.3259   18.4340
  590.3613  653.2711    6.9483  459.6597    7.8665  182.8457   25.4200   22.6034
  720.4838  806.2868    5.7877  457.8624    7.7811  218.7911   34.3357   23.6180
  648.5544  809.8824    7.1021  459.6370    7.7878  222.6042   27.8939   18.1628
  539.9638  764.9282    2.4764  457.8482    7.7716  181.0206   25.9262   20.6518
  539.5158  644.1088    6.0265  457.8727    8.1254  190.1524   31.9424   19.1781
  717.2670  668.7456    6.8759  459.3794    8.1426  193.9583   31.4895   24.2963
  532.0260  657.0039    7.1272  459.6785    7.7946  237.3212   33.7134   21.3380
  717.7133  814.4170    6.2213  458.0229    7.7697  232.6920   27.6972   19.2873
  629.3756  645.9560    6.2907  457.9982    7.8717  242.3155   33.3423   23.7655
  574.9957  649.3854    3.3928  459.6187    7.8140  241.4513   34.1215   23.0929
  528.7844  758.8774    6.2865  458.9144    7.8192  178.6843   25.3852   22.8561
  706.9700  646.8241    5.8021  459.3077    7.8680  240.7703   34.6223   18.1482
  608.4099  799.9836    6.2091  457.8755    8.1863  195.7223   33.3549   24.2863
  718.4943  795.5942    2.7893  458.1188    8.1533  233.8071   26.1774   18.1434
  529.8809  722.2985    4.4748  459.4717    7.7955  239.8618   25.7878   24.3134
  720.9837  801.1540    6.1360  457.8606    8.1640  185.5539   29.1658   17.8542
  526.1791  790.8692    7.0843  459.6036    8.1154  190.4614   33.4006   20.8267
  686.0528  794.1031    3.5176  459.5770    8.0850  236.1995   26.5392   18.1296
  593.1291  820.0314    3.9304  459.2087    7.7724  179.5751   27.3313   23.9348
  718.2144  816.6599    3.1226  458.3310    7.8567  185.3947   25.9275   24.1430
  528.5661  673.8380    6.5469  458.7158    7.7992  231.3066   32.4389   17.8377
  668.9718  762.6345    3.4068  459.6271    8.0133  178.0752   32.4409   18.4581
  722.3869  817.0348    6.9629  459.5679    7.9952  183.6812   33.8169   18.0835
  632.6075  646.4024    4.7424  458.2445    8.1858  181.2783   34.4373   17.9140
  638.0669  654.2135    3.3590  459.3990    7.7943  185.3341   34.2840   22.7055
  546.0827  741.4131    5.9887  459.6626    8.1561  179.8786   34.5311   24.1675
  704.7167  814.6539    6.2186  458.5611    8.0485  241.2956   33.9203   24.2525
  562.7706  657.4383    3.6351  459.1968    7.8359  183.5139   26.5660   23.8118
  543.7676  819.2638    3.7106  458.1378    8.1897  183.5688   28.2657   18.7068
  532.2066  819.1946    6.4671  458.3287    8.0837  222.0184   26.9854   18.1027
  537.0014  679.5736    3.3604  459.1841    8.0327  236.6204   34.5237   18.2371
  603.9031  659.5710    2.6235  457.9924    7.9266  240.8421   26.8233   22.8577
  661.6772  799.1272    6.1074  457.9524    7.8043  187.8514   33.9010   22.0832
  552.0477  672.8986    2.9688  458.0634    7.9129  181.1675   27.2005   19.1746
  545.5117  797.5641    6.6023  458.9409    7.7829  242.2753   27.5599   23.8258
  627.4937  662.6770    3.7246  458.0125    7.7849  241.3310   34.5519   18.4896
  709.8969  689.8087    4.6547  457.9094    7.8405  237.9741   27.2529   18.0405
  725.1462  753.9259    3.2335  459.1516    7.7835  240.9870   30.3685   23.6061
  523.3948  664.8445    3.7842  459.5416    7.7907  234.7535   28.7888   19.5948
  556.0656  676.5650    6.1448  457.8529    8.1257  177.8827   26.5347   20.7246
  522.9301  723.2564    6.1997  457.8934    7.8115  179.4860   29.0727   23.5873
  563.3262  707.3386    2.9000  457.8767    8.1444  237.4042   34.3572   19.8453
  701.1408  657.7691    6.6004  459.2957    7.8166  241.7922   25.7202   18.0041
  625.4683  676.1480    3.6594  459.6025    8.1408  182.6242   28.1042   24.0932
  523.3777  654.4564    5.9271  458.3457    8.1465  242.4074   31.9277   23.7591
  619.5433  818.7671    6.2006  457.9369    7.8889  238.8126   25.8782   23.4497
  683.8594  796.7294    6.1927  459.6530    8.1784  216.1805   25.4264   24.0267
  703.5035  647.4552    6.2238  458.4030    8.1400  204.6675   33.5913   18.0285
  725.6421  655.7271    5.2514  458.2393    8.1094  188.6630   27.1857   24.0647
  523.9004  655.4272    3.7844  459.3104    8.1770  202.0451   28.8740   18.9470
  532.3802  750.7303    5.1529  457.8540    8.1931  242.3025   34.4913   20.6519
  630.0947  733.8724    2.9945  458.9890    8.1848  182.5726   33.5555   24.1104
  540.9329  669.0313    5.0286  459.6752    8.1298  186.2642   28.7809   18.9156
  721.0326  655.5922    2.9157  458.3565    7.7902  180.7287   29.1263   20.2243
  626.3329  818.1844    5.8305  458.2816    8.1603  240.0699   32.8840   23.5367
  603.3794  686.3673    6.2881  459.5449    7.8511  186.6456   34.3974   17.8626
  540.0848  784.7746    2.9013  458.0797    7.8278  191.0617   26.6923   23.4786
  532.1533  652.1136    2.6401  458.1769    8.0720  200.6744   33.0732   22.1071
  536.6778  674.7118    6.9031  459.2547    8.1051  236.8954   26.7643   18.2343
  704.2120  665.6785    2.9564  458.6914    7.8441  236.8908   27.6985   24.1149
  669.6869  653.0226    3.7074  459.6469    7.8837  211.2295   25.5060   17.9770
  715.7288  661.4132    6.2710  457.9875    8.1953  179.3285   26.2212   21.6025
  564.4335  719.3610    5.7436  458.0627    7.8894  182.1670   26.0920   17.8325
  662.1797  816.5426    3.5508  459.2709    8.1603  235.4528   25.7049   19.4878
  551.5562  819.3402    3.0981  458.5610    8.0423  212.8235   34.3919   23.6987
  638.4686  679.2542    3.7626  459.6822    8.1963  240.2018   33.9796   20.6318
  538.8077  796.0909    6.2861  459.3180    8.1561  225.3874   34.2537   17.9960
  725.2775  805.7415    3.3451  459.5151    7.9351  222.9872   33.3272   20.7260
  710.3250  652.0959    6.2207  458.5957    7.8358  236.2527   25.7595   24.1594];

Y =  [747.9717
  754.7873
  752.2763
  754.4523
  755.6151
  749.1570
  755.6212
  754.8222
  762.5214
  751.3284
  754.5186
  758.8237
  755.5034
  749.1519
  755.1640
  761.5639
  758.3229
  753.5582
  747.9715
  751.5875
  759.9098
  750.2409
  757.9070
  750.2529
  784.6378
  749.2518
  751.9813
  760.6246
  754.7947
  760.1022
  759.5814
  756.7536
  756.6386
  754.4354
  757.5633
  754.8316
  761.6923
  753.3890
  749.4962
  758.6789
  752.0758
  752.1484
  748.9819
  750.3112
  754.7273
  757.5965
  749.7142
  754.5629
  759.5404
  754.5776
  756.4280
  760.0264
  754.6123
  750.0559
  784.9682
  754.7555
  749.7714
  756.0751
  747.9723
  760.0324
  755.1353
  758.9707
  750.1977
  754.7864
  754.6563
  786.3628
  754.8562
  754.7002
  750.2268
  754.9074
  762.8878
  761.1139
  759.8734
  749.7584
  749.8581
  756.7556
  750.6580
  754.6019
  747.9725
  749.2828
  751.6547];


Xval =  [725.0685  781.7816    6.9384  459.4912    7.7733  237.5589   33.1111   22.2832
  549.5551  774.2057    3.1280  458.0601    7.7784  239.2220   31.4810   22.3404
  662.4716  746.9819    6.4288  458.3088    8.1694  241.7961   27.5912   24.2198
  610.8708  787.0258    3.0931  458.2576    8.0310  241.4077   25.3918   21.0825
  543.0488  818.1634    6.8525  459.6527    8.1530  184.9194   27.5912   24.0002
  636.0680  815.6165    4.2852  458.4441    7.7875  184.6937   25.6384   18.4988
  714.7678  725.6351    4.8628  459.4133    7.7922  186.0720   34.5718   24.3214
  683.7380  687.5032    6.2537  457.9840    7.7775  179.9988   26.9656   23.8493
  687.1784  794.4096    2.8457  458.0223    8.1918  180.8782   32.2755   22.5621
  624.2650  818.5998    3.1922  459.2490    7.7646  181.7584   27.1295   18.7087
  718.3192  795.7153    3.3323  458.6267    8.1879  196.3894   26.0731   18.6682
  530.8525  668.2901    2.7510  458.0610    7.7689  216.1771   33.8673   22.6854
  590.3043  689.3204    6.7330  459.4536    7.9290  205.3067   32.7507   18.4200
  561.6153  767.1858    6.5974  459.5537    8.1307  232.3361   33.1744   23.4750
  686.9034  802.3570    6.1059  458.3425    7.8231  185.2575   31.1950   18.1270
  723.3648  660.6639    5.4002  458.2018    8.0062  194.3424   33.9041   24.0593];

Yval = [757.7108
  749.0831
  760.0174
  754.4514
  755.9609
  754.0085
  750.7095
  751.3386
  750.2988
  748.8666
  760.1596
  748.1521
  790.0111
  750.0867
  757.4803
  756.1016];


% Specify the marginals as follows:The probabilistic input model consists of eleven (8) independent random variables.
Input.Marginals(1).Name = 'SMYS';  % Tensile Yield Strength
Input.Marginals(1).Type = 'Gaussian';
Input.Marginals(1).Parameters = [624.48 32.97];  % (MPa)

Input.Marginals(2).Name = 'UTS';  % Ultimate Tensile Strength
Input.Marginals(2).Type = 'Gaussian';
Input.Marginals(2).Parameters = [732 28.52];  % (MPa)

Input.Marginals(3).Name = 'dc';  %  dmax or dcluster
Input.Marginals(3).Type = 'Gaussian';
Input.Marginals(3).Parameters = [4.790 1.544];  % (mm)

Input.Marginals(4).Name = 'D';  %  Outer Diameter
Input.Marginals(4).Type = 'Gaussian';
Input.Marginals(4).Parameters = [458.6 0.691];  % (mm)

Input.Marginals(5).Name = 't';  % Pipe wall thickness
Input.Marginals(5).Type = 'Gaussian';
Input.Marginals(5).Parameters = [7.98 0.071];  % (mm)

Input.Marginals(6).Name = 'Lc';  % Parameter 2 Effective corrosion cluster length
Input.Marginals(6).Type = 'Gaussian';
Input.Marginals(6).Parameters = [210 10.5];  % (mm)

Input.Marginals(7).Name = 'SL';  % Parameter 3 Interacting longitudinal distance
Input.Marginals(7).Type = 'Gaussian';
Input.Marginals(7).Parameters = [30 1.5];  % (mm)

Input.Marginals(8).Name = 'Sc';  % Parameter 4 Interacting circumferential distance
Input.Marginals(8).Type = 'Gaussian';
Input.Marginals(8).Parameters = [21.08 1.054];  % (mm)
%Create an INPUT object based on the specified marginals:
myInput = uq_createInput(Input);
%Specify an INPUT object: an INPUT object myInput can be specified with:
%MetaOpts.Input = myInput;
%  Select PCE as the metamodeling tool: Sparse regression methods (LARS, OMP, SP, BCS)
MetaOpts.Type = 'Metamodel';
MetaOpts.MetaType = 'PCE';
%MetaOpts.Display = 'quiet';
%MetaOpts.TruncOptions.qNorm = 0.75;
%MetaOpts.TruncOptions.MaxInteraction = 3;
% Specification of 
%Specify the range of the degrees to be compared by the adaptive algorithm:
MetaOpts.Degree = 10:30;
%MetaOpts.Degree = 14;
MetaOpts.Method = 'LARS';
MetaOpts.LARS.KeepIterations = 1;
%Use experimental design loaded from the data files:
MetaOpts.ExpDesign.X = X;
MetaOpts.ExpDesign.Y = Y;
%Set the maximum polynomial degree to 20:
%MetaOpts.Degree = 1:20;
%Provide the validation data set to get the validation error:
MetaOpts.ValidationSet.X = Xval;
MetaOpts.ValidationSet.Y = Yval;
%Create the metamodel object and add it to UQLab:
myPCE_LARSa = uq_createModel(MetaOpts);
%Specify a MODEL object: a MODEL object myModel can be specified with:
MetaOpts.FullModel = myPCE_LARSa;
%Print a summary of the resulting PCE metamodel:
%uq_print(myPCE_LARSa)
%5 - VALIDATION
%Evaluate the PCE metamodel at the validation set points:
%YPCE_LARS = uq_evalModel(myPCE_LARSa,Xval);
%Print the validation and leave-one-out (LOO) cross-validation errors:
%fprintf('PCE metamodel validation error: %5.4e\n', myPCE_LARSa.Error.Val)
%fprintf('PCE metamodel LOO error:        %5.4e\n', myPCE_LARSa.Error.LOO)
%save myPCE_LARSa

% Merge scalar inputs and random fields
myFullRFInput1 = uq_mergeInputs(myRFInput1.UnderlyingGaussian, myRFInput2.UnderlyingGaussian) ;
myFullRFInput2 = uq_mergeInputs(myRFInput3.UnderlyingGaussian, myRFInput4.UnderlyingGaussian) ;
myFullRFInput = uq_mergeInputs(myFullRFInput1, myFullRFInput2) ;

myFullInput = uq_mergeInputs(myInput, myFullRFInput) ;

% wrap the  model as a function handle 
ModelOpts.mHandle = @(X) uq_evalModel(myPCE_LARSa, X);
%ModelOpts.isVectorized = true;
ModelOpts.isVectorized = false;
ModelOpts.Parameters.SMYS = myRFInput1 ;
ModelOpts.Parameters.UTS = myRFInput2 ;
ModelOpts.Parameters.dc = myRFInput3 ;
ModelOpts.Parameters.Lc = myRFInput4 ;
%ModelOpts.Parameters(1) = myRFInput1 ;
%ModelOpts.Parameters(2) = myRFInput2 ;
%ModelOpts.Parameters(3) = myRFInput3 ;
%ModelOpts.Parameters(6) = myRFInput4 ;
myModel = uq_createModel(ModelOpts);
%MetaOpts.FullModel = myModel;



%4. EXPERMENTAL DESIGN OVER DESCRITISATION DOMAIN (360mobths)
%Construct a surrogate based on 150 model runs
X = uq_getSample(myFullInput, 360, 'LHS') ;
%Xm1 = uq_getSample(myFullInput, 360, 'MC') ;
%Xm = uq_getSample(myFullInput, 360, 'Sobol') ;
%Xm = uq_getSample(myFullInput, 360, 'Halton') ;
%Evaluate the computational model
fprintf('Evaluating the experimental design...')
Y = uq_evalModel(myModel, X) ;
fprintf(' done!\n\n')




