Using ParamFitter
Once a sample has been created, one may want to fit the parameters to match some observed data. This is facilitated by ParamFitter class. It uses scipy’s least_square function to fit. However, internally it takes care of converting between the ParamDict and numpy vector as well as only varying the non-fixed parameters.
[1]:
import pymecht as pmt
import numpy as np
from matplotlib import pyplot as plt
pmt.ParamFitter?
Init signature: pmt.ParamFitter(sim_func, output, params=None)
Docstring:
A class to fit parameters of a simulation function to a given output
Parameters
----------
sim_func : function
A function that takes in a dictionary of parameters and returns a numpy array of the same shape as the output
output : numpy array
A numpy array of the same shape as the output of the sim_func to be fitted
params : ParamDict
A dictionary of parameters to be varied
File: /usr/local/lib/python3.9/site-packages/pymecht/ParamFitter.py
Type: type
Subclasses:
As we can see, creating an instance of ParamFitter requires three main things:
a Python function that takes in parameters, perform the simulation and returns the modeled output
the observed output array
parameters
The ability to provide a function for the simulation provides flexibility in how the problem can be setup. In a simple setting it could be just performing a single simulation of a single sample. In more complex scenarios, one can perform multiple simulations of multiple samples and combine their output, potentially through a weighting function. We can start with a simple case
Simple fitting
Assume a uniaxial extension experiment, which was force controlled. That is the applied stress was varied linearly and the stretch was observed.
[2]:
mat = pmt.MatModel('yeoh')
sample = pmt.UniaxialExtension(mat, force_measure='cauchy')
params = sample.parameters
print(sample, params)
applied_stress = np.linspace(0,50,10)
observed_stretch = np.array([1., 1.09144422, 1.17650268, 1.25077441, 1.31313774, 1.36499585,
1.40850493, 1.4455964, 1.4777529, 1.50606163])
def simulate(param_vary):
mod_stretch = sample.force_controlled(applied_stress,param_vary) #make sure that the argument param_vary is used
return mod_stretch
An object of type UniaxialExtensionwith stretch as input, cauchy as output, and the following material
Material model with 1 component:
Component1: YEOH
------------------------------------------------------------------
Keys Value Fixed? Lower bound Upper bound
------------------------------------------------------------------
L0 1.00 No 1.00e-04 1.00e+03
A0 1.00 No 1.00e-04 1.00e+03
c1_0 1.00 No 1.00e-04 1.00e+02
c2_0 1.00 No 0.00 1.00e+02
c3_0 1.00 No 0.00 1.00e+02
c4_0 0.00 No 0.00 1.00e+02
------------------------------------------------------------------
For stress and stretch, the length and cross-sectional area would not be relevant. Thus, it is best to fix them. We will also fix c4_0. Also, note that when creating the instance of ParamFitter it will call the simulation function once to check that it returns an array of same size as the observations. Otherwise, it will raise an error.
[3]:
params.fix('L0')
params.fix('A0')
params.fix('c4_0')
try:
param_fitter = pmt.ParamFitter(simulate,observed_stretch[:-1],params)
except ValueError as e:
print("Error occured", e)
param_fitter = pmt.ParamFitter(simulate,observed_stretch,params)
Error occured ('Output and simulation results are of different shape', (10,), (9,))
Parameter fitting instance created with the following settings
------------------------------------------------------------------
Keys Value Fixed? Lower bound Upper bound
------------------------------------------------------------------
L0 1.00 Yes - -
A0 1.00 Yes - -
c1_0 1.00 No 1.00e-04 1.00e+02
c2_0 1.00 No 0.00 1.00e+02
c3_0 1.00 No 0.00 1.00e+02
c4_0 0.00 Yes - -
------------------------------------------------------------------
3 parameters will be fitted.
Once the ParamFitter instance has been created successfully, one can just do .fit() to perform the fitting.
[4]:
param_fitter.fit()
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 4.5172e-01 2.56e+01
1 2 1.6712e-01 2.85e-01 1.73e+00 4.88e+00
2 3 7.2197e-02 9.49e-02 1.72e+00 2.42e+00
3 4 8.9118e-03 6.33e-02 3.43e+00 5.76e-01
4 5 3.2355e-04 8.59e-03 4.96e+00 2.38e-03
5 6 6.3679e-05 2.60e-04 1.34e+00 1.68e-03
6 7 4.2171e-05 2.15e-05 7.10e-01 2.86e-03
7 8 1.1726e-06 4.10e-05 2.70e+00 6.08e-04
8 9 1.7271e-10 1.17e-06 5.48e-01 2.90e-05
9 10 3.0903e-16 1.73e-10 5.55e-03 8.86e-09
`gtol` termination condition is satisfied.
Function evaluations 10, initial cost 4.5172e-01, final cost 3.0903e-16, first-order optimality 8.86e-09.
Fitting completed, with the following results
------------------------------------------------------------------
Keys Value Fixed? Lower bound Upper bound
------------------------------------------------------------------
L0 1.00 Yes - -
A0 1.00 Yes - -
c1_0 10.00 No 1.00e-04 1.00e+02
c2_0 2.00 No 0.00 1.00e+02
c3_0 3.00 No 0.00 1.00e+02
c4_0 0.00 Yes - -
------------------------------------------------------------------
[4]:
message: `gtol` termination condition is satisfied.
success: True
status: 1
fun: [ 0.000e+00 5.430e-09 3.199e-09 3.850e-10 -2.325e-09
-1.012e-08 -5.783e-09 2.076e-08 -2.431e-09 3.523e-10]
x: [ 1.000e+01 2.000e+00 3.000e+00]
cost: 3.0902897658340763e-16
jac: [[ 0.000e+00 0.000e+00 0.000e+00]
[-8.833e-03 -4.184e-04 -1.487e-05]
...
[-1.632e-02 -1.753e-02 -1.413e-02]
[-1.540e-02 -1.837e-02 -1.642e-02]]
grad: [-8.290e-11 -8.063e-11 -9.129e-11]
optimality: 8.855302470909847e-09
active_mask: [0 0 0]
nfev: 10
njev: 10
The results have converged, and we can plot to make sure that it is a good fit.
[5]:
fig,ax = plt.subplots(1,1,figsize=(4*0.7,3*0.7))
ax.plot(applied_stress, observed_stretch,'o',label='Observed data')
ax.plot(applied_stress, simulate(params),'-',label='Fitted model')
ax.set_xlabel('Stress')
ax.set_ylabel('Stretch')
plt.show()
Adding fiber direction
For an anisotropic material, one has to specify the fiber direction. It could be included as a fitting parameter as shown below.
[6]:
mat = pmt.MatModel('goh','nh')
sample = pmt.UniaxialExtension(mat, force_measure='cauchy')
params = sample.parameters
params['phi'] = pmt.Param(30,0,90) #Add fiber angle as a parameter
params.fix('L0')
params.fix('A0')
print(sample, params)
applied_stress = np.linspace(0,50,10)
observed_stretch = np.array([1., 1.09144422, 1.17650268, 1.25077441, 1.31313774, 1.36499585,
1.40850493, 1.4455964, 1.4777529, 1.50606163])
def simulate(param_vary):
pmt.specify_two_fibers(sample,param_vary['phi'],verbose=False)
mod_stretch = sample.force_controlled(applied_stress, param_vary)
return mod_stretch
param_fitter = pmt.ParamFitter(simulate,observed_stretch,params)
param_fitter.fit()
An object of type UniaxialExtensionwith stretch as input, cauchy as output, and the following material
Material model with 2 components:
Component1: GOH
Component2: NH
------------------------------------------------------------------
Keys Value Fixed? Lower bound Upper bound
------------------------------------------------------------------
L0 1.00 Yes - -
A0 1.00 Yes - -
k1_0 10.00 No 0.10 30.00
k2_0 10.00 No 0.10 30.00
k3_0 0.10 No 0.00 0.33
mu_1 1.00 No 1.00e-04 1.00e+02
phi 30.00 No 0.00 90.00
------------------------------------------------------------------
Parameter fitting instance created with the following settings
------------------------------------------------------------------
Keys Value Fixed? Lower bound Upper bound
------------------------------------------------------------------
L0 1.00 Yes - -
A0 1.00 Yes - -
k1_0 10.00 No 0.10 30.00
k2_0 10.00 No 0.10 30.00
k3_0 0.10 No 0.00 0.33
mu_1 1.00 No 1.00e-04 1.00e+02
phi 30.00 No 0.00 90.00
------------------------------------------------------------------
5 parameters will be fitted.
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 6.4881e-02 3.41e-01
1 2 2.0477e-02 4.44e-02 7.98e+00 2.23e-01
2 3 3.1627e-03 1.73e-02 1.30e+01 7.88e-02
3 4 2.3278e-04 2.93e-03 3.78e+00 1.56e-02
4 5 1.0941e-05 2.22e-04 1.74e+00 2.89e-03
5 6 5.6819e-06 5.26e-06 3.31e-01 3.30e-05
6 7 5.4415e-06 2.40e-07 9.64e-01 1.49e-03
7 8 5.0847e-06 3.57e-07 1.09e-01 8.32e-06
8 12 5.0731e-06 1.15e-08 2.76e-01 1.21e-04
9 13 5.0615e-06 1.16e-08 5.23e-01 3.40e-04
10 14 5.0293e-06 3.22e-08 4.37e-01 2.27e-04
11 15 5.0103e-06 1.90e-08 5.23e-01 3.28e-04
12 16 4.9793e-06 3.10e-08 4.55e-01 2.47e-04
13 17 4.9674e-06 1.19e-08 6.03e-01 4.46e-04
14 18 4.9222e-06 4.52e-08 3.37e-01 1.36e-04
15 20 4.9047e-06 1.75e-08 2.52e-01 8.61e-05
16 21 4.8908e-06 1.39e-08 5.08e-01 3.65e-04
17 22 4.8544e-06 3.64e-08 4.17e-01 2.26e-04
18 23 4.8369e-06 1.75e-08 5.31e-01 3.88e-04
19 24 4.7974e-06 3.95e-08 3.92e-01 2.08e-04
20 25 4.7828e-06 1.46e-08 5.32e-01 4.19e-04
21 26 4.7395e-06 4.33e-08 3.62e-01 1.86e-04
22 27 4.7291e-06 1.04e-08 5.32e-01 4.59e-04
23 28 4.6809e-06 4.82e-08 3.30e-01 1.60e-04
24 29 4.6769e-06 3.97e-09 5.30e-01 5.14e-04
25 30 4.6277e-06 4.93e-08 1.32e-01 3.23e-05
26 31 4.6146e-06 1.31e-08 2.25e-01 3.44e-04
27 32 4.5798e-06 3.48e-08 2.68e-01 1.35e-04
28 34 4.5626e-06 1.72e-08 1.17e-01 7.79e-05
29 35 4.5445e-06 1.82e-08 2.37e-01 3.01e-04
30 36 4.5260e-06 1.84e-08 5.02e-01 4.43e-04
31 37 4.4804e-06 4.57e-08 3.37e-01 1.96e-04
32 39 4.4599e-06 2.04e-08 1.22e-01 7.61e-05
33 40 4.4451e-06 1.48e-08 2.38e-01 3.68e-04
34 41 4.4084e-06 3.66e-08 2.72e-01 1.51e-04
35 43 4.3891e-06 1.94e-08 1.19e-01 9.75e-05
36 44 4.3717e-06 1.74e-08 2.48e-01 3.36e-04
37 45 4.3384e-06 3.33e-08 2.72e-01 1.60e-04
38 47 4.3186e-06 1.98e-08 1.21e-01 9.85e-05
39 48 4.3019e-06 1.67e-08 2.51e-01 3.44e-04
40 49 4.2686e-06 3.33e-08 2.72e-01 1.61e-04
41 51 4.2489e-06 1.97e-08 1.23e-01 1.01e-04
42 52 4.2325e-06 1.64e-08 2.54e-01 3.40e-04
43 53 4.2003e-06 3.22e-08 2.72e-01 1.61e-04
44 55 4.1809e-06 1.94e-08 1.24e-01 1.00e-04
45 56 4.1650e-06 1.59e-08 2.55e-01 3.35e-04
46 57 4.1341e-06 3.09e-08 2.71e-01 1.59e-04
47 59 4.1153e-06 1.88e-08 1.24e-01 9.79e-05
48 60 4.0999e-06 1.54e-08 2.56e-01 3.27e-04
49 61 4.0705e-06 2.94e-08 2.70e-01 1.56e-04
50 63 4.0525e-06 1.80e-08 1.25e-01 9.43e-05
51 64 4.0376e-06 1.49e-08 2.56e-01 3.25e-04
52 65 4.0099e-06 2.77e-08 2.69e-01 1.60e-04
53 67 3.9928e-06 1.71e-08 1.25e-01 9.61e-05
54 68 3.9784e-06 1.44e-08 2.56e-01 3.40e-04
55 69 3.9525e-06 2.59e-08 2.67e-01 1.68e-04
56 71 3.9364e-06 1.61e-08 1.25e-01 9.76e-05
57 72 3.9225e-06 1.40e-08 2.54e-01 3.52e-04
58 73 3.9023e-06 2.01e-08 3.95e-01 3.56e-04
59 74 3.8828e-06 1.95e-08 3.94e-01 3.60e-04
60 75 3.8636e-06 1.92e-08 3.92e-01 3.63e-04
61 76 3.8452e-06 1.85e-08 3.93e-01 3.69e-04
62 77 3.8268e-06 1.84e-08 3.88e-01 3.66e-04
63 78 3.8096e-06 1.72e-08 3.94e-01 3.82e-04
64 79 3.7915e-06 1.81e-08 3.80e-01 3.58e-04
65 80 3.7767e-06 1.48e-08 4.07e-01 4.18e-04
66 81 3.7571e-06 1.96e-08 3.51e-01 3.11e-04
67 82 3.7509e-06 6.20e-09 4.71e-01 5.79e-04
68 83 3.7238e-06 2.71e-08 2.56e-01 1.64e-04
69 85 3.7113e-06 1.24e-08 1.15e-01 1.37e-04
70 86 3.6987e-06 1.26e-08 2.43e-01 3.34e-04
71 87 3.6902e-06 8.55e-09 4.47e-01 5.38e-04
72 88 3.6686e-06 2.15e-08 2.80e-01 2.05e-04
73 90 3.6577e-06 1.09e-08 1.16e-01 1.07e-04
74 91 3.6462e-06 1.15e-08 2.34e-01 4.02e-04
75 92 3.6328e-06 1.34e-08 3.76e-01 3.87e-04
76 93 3.6209e-06 1.20e-08 3.91e-01 4.22e-04
77 94 3.6073e-06 1.36e-08 3.58e-01 3.53e-04
78 95 3.5983e-06 8.93e-09 4.26e-01 5.08e-04
79 96 3.5824e-06 1.59e-08 2.97e-01 2.42e-04
80 97 3.5794e-06 3.00e-09 4.70e-01 6.71e-04
81 98 3.5595e-06 1.99e-08 2.24e-01 1.35e-04
82 100 3.5498e-06 9.67e-09 1.04e-01 1.56e-04
83 101 3.5399e-06 9.91e-09 2.27e-01 2.72e-04
84 102 3.5355e-06 4.40e-09 4.60e-01 6.28e-04
85 103 3.5193e-06 1.62e-08 2.39e-01 1.55e-04
86 104 3.5170e-06 2.28e-09 4.52e-01 7.25e-04
87 105 3.5006e-06 1.64e-08 1.12e-01 3.75e-05
88 106 3.4970e-06 3.68e-09 1.77e-01 8.65e-04
89 107 3.4759e-06 2.11e-08 1.73e-01 7.87e-05
90 108 3.4665e-06 9.36e-09 2.11e-01 3.57e-04
91 109 3.4597e-06 6.81e-09 4.13e-01 4.83e-04
92 110 3.4498e-06 9.94e-09 3.04e-01 2.57e-04
93 111 3.4450e-06 4.73e-09 4.36e-01 5.59e-04
94 112 3.4344e-06 1.06e-08 2.59e-01 1.83e-04
95 113 3.4300e-06 4.36e-09 4.28e-01 5.75e-04
96 114 3.4198e-06 1.03e-08 2.48e-01 1.66e-04
97 115 3.4153e-06 4.44e-09 4.21e-01 5.62e-04
98 116 3.4059e-06 9.45e-09 2.48e-01 1.66e-04
99 117 3.4014e-06 4.52e-09 4.14e-01 5.36e-04
100 118 3.3928e-06 8.53e-09 2.55e-01 1.74e-04
101 119 3.3883e-06 4.55e-09 4.07e-01 5.04e-04
102 120 3.3806e-06 7.66e-09 2.65e-01 1.87e-04
103 122 3.3759e-06 4.76e-09 1.71e-01 1.11e-04
104 123 3.3704e-06 5.50e-09 3.41e-01 4.04e-04
105 124 3.3643e-06 6.04e-09 3.21e-01 2.73e-04
106 125 3.3612e-06 3.17e-09 4.64e-01 5.77e-04
107 126 3.3538e-06 7.39e-09 2.19e-01 1.23e-04
108 128 3.3493e-06 4.45e-09 1.62e-01 1.13e-04
109 129 3.3445e-06 4.81e-09 3.27e-01 3.41e-04
110 130 3.3398e-06 4.69e-09 3.56e-01 3.29e-04
111 131 3.3354e-06 4.45e-09 3.62e-01 3.39e-04
112 132 3.3309e-06 4.47e-09 3.46e-01 3.06e-04
113 133 3.3269e-06 4.02e-09 3.75e-01 3.59e-04
114 134 3.3225e-06 4.40e-09 3.14e-01 2.49e-04
115 135 3.3196e-06 2.88e-09 4.41e-01 4.91e-04
116 136 3.3147e-06 4.96e-09 2.22e-01 1.20e-04
117 138 3.3116e-06 3.10e-09 1.48e-01 7.83e-05
118 139 3.3075e-06 4.06e-09 2.97e-01 2.71e-04
119 140 3.3044e-06 3.16e-09 3.81e-01 3.57e-04
120 141 3.3007e-06 3.64e-09 2.84e-01 1.96e-04
121 142 3.2995e-06 1.17e-09 5.00e-01 6.09e-04
122 143 3.2944e-06 5.18e-09 1.59e-01 5.92e-05
123 145 3.2913e-06 3.06e-09 1.35e-01 7.70e-05
124 146 3.2883e-06 3.00e-09 2.78e-01 2.09e-04
125 147 3.2864e-06 1.92e-09 4.35e-01 4.48e-04
126 148 3.2830e-06 3.40e-09 2.00e-01 9.24e-05
127 150 3.2811e-06 1.86e-09 1.33e-01 5.05e-05
128 151 3.2784e-06 2.74e-09 2.65e-01 1.88e-04
129 152 3.2770e-06 1.44e-09 4.44e-01 4.57e-04
130 153 3.2739e-06 3.07e-09 1.80e-01 7.31e-05
131 154 3.2736e-06 3.01e-10 5.11e-01 6.10e-04
132 155 3.2695e-06 4.10e-09 1.24e-01 3.39e-05
133 156 3.2673e-06 2.21e-09 2.45e-01 1.53e-04
134 157 3.2666e-06 6.74e-10 4.69e-01 4.94e-04
135 158 3.2637e-06 2.88e-09 1.44e-01 4.50e-05
136 159 3.2631e-06 6.26e-10 4.71e-01 5.04e-04
137 160 3.2603e-06 2.81e-09 1.31e-01 3.72e-05
138 161 3.2596e-06 6.72e-10 4.54e-01 4.62e-04
139 162 3.2572e-06 2.38e-09 1.33e-01 3.79e-05
140 163 3.2566e-06 6.23e-10 4.37e-01 4.21e-04
141 164 3.2546e-06 2.01e-09 1.35e-01 3.87e-05
142 165 3.2540e-06 5.78e-10 4.20e-01 3.84e-04
143 166 3.2523e-06 1.70e-09 1.37e-01 3.95e-05
144 167 3.2518e-06 5.35e-10 4.03e-01 3.49e-04
145 168 3.2504e-06 1.43e-09 1.38e-01 4.00e-05
146 169 3.2499e-06 4.94e-10 3.85e-01 3.16e-04
147 170 3.2487e-06 1.20e-09 1.38e-01 4.01e-05
148 171 3.2482e-06 4.53e-10 3.67e-01 2.84e-04
149 172 3.2472e-06 1.00e-09 1.39e-01 4.03e-05
150 173 3.2468e-06 4.11e-10 3.48e-01 2.55e-04
151 174 3.2460e-06 8.32e-10 1.39e-01 4.01e-05
152 175 3.2456e-06 3.68e-10 3.30e-01 2.27e-04
153 176 3.2449e-06 6.85e-10 1.39e-01 4.00e-05
154 177 3.2446e-06 3.25e-10 3.11e-01 2.01e-04
155 178 3.2440e-06 5.60e-10 1.40e-01 4.03e-05
156 179 3.2437e-06 2.83e-10 2.91e-01 1.76e-04
157 180 3.2433e-06 4.49e-10 1.38e-01 3.93e-05
158 181 3.2431e-06 2.42e-10 2.71e-01 1.52e-04
159 182 3.2427e-06 3.55e-10 1.36e-01 3.82e-05
160 184 3.2425e-06 1.58e-10 1.11e-01 2.52e-05
161 185 3.2423e-06 2.00e-10 2.16e-01 9.63e-05
162 186 3.2421e-06 2.31e-10 1.71e-01 6.01e-05
163 187 3.2420e-06 1.58e-10 2.38e-01 1.16e-04
164 188 3.2418e-06 2.00e-10 1.14e-01 2.66e-05
165 190 3.2417e-06 8.07e-11 8.92e-02 1.63e-05
166 191 3.2416e-06 1.07e-10 1.73e-01 6.16e-05
167 192 3.2415e-06 1.14e-10 1.53e-01 4.81e-05
168 193 3.2414e-06 8.98e-11 1.69e-01 5.83e-05
169 194 3.2413e-06 8.15e-11 1.27e-01 3.32e-05
170 195 3.2412e-06 5.22e-11 1.74e-01 6.19e-05
171 196 3.2412e-06 5.80e-11 8.65e-02 1.53e-05
172 198 3.2412e-06 1.82e-11 5.17e-02 5.48e-06
173 199 3.2411e-06 2.43e-11 9.89e-02 2.00e-05
174 200 3.2411e-06 2.05e-11 1.34e-01 3.69e-05
175 201 3.2411e-06 1.83e-11 5.13e-02 5.39e-06
176 202 3.2411e-06 3.40e-12 1.38e-01 3.88e-05
177 203 3.2411e-06 1.28e-11 1.87e-02 7.20e-07
178 204 3.2411e-06 1.21e-12 8.88e-02 1.61e-05
179 205 3.2411e-06 1.49e-12 2.58e-02 1.36e-06
180 206 3.2411e-06 5.34e-13 4.09e-02 3.42e-06
181 211 3.2411e-06 2.51e-14 9.86e-03 2.01e-07
182 212 3.2411e-06 6.26e-16 3.41e-05 1.02e-09
`gtol` termination condition is satisfied.
Function evaluations 212, initial cost 6.4881e-02, final cost 3.2411e-06, first-order optimality 1.02e-09.
Fitting completed, with the following results
------------------------------------------------------------------
Keys Value Fixed? Lower bound Upper bound
------------------------------------------------------------------
L0 1.00 Yes - -
A0 1.00 Yes - -
k1_0 1.07 No 0.10 30.00
k2_0 0.42 No 0.10 30.00
k3_0 4.90e-14 No 0.00 0.33
mu_1 16.24 No 1.00e-04 1.00e+02
phi 0.02 No 0.00 90.00
------------------------------------------------------------------
[6]:
message: `gtol` termination condition is satisfied.
success: True
status: 1
fun: [ 2.220e-16 1.545e-03 1.479e-04 -9.979e-04 -8.545e-04
-3.821e-06 7.319e-04 8.811e-04 3.014e-04 -9.714e-04]
x: [ 1.072e+00 4.182e-01 4.898e-14 1.624e+01 1.712e-02]
cost: 3.2410696722792836e-06
jac: [[ 0.000e+00 0.000e+00 ... 0.000e+00 0.000e+00]
[-1.507e-02 -6.118e-04 ... -4.462e-03 4.172e-07]
...
[-9.765e-02 -1.468e-01 ... -7.910e-03 4.098e-06]
[-1.011e-01 -1.735e-01 ... -7.228e-03 4.470e-06]]
grad: [-3.516e-11 -4.143e-12 3.048e-06 2.088e-12 1.666e-11]
optimality: 1.0171337805176332e-09
active_mask: [ 0 0 -1 0 0]
nfev: 212
njev: 183
It took higher number of iterations, but fitting did work. When plotting the results, if we call the simulate function with params, we will get an error. This can solved by getting a normal dictionary from params and passing that instead.
[7]:
fig,ax = plt.subplots(1,1,figsize=(4*0.7,3*0.7))
ax.plot(applied_stress, observed_stretch,'o',label='Observed data')
ax.plot(applied_stress, simulate(params),'-',label='Fitted model')
ax.set_xlabel('Stress')
ax.set_ylabel('Stretch')
plt.show()
Combining multiple experiments
One can combine multiple experiments. For example, if two samples were taken at 90 degrees angles from the sample tissue and a uniaxial extension experiment was performed on each, we may want to have the same set of parameters fit both. This is demonstrated below.
[8]:
#since the fiber directions are stored in material, it is better to create separate instances for the two
matC = pmt.MatModel('goh','nh')
matL = pmt.MatModel('goh','nh')
sample_C = pmt.UniaxialExtension(matC, force_measure='cauchy')
sample_L = pmt.UniaxialExtension(matL, force_measure='cauchy')
params = sample_L.parameters
#add fiber angle
params['phi'] = pmt.Param(30,0,90)
stretch_C = np.linspace(1,1.2,20)
stretch_L = np.linspace(1,1.5,20)
obs_stress_C = np.array([0. , 1.52571403, 3.07863692, 4.66062427, 6.27372501,
7.92021013, 9.60260466, 11.32372343, 13.08671143, 14.89508926,
16.75280476, 18.66429179, 20.63453731, 22.66915826, 24.77448991,
26.95768758, 29.22684422, 31.59112643, 34.06093242, 36.64807567,])
obs_stress_L = np.array([0. , 3.15860501, 6.32132964, 9.49192622, 12.67378974,
15.86999936, 19.08335432, 22.31640505, 25.57148019, 28.8507102 ,
32.15604801, 35.48928713, 38.85207756, 42.24593993, 45.67227786,
49.13238907, 52.62747512, 56.15865021, 59.72694895, 63.40018943])
observed_stresses = np.array([obs_stress_C,obs_stress_L])
def simulate_both(param_vary):
phi = param_vary['phi']
if type(phi) is pmt.Param:
phi = phi.value
pmt.specify_two_fibers(sample_C,phi,verbose=False)
stress_C = sample_C.disp_controlled(stretch_C,param_vary)
pmt.specify_two_fibers(sample_L,90-phi,verbose=False)
stress_L = sample_L.disp_controlled(stretch_L,param_vary)
return np.array([stress_C,stress_L])
[9]:
params.set('k1_0',20)
params.set('k2_0',5)
params.set('k3_0',0.2)
params.set('mu_1',40)
params.set('phi',10)
observed_stress = simulate_both(params)
[10]:
observed_stress
[10]:
array([[ 0. , 1.52571403, 3.07863692, 4.66062427, 6.27372501,
7.92021013, 9.60260466, 11.32372343, 13.08671143, 14.89508926,
16.75280476, 18.66429179, 20.63453731, 22.66915826, 24.77448991,
26.95768758, 29.22684422, 31.59112643, 34.06093242, 36.64807567],
[ 0. , 3.15860501, 6.32132964, 9.49192622, 12.67378974,
15.86999936, 19.08335432, 22.31640505, 25.57148019, 28.8507102 ,
32.15604801, 35.48928713, 38.85207756, 42.24593993, 45.67227786,
49.13238907, 52.62747512, 56.15865021, 59.72694895, 63.40018943]])
[11]:
params.fix('L0')
params.fix('A0')
params.fix('k2_0',5)
params.set('k1_0',10)
params.set('k3_0',0.1)
params.set('mu_1',10)
params.set('phi',30)
param_fitter = pmt.ParamFitter(simulate_both,observed_stress,params)
Parameter fitting instance created with the following settings
------------------------------------------------------------------
Keys Value Fixed? Lower bound Upper bound
------------------------------------------------------------------
L0 1.00 Yes - -
A0 1.00 Yes - -
k1_0 10.00 No 0.10 30.00
k2_0 5.00 Yes - -
k3_0 0.10 No 0.00 0.33
mu_1 10.00 No 1.00e-04 1.00e+02
phi 30.00 No 0.00 90.00
------------------------------------------------------------------
4 parameters will be fitted.
[12]:
param_fitter.fit()
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 8.5494e+03 5.13e+04
1 2 3.1698e+02 8.23e+03 2.25e+01 5.96e+03
2 3 1.2913e+01 3.04e+02 7.48e+00 2.20e+02
3 4 1.1669e+01 1.24e+00 5.30e+00 5.89e+02
4 5 1.6546e+00 1.00e+01 1.40e+00 4.46e+01
5 6 6.1012e-01 1.04e+00 2.48e+00 6.66e+01
6 7 1.3591e-01 4.74e-01 1.65e+00 2.54e+01
7 8 5.7974e-02 7.79e-02 7.77e-01 9.42e+00
8 9 4.4427e-02 1.35e-02 1.21e+00 1.07e+01
9 10 3.0427e-02 1.40e-02 2.97e-01 3.80e-01
10 12 2.5246e-02 5.18e-03 6.58e-01 4.15e+00
11 13 2.0088e-02 5.16e-03 5.80e-01 2.69e+00
12 14 1.6521e-02 3.57e-03 7.55e-01 4.67e+00
13 15 1.3138e-02 3.38e-03 4.09e-01 1.22e+00
14 17 1.1830e-02 1.31e-03 2.85e-01 6.23e-01
15 18 1.0091e-02 1.74e-03 5.70e-01 2.33e+00
16 19 8.3100e-03 1.78e-03 5.59e-01 2.02e+00
17 20 6.9170e-03 1.39e-03 5.66e-01 2.02e+00
18 21 5.7624e-03 1.15e-03 5.02e-01 1.50e+00
19 22 4.8859e-03 8.77e-04 5.82e-01 1.93e+00
20 23 4.0831e-03 8.03e-04 4.12e-01 9.04e-01
21 24 3.7214e-03 3.62e-04 7.23e-01 2.73e+00
22 25 2.9235e-03 7.98e-04 2.39e-01 2.59e-01
23 27 2.6770e-03 2.46e-04 2.64e-01 3.44e-01
24 28 2.3590e-03 3.18e-04 5.28e-01 1.29e+00
25 29 1.9873e-03 3.72e-04 3.71e-01 5.89e-01
26 30 1.8454e-03 1.42e-04 6.80e-01 1.97e+00
27 31 1.4282e-03 4.17e-04 2.07e-01 1.59e-01
28 33 1.3016e-03 1.27e-04 2.54e-01 2.61e-01
29 34 1.1444e-03 1.57e-04 5.06e-01 9.86e-01
30 35 9.4888e-04 1.95e-04 3.14e-01 3.51e-01
31 37 8.7266e-04 7.62e-05 1.78e-01 1.15e-01
32 38 7.6457e-04 1.08e-04 3.56e-01 4.47e-01
33 39 6.5920e-04 1.05e-04 5.21e-01 9.15e-01
34 40 5.2550e-04 1.34e-04 2.49e-01 1.95e-01
35 42 4.7924e-04 4.63e-05 1.73e-01 9.68e-02
36 43 4.1040e-04 6.88e-05 3.45e-01 3.73e-01
37 44 3.3917e-04 7.12e-05 4.49e-01 6.07e-01
38 45 2.6437e-04 7.48e-05 2.61e-01 1.95e-01
39 47 2.3584e-04 2.85e-05 1.56e-01 7.01e-02
40 48 1.9556e-04 4.03e-05 3.10e-01 2.72e-01
41 49 1.5507e-04 4.05e-05 4.09e-01 4.55e-01
42 50 1.1338e-04 4.17e-05 2.27e-01 1.34e-01
43 52 9.8195e-05 1.52e-05 1.38e-01 4.99e-02
44 53 7.6953e-05 2.12e-05 2.75e-01 1.95e-01
45 54 5.6112e-05 2.08e-05 3.41e-01 2.89e-01
46 55 3.7182e-05 1.89e-05 2.03e-01 9.85e-02
47 56 3.1294e-05 5.89e-06 4.00e-01 3.84e-01
48 57 1.4213e-05 1.71e-05 9.76e-02 2.10e-02
49 59 1.0261e-05 3.95e-06 1.31e-01 3.99e-02
50 60 6.4072e-06 3.85e-06 2.61e-01 1.55e-01
51 61 2.6049e-06 3.80e-06 9.99e-02 2.18e-02
52 62 1.8644e-06 7.40e-07 2.39e-01 1.27e-01
53 63 1.8686e-07 1.68e-06 3.36e-02 2.32e-03
54 64 3.8722e-08 1.48e-07 9.57e-02 1.98e-02
55 65 7.7660e-11 3.86e-08 4.40e-03 3.89e-05
56 66 8.6566e-15 7.77e-11 2.09e-03 9.42e-06
57 67 5.2376e-24 8.66e-15 2.46e-06 1.49e-11
`gtol` termination condition is satisfied.
Function evaluations 67, initial cost 8.5494e+03, final cost 5.2376e-24, first-order optimality 1.49e-11.
Fitting completed, with the following results
------------------------------------------------------------------
Keys Value Fixed? Lower bound Upper bound
------------------------------------------------------------------
L0 1.00 Yes - -
A0 1.00 Yes - -
k1_0 20.00 No 0.10 30.00
k2_0 5.00 Yes - -
k3_0 0.20 No 0.00 0.33
mu_1 40.00 No 1.00e-04 1.00e+02
phi 10.00 No 0.00 90.00
------------------------------------------------------------------
[12]:
message: `gtol` termination condition is satisfied.
success: True
status: 1
fun: [ 0.000e+00 4.063e-13 ... -4.974e-14 1.918e-13]
x: [ 2.000e+01 2.000e-01 4.000e+01 1.000e+01]
cost: 5.2375861683058635e-24
jac: [[ 0.000e+00 0.000e+00 0.000e+00 0.000e+00]
[ 1.313e-02 -3.861e+00 3.158e-02 -4.021e-03]
...
[ 0.000e+00 0.000e+00 1.493e+00 0.000e+00]
[ 3.343e-03 3.988e+01 1.583e+00 1.050e-01]]
grad: [-2.805e-13 7.453e-11 -1.090e-13 1.037e-13]
optimality: 1.4906982976095317e-11
active_mask: [0 0 0 0]
nfev: 67
njev: 58