1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
| from __future__ import division
from os import getcwd, path as op
from numpy import array, sqrt, mean, hstack, max, zeros
from anypytools.abcutils import AnyPyProcess
from scipy.optimize import fmin_l_bfgs_b, fmin_tnc
global iteration_counter
iteration_counter = 0
def setup_study(designvars):
if designvars is None:
designvars = array([1.0]*5)
basepath = op.join( getcwd(), '..')
app = AnyPyProcess(basepath, anybodycon_path = "C:/Program Files/AnyBody Technology/AnyBody.6.0/anybodycon.exe",
num_processes = 3, disp = False)
loadmacro = 'load "ArmStrength.main.any" -def WriteOutputFiles=0'
macrocmds = ['operation RunApplication', 'run']
inputs = [('Main.DesignVars.ElbowFlexionCalFactor', designvars[0]),
('Main.DesignVars.ElbowExtensionCalFactor', designvars[1]),
('Main.DesignVars.ElbowExtension',designvars[2]),
('Main.DesignVars.ElbowPronation',designvars[3]),
('Main.DesignVars.ElbowSupination',designvars[4])]
outputs = ['Main.Studies.ElbowFlexionParamStudy.Output.Strength.Val',
'Main.Studies.ElbowExtensionParamStudy.Output.Strength.Val',
'Main.Studies.ElbowPronationParamStudy.Output.Strength.Val',
'Main.Studies.ElbowSupinationParamStudy.Output.Strength.Val']
return (app, loadmacro, macrocmds, inputs,outputs)
def evalualte(designvars,*arg):
global iteration_counter
iteration_counter += 1
print iteration_counter, 'Design vars:', designvars
#Setup and evalute objective function
(app, loadmacro, macrocmds, inputs,outputs) = setup_study(designvars)
res = app.start_param_job(loadmacro, macrocmds, inputs, outputs)
# Keys
key_flex = outputs[0]; key_ext = outputs[1];
key_pro = outputs[2];key_sup = outputs[3]
objective = metric(res[key_flex],res[key_ext],res[key_pro],res[key_sup])
print 'Objective:', objective
return objective
def eval_w_gradient(designvars,*arg):
pertubation_factor = arg[0]
global iteration_counter
iteration_counter += 1
print iteration_counter, 'Design vars:', designvars
#Setup ond run study
(app, loadmacro, macrocmds, inputs,outputs) = setup_study(designvars)
(res, pert) = app.start_pertubation_job(loadmacro, macrocmds, inputs, outputs,
perturb_factor = pertubation_factor)
# Keys
key_flex = outputs[0]; key_ext = outputs[1];
key_pro = outputs[2];key_sup = outputs[3]
pertubations = zeros((len(designvars,)) )
for i in range(len(designvars)):
pertubations[i] = metric(pert[key_flex][i], pert[key_ext][i],
pert[key_pro][i], pert[key_sup][i])
objective = metric(res[key_flex],res[key_ext],res[key_pro],res[key_sup])
gradient = (pertubations-objective)/pertubation_factor
print 'Objective:', objective
print 'Gradient:', gradient
return objective, gradient
def gradient(designvars, *arg):
res, grad = eval_w_gradient(designvars, arg)
return grad
def metric(c_flex, c_ext, c_pro, c_sup):
"""Calculate the metric (objective value) based on Fabians experimental
data
Keyword arguments:
c_flex -- computed result for elbow flexion strength
c_ext -- computed result for elbow extension strength
c_pro -- computed result for elbow pronation strength
c_sup -- computed result for elbow supination strength
"""
if isinstance(c_flex, list):
c_flex = c_flex[0]
if isinstance(c_ext, list):
c_ext = c_ext[0]
if isinstance(c_pro, list):
c_pro = c_pro[0]
if isinstance(c_sup, list):
c_sup = c_sup[0]
import experimental_data as data
# Get mean normalized experimental data
m_flex = mean(data.flexion_normalized,1)
m_ext = mean(data.extension_normalized,1)
m_pro = mean(data.pronation_normalized,1)
m_sup = mean(data.supination_normalized,1)
# Normalize computed strength arrays by strength in 90 e.flex, 0deg e. pro
norm_factor = c_flex[6]
c_flex = c_flex/norm_factor
c_ext = c_ext/norm_factor
c_pro = c_pro/norm_factor
c_sup = c_sup/norm_factor
# Calculate the RMS errors
try:
rms_flex = sqrt( sum(((c_flex-m_flex)/m_flex)**2))
rms_ext = sqrt( sum(((c_ext-m_ext)/m_ext)**2) )
rms_pro = sqrt( sum(((c_pro-m_pro)/m_pro)**2))
rms_sup = sqrt( sum(((c_sup-m_sup)/m_sup)**2))
except TypeError:
print 'Error in metric. Return 100'
return 100
# Calucate Weight for the metrics
N_total = len(m_pro)+len(m_sup)+len(m_ext)+len(m_flex)
weight = float(N_total) / array([len(m_flex), len(m_ext), len(m_pro),
len(m_sup) ] )
return rms_flex*weight[0] + rms_ext*weight[1] + rms_sup*weight[2] + rms_pro*weight[3]
if __name__ == '__main__':
# e1 = eval_strength(array([1.2,1.,1.,1.,1.]),(1e-3))
# e2 = eval_strength(array([1.,1.,1.,1.,1.]))
x0 = array([1.0, 1.0, 1.0, 1.0, 1.0])
epsilon = 1e-4 # pertubation to find gradient
(x,f,d) = fmin_l_bfgs_b(func = eval_w_gradient,
x0 = array([1.0, 1.0, 1.0, 1.0, 1.0]),
args=(epsilon,),
bounds = [(0.01,1.8),(0.01,1.8),(0.1,4), (0.1,4), (0.1,4)] ) |
Procédure d'optimisation Python
Répondre avec citation
Partager