MATH920 Lesson03: DNN-ELS, Neural net classifier of elastic moduli

Data input

Read data from solution of elasticity equations: get data file names, extract grid info, find number of time slices

SetDirectory["/home/student/courses/MATH920/models/DNN-ELS/inhomogeneous/data"];
qFiles=FileNames["q*.hdf"];  tFiles=FileNames["t*.hdf"];
tInfo = Map[Import[#,{"HDF","Data"}][[2]]&,tFiles];
nt = Length[tFiles]

Extract time values at each slice

t = tInfo[[;;,3]]

Extract grid step size (assume dx=dy=h, use grid spacing along x)

h =(tInfo[[1,14]]-tInfo[[1,10]])/tInfo[[1,6]]

Read data file with material properties used in BEARCLAW simulation.

setprob=ImportString[Import["../setprob.data","Text"],"Table"]

Extract numerical values of interest.

{cp1,cs1,ρ1,cp2,cs2,ρ2,sig0,di,dcp,dcs,dρ} = Map[#[[1]]&,setprob]

Compute elastic moduli used for numerical simulation.

μ1=ρ1 (cs1^2); μ2=ρ2 (cs2^2); λ1=ρ1 (cp1^2)-2 μ1; ; λ2=ρ2 (cp2^2)-2 μ2;

Import the data files and extract number of time steps, components, grid sizes

q=Transpose[Table[Import[qFiles[[k]],{"HDF","Data"}],{k,nt}],{2,1,4,3}];
{nq,nt,m1,m2}=Dimensions[q]

Extract components of q into standard physics notation

{(Subscript[σ, 11]),(Subscript[σ, 22]),(Subscript[σ, 33]),(Subscript[σ, 12]),(Subscript[v, 1]),(Subscript[v, 2]),p}=q[[;;,;;,di+1;;m1-di,di+1;;m2-di]];
{nt,m1,m2}=Dimensions[p]

Define a subset of displacement values for visualization

xyuv=Table[Flatten[Table[{i1, i2,(Subscript[v, 1])[[k,i1,i2]], (Subscript[v, 2])[[k,i1,i2]]} ,{i1,1,m1,5},{i2,1,m2,5}],1],{k,nt}];

Visualize by plotting p=Tr(σ)/3 and displacement vectors

pdPlot[k_]:=Show[{ArrayPlot[Transpose[ p[[k]]],DataReversed->True ,PlotLegends->Automatic,ColorFunction->"Rainbow"],ListVectorPlot[xyuv[[k]]]}]

Investigate stress, displacement fields. Initial condition corresponds to an incoming plane wave from the left.

pdPlot[8]

pdPlot[20]

Process input data

(Subscript[X, 1])=Table[i1 h,{i1,m1},{i2,m2}];
(Subscript[X, 2])=Table[i2 h,{i1,m1},{i2,m2}];

Find displacements from velocity

zero=Table[0.,{i1,m1},{i2,m2}]; (Subscript[u, 1])={zero}; (Subscript[u, 2])={zero};
Monitor[For[k=2,k≤nt,k++,
         dt=(t[[k]]-t[[k-1]]);
         (Subscript[du, 1])= ((Subscript[v, 1])[[k-1]]+(Subscript[v, 1])[[k]]) dt/2; 
         (Subscript[du, 2])= ((Subscript[v, 2])[[k-1]]+(Subscript[v, 2])[[k]]) dt/2;
         AppendTo[(Subscript[u, 1]),(Subscript[u, 1])[[k-1]]+(Subscript[du, 1])];
         AppendTo[(Subscript[u, 2]),(Subscript[u, 2])[[k-1]]+(Subscript[du, 2])];
],k]

Define deformed coordinates, scaled for visualization

(Subscript[x, 1])={(Subscript[X, 1])}; (Subscript[x, 2])={(Subscript[X, 2])};
sc=25h/(Max[(Subscript[u, 1])[[nt]]]-Min[(Subscript[u, 1])[[nt]]]);
Monitor[
For[k=2,k≤nt,k++,
         AppendTo[(Subscript[x, 1]),(Subscript[X, 1])+sc*(Subscript[u, 1])[[k]]];
         AppendTo[(Subscript[x, 2]),(Subscript[X, 2])+sc*(Subscript[u, 2])[[k]]]; 
],k]

Define a function to plot grid lines

GridPlot[x_,y_]:= Module[{pts,xlines,ylines,mx,my,grid},
{mx,my}=Dimensions[x];
pts =Table[{x[[i,j]],y[[i,j]]},{i,1,mx,5},{j,1,my,5}]; 
xlines =Map[Line[#]&,pts];
ylines= Map[Line[#]&,Transpose[pts]];
grid=Graphics[{xlines,ylines}]
]

Animate[GridPlot[(Subscript[x, 1])[[k]],(Subscript[x, 2])[[k]]],{k,1,nt,1}]

Feedforward neural network model of material properties

Define the material properties

cpM=cp1;csM=cs1;ρM=ρ1; μ1= ρ1 (cs1^2);
μM=ρM (csM^2);λM=ρM (cpM^2)-2μM;
λdata=Table[λM,{i1,m1},{i2,m2}];
μdata=Table[μM,{i1,m1},{i2,m2}];
For[i1=1,i1≤m1,i1+=di,
For[i2=1,i2≤m2,i2+=di,
μM=ρM (csM^2);λM=ρM (cpM^2)-2μ1;
λdata[[i1;;i1+di-1,i2;;i2+di-1]]=λM;
μdata[[i1;;i1+di-1,i2;;i2+di-1]]=μM;
ρM+=dρ; cpM+=dcp;  csM+=dcs;
]
];

ArrayPlot[Transpose[ λdata],DataReversed->True ,PlotLegends->Automatic,ColorFunction->"Rainbow"]

ArrayPlot[Transpose[ μdata],DataReversed->True ,PlotLegends->Automatic,ColorFunction->"Rainbow"]

Define training, test data

Choose one time slice as training data, and another as test data

kTrain=20; kTest=30;

Compute strain and stress tensors used in DNN training

k=kTrain; d=2;
ϵ=Table[0,{j,d},{k,d},{i1,m1},{i2,m2}];
For[i1=2,i1<m1,i1++,
For[i2=2,i2<m2, i2++,
ϵ[[1,1,i1,i2]]=((Subscript[u, 1])[[k,i1+1,i2]]-(Subscript[u, 1])[[k,i1-1,i2]])/(2h);
ϵ[[1,2,i1,i2]]=((Subscript[u, 1])[[k,i1,i2+1]]-(Subscript[u, 1])[[k,i1,i2-1]]+
                                  (Subscript[u, 2])[[k,i1+1,i2]]-(Subscript[u, 2])[[k,i1-1,i2]])/(4h);
ϵ[[2,1,i1,i2]]= ϵ[[1,2,i1,i2]];
ϵ[[2,2,i1,i2]]=((Subscript[u, 2])[[k,i1,i2+1]]-(Subscript[u, 2])[[k,i1,i2-1]])/(2h);
]];
sig={{(Subscript[σ, 11])[[k]],(Subscript[σ, 12])[[k]]},{(Subscript[σ, 12])[[k]],(Subscript[σ, 22])[[k]]}};
TrainingInput = Flatten[Table[{(Subscript[σ, 11])[[kTrain,i1,i2]],(Subscript[σ, 22])[[kTrain,i1,i2]],(Subscript[σ, 12])[[kTrain,i1,i2]],ϵ[[1,1,i1,i2]],ϵ[[2,2,i1,i2]],ϵ[[1,2,i1,i2]]}
,{i1,m1},{i2,m2}],1];
TrainingOutput=Flatten[Table[λdata[[i1,i2]],{i1,m1},{i2,m2}]];
nTrainingCases=Length[TrainingInput]

Normalize the training data

{MeanInput,SDInput}={Mean[TrainingInput],StandardDeviation[TrainingInput]};
TableForm[{MeanInput,SDInput}]

{MeanOutput,SDOutput}={Mean[TrainingOutput],StandardDeviation[TrainingOutput]};
TableForm[{MeanOutput,SDOutput}]

StandardizedTrainingInput =Table[ (TrainingInput[[k]]-MeanInput)/SDInput,{k,nTrainingCases}];
StandardizedTrainingOutput=(TrainingOutput-MeanOutput)/SDOutput;

TrainingData = Table[StandardizedTrainingInput[[k]]->StandardizedTrainingOutput[[k]],{k,nTrainingCases}];

Sample of training data

TableForm[RandomSample[TrainingData,3]]

Compute strain and stress tensors used in DNN testing

k=kTest; d=2;
ϵ=Table[0,{j,d},{k,d},{i1,m1},{i2,m2}];
For[i1=2,i1<m1,i1++,
For[i2=2,i2<m2, i2++,
ϵ[[1,1,i1,i2]]=((Subscript[u, 1])[[k,i1+1,i2]]-(Subscript[u, 1])[[k,i1-1,i2]])/(2h);
ϵ[[1,2,i1,i2]]=((Subscript[u, 1])[[k,i1,i2+1]]-(Subscript[u, 1])[[k,i1,i2-1]]+
                                  (Subscript[u, 2])[[k,i1+1,i2]]-(Subscript[u, 2])[[k,i1-1,i2]])/(4h);
ϵ[[2,1,i1,i2]]= ϵ[[1,2,i1,i2]];
ϵ[[2,2,i1,i2]]=((Subscript[u, 2])[[k,i1,i2+1]]-(Subscript[u, 2])[[k,i1,i2-1]])/(2h);
]];
sig={{(Subscript[σ, 11])[[k]],(Subscript[σ, 12])[[k]]},{(Subscript[σ, 12])[[k]],(Subscript[σ, 22])[[k]]}};
TestInput = Flatten[Table[{(Subscript[σ, 11])[[kTrain,i1,i2]],(Subscript[σ, 22])[[kTrain,i1,i2]],(Subscript[σ, 12])[[kTrain,i1,i2]],ϵ[[1,1,i1,i2]],ϵ[[2,2,i1,i2]],ϵ[[1,2,i1,i2]]}
,{i1,m1},{i2,m2}],1];
TestOutput=Flatten[Table[λdata[[i1,i2]],{i1,m1},{i2,m2}]];
nTestCases=Length[TestInput]

Transform the test data according to standardization of training data

StandardizedTestInput =Table[ (TestInput[[k]]-MeanInput)/SDInput,{k,nTestCases}];
StandardizedTestOutput=(TestOutput-MeanOutput)/SDOutput;

TestData = Table[StandardizedTestInput[[k]]->StandardizedTestOutput[[k]],{k,nTestCases}];

Single linear layer model

Extract the dimensions of input, output data of the training data

in=First[TrainingData[[1]]] ;nIn=If[NumberQ[in],1,Length[in]];
out=Last[TrainingData[[1]]];
nOut=If[NumberQ[out],1,Length[out]];
{nIn,nOut}

Train the layer for a subset of the available training data (to reduce computational time), which for a linear layer is simply the linear least squares approximation

nSubset=200;
SLM=NetTrain[LinearLayer[],RandomSample[TrainingData,nSubset]]

Test the accuracy of the linear regression on another subset of the training data

subset=RandomSample[TrainingData,nSubset];
exact = Map[Last[#]&,subset];
apprx=Map[SLM[First[#]]&,subset];
error=Norm[exact-apprx]/Norm[exact]

Try a larger subset

nSubset=2000;
SLM=NetTrain[LinearLayer[],RandomSample[TrainingData,nSubset]]
subset=RandomSample[TrainingData,nSubset];
exact = Map[Last[#]&,subset];
apprx=Map[SLM[First[#]]&,subset];
error=Norm[exact-apprx]/Norm[exact]

Try the full training data

SLM=NetTrain[LinearLayer[],TrainingData]
subset=RandomSample[TrainingData,nSubset];
exact = Map[Last[#]&,subset];
apprx=Map[SLM[First[#]]&,subset];
error=Norm[exact-apprx]/Norm[exact]

Evaluate the error in prediction of the test data

subset=RandomSample[TestData,nSubset];
exact = Map[Last[#]&,subset];
apprx=Map[SLM[First[#]]&,subset];
error=Norm[exact-apprx]/Norm[exact]

Two layer model

tlm=NetChain[{nIn,LogisticSigmoid,nIn,LogisticSigmoid,nOut}]

nSubset=200;
subset=RandomSample[TrainingData,nSubset];
TLM=NetTrain[tlm,subset];
exact = Map[Last[#]&,subset];
apprx=Map[TLM[First[#]]&,subset];
error=Norm[exact-apprx]/Norm[exact]

n-layer model

nLayers=6;
nlm=NetChain[Flatten[{Table[{nIn,LogisticSigmoid},{n,nLayers}],nOut}]]

nSubset=200;
subset=RandomSample[TrainingData,nSubset];
NLM=NetTrain[nlm,subset];
exact = Map[Last[#]&,subset];
apprx=Map[NLM[First[#]]&,subset];
error=Norm[exact-apprx]/Norm[exact]

Evaluate on a different data subset

nSubset=200;
subset=RandomSample[TrainingData,nSubset];
exact = Map[Last[#]&,subset];
apprx=Map[NLM[First[#]]&,subset];
error=Norm[exact-apprx]/Norm[exact]

Train on the full data set

NLM=NetTrain[nlm,TrainingData];
exact = Map[Last[#]&,TrainingData];
apprx=Map[NLM[First[#]]&,TrainingData];
error=Norm[exact-apprx]/Norm[exact]

Evaluate on the test data (different simulation data time slice)

exact = Map[Last[#]&,TestData];
apprx=Map[NLM[First[#]]&,TestData];
error=Norm[exact-apprx]/Norm[exact]

Higher dimension hidden layers

nLayers=3;
dlm=NetChain[Flatten[{Table[{nIn+n-1,LogisticSigmoid},{n,nLayers}],Table[{nIn+n-1,LogisticSigmoid},{n,nLayers,1,-1}],nOut}]]

nSubset=200;
subset=RandomSample[TrainingData,nSubset];
DLM=NetTrain[dlm,subset];
exact = Map[Last[#]&,subset];
apprx=Map[DLM[First[#]]&,subset];
error=Norm[exact-apprx]/Norm[exact]

Train on the full data set

DLM=NetTrain[dlm,TrainingData];
exact = Map[Last[#]&,TrainingData];
apprx=Map[DLM[First[#]]&,TrainingData];
error=Norm[exact-apprx]/Norm[exact]

Evaluate on the test data (different simulation data time slice)

exact = Map[Last[#]&,TestData];
apprx=Map[DLM[First[#]]&,TestData];
error=Norm[exact-apprx]/Norm[exact]

Classifier neural network model of material properties

Extract the finite set of possible moduli

λclasses = (DeleteDuplicates[Flatten[λdata]]-MeanOutput)/SDOutput;
nclasses=Length[λclasses]

TrainingClasses=Table[TrainingData[[k,1]]->Position[λclasses,Last[TrainingData[[k]]]][[1,1]],{k,nTrainingCases}];

Shallow network

dHidden=2nclasses;
cls=NetChain[{dHidden,Ramp,nclasses,SoftmaxLayer[]},"Input"->nIn,"Output"->NetDecoder[{"Class",Table[k,{k,nclasses}]}]]

Train on the full data set

CLS=NetTrain[cls,TrainingClasses];

sample=RandomSample[TrainingClasses,10]

Map[CLS[#[[1]]]&,sample]

exact = Map[Last[#]&,TrainingClasses];
apprx=Map[CLS[First[#]]&,TrainingClasses];
error=N[Norm[exact-apprx]/Norm[exact]]