Diff of /arm_model/analysis.py [000000] .. [17a672]

Switch to unified view

a b/arm_model/analysis.py
1
import numpy as np
2
import pylab as plt
3
import matplotlib as mpl
4
# from mpl_toolkits.mplot3d import Axes3D
5
import pandas as pd
6
from util import to_np_mat, plot_corr_ellipses,  \
7
    convex_bounded_vertex_enumeration, nullspace
8
from logger import Logger
9
10
11
# ------------------------------------------------------------------------
12
# FeasibleMuscleSetAnalysis
13
# ------------------------------------------------------------------------
14
15
def construct_muscle_space_inequality(NR, fm_par, Fmax):
16
    """Construct the feasible muscle space Z f_m0 <= B .
17
18
    Parameters
19
    ----------
20
21
    NR: moment arm null space matrix
22
23
    fm_par: particular muscle forces
24
25
    Fmax: maximum muscle force
26
27
    """
28
    Z0 = -NR
29
    Z1 = NR
30
    B0 = fm_par
31
    B1 = np.asmatrix(np.diag(Fmax)).reshape(Fmax.shape[0], 1) - fm_par
32
    Z = np.concatenate((Z0, Z1), axis=0)
33
    B = np.concatenate((B0, B1), axis=0)
34
    return Z, B
35
36
37
class FeasibleMuscleSetAnalysis:
38
    """Feasible muscle set analysis.
39
40
    The required command along with the state of the system are recorded. Then
41
    this information is used to compute the feasible muscle null space and
42
    visualize it.
43
44
    """
45
46
    def __init__(self, model, simulation_reporter):
47
        """
48
        """
49
        self.logger = Logger('FeasibleMuscleSetAnalsysis')
50
        self.model = model
51
        self.simulation_reporter = simulation_reporter
52
53
    def visualize_simple_muscle(self, t, ax=None):
54
        """Visualize the feasible force set at a particular time instance for a linear
55
        muscle.
56
57
        Parameters
58
        ----------
59
60
        t: time
61
62
        ax: 1 x 3 axis
63
64
        """
65
        m = self.model.md
66
        q, Z, B, NR, fm_par = self.calculate_simple_muscles(t)
67
        x_max = np.max(to_np_mat(self.model.Fmax))
68
        fm_set = self.generate_solutions(Z, B, NR, fm_par)
69
        dataframe = pd.DataFrame(fm_set, columns=['$m_' + str(i) + '$' for i in
70
                                                  range(1, m + 1)])
71
72
        # box plot
73
        if ax is None or ax.shape[0] < 3:
74
            fig, ax = plt.subplots(1, 3, figsize=(15, 5))
75
76
        # box plot
77
        dataframe.plot.box(ax=ax[0])
78
        ax[0].set_xlabel('muscle id')
79
        ax[0].set_ylabel('force $(N)$')
80
        ax[0].set_title('Muscle-Force Box Plot')
81
        ax[0].set_ylim([0, 1.1 * x_max])
82
83
        # correlation matrix
84
        cmap = mpl.cm.jet
85
        norm = mpl.colors.Normalize(vmin=-1, vmax=1)
86
        corr = dataframe.corr()
87
        m = plot_corr_ellipses(corr, ax=ax[1], norm=norm, cmap=cmap)
88
        cb = plt.colorbar(m, ax=ax[1], orientation='vertical', norm=norm,
89
                          cmap=cmap)
90
        cb.set_label('Correlation Coefficient')
91
        ax[1].margins(0.1)
92
        ax[1].set_xlabel('muscle id')
93
        ax[1].set_ylabel('muscle id')
94
        ax[1].set_title('Correlation Matrix')
95
        ax[1].axis('equal')
96
97
        # draw model
98
        self.model.draw_model(q, False, ax[2], scale=0.7, text=False)
99
100
    def calculate_simple_muscles(self, t):
101
        """Construct Z f_m0 <= B for the case of a linear muscle model for a particular
102
        time instance.
103
104
        Parameters
105
        ----------
106
107
        t: time
108
109
        """
110
        # find nearesrt index corresponding to t
111
        idx = np.abs(np.array(self.simulation_reporter.t) - t).argmin()
112
        t = self.simulation_reporter.t[idx]
113
        q = self.simulation_reporter.q[idx]
114
        u = self.simulation_reporter.u[idx]
115
        tau = self.simulation_reporter.tau[idx]
116
        pose = self.model.model_parameters(q=q, u=u)
117
        n = self.model.nd
118
119
        # calculate required variables
120
        R = to_np_mat(self.model.R.subs(pose))
121
        RBarT = np.asmatrix(np.linalg.pinv(R.T))
122
        # reduce to independent columns to avoid singularities (proposition 3)
123
        NR = nullspace(R.transpose())
124
        fm_par = np.asmatrix(-RBarT * tau.reshape((n, 1)))
125
        Fmax = to_np_mat(self.model.Fmax)
126
127
        Z, B = construct_muscle_space_inequality(NR, fm_par, Fmax)
128
129
        return q, Z, B, NR, fm_par
130
131
    def generate_solutions(self, A, b, NR, fm_par):
132
        """Sample the solution space that satisfy A x <= b.
133
134
        Parameters
135
        ----------
136
137
        A: matrix A
138
139
        b: column vector
140
141
        NR: moment arm nullspace
142
143
        fm_par: particular solution
144
145
        Returns
146
        -------
147
148
        muscle forces: a set of solutions that satisfy the problem
149
150
        """
151
        feasible_set = []
152
        fm0_set = convex_bounded_vertex_enumeration(np.array(A),
153
                                                    np.array(b).flatten(), 0)
154
        n = fm0_set.shape[0]
155
        for i in range(0, n):
156
            fm = fm_par + NR * np.matrix(fm0_set[i, :]).reshape(-1, 1)
157
            feasible_set.append(fm)
158
159
        return np.array(feasible_set).reshape(n, -1)
160
161
162
def test_feasible_set(model):
163
    feasible_set = FeasibleMuscleSetAnalysis(model)
164
    n = model.nd
165
    m = model.md
166
    feasible_set.record(1,
167
                        np.random.random((m, 1)),
168
                        np.random.random((m, m)),
169
                        np.random.random((n, 1)))
170
    fig, ax = plt.subplots(2, 3, figsize=(10, 10))
171
    feasible_set.visualize_simple_muscle(1, ax[0])
172
    feasible_set.visualize_simple_muscle(1, ax[1])
173
    plt.show()