From a2111ac26700b95cbf44fef4b297a1b44f7c7987 Mon Sep 17 00:00:00 2001 From: stwilliams333 Date: Mon, 12 Jun 2023 11:57:04 -0500 Subject: [PATCH 1/5] Update SynX bounds --- .../calcMuscleExcitationsSynX.m | 90 +++++++++---------- .../getLinearInequalityConstraints.m | 33 +++++-- 2 files changed, 69 insertions(+), 54 deletions(-) diff --git a/src/MuscleTendonPersonalization/SynergyExtrapolation/calcMuscleExcitationsSynX.m b/src/MuscleTendonPersonalization/SynergyExtrapolation/calcMuscleExcitationsSynX.m index 2cf57d5e..eae6b240 100644 --- a/src/MuscleTendonPersonalization/SynergyExtrapolation/calcMuscleExcitationsSynX.m +++ b/src/MuscleTendonPersonalization/SynergyExtrapolation/calcMuscleExcitationsSynX.m @@ -7,11 +7,11 @@ % data: % time - timestamp vectors % timeDelay - electromechanical delay -% emgScalingFactor - EMG scaling factors +% emgScalingFactor - EMG scaling factors % synergyExtrapolation - vector for SynX (including both SynX weights and residual weights) % experimentalData.extrapolationCommands - synergy excitations extracted from measured muscle % excitations for SynX -% experimentalData.residualCommands - synergy excitations extracted from measured muscle excitations +% experimentalData.residualCommands - synergy excitations extracted from measured muscle excitations % for residual muscle excitations % params.missingEmgChannelGroups - Index for expanding SynX EMG to unmeasured muscles % MeasuredMuscIndexExpand - Index for expanding residual EMG to measured muscles @@ -22,9 +22,9 @@ % size(params.taskNames,2) - number of tasks - double % experimentalData.emgData - number of time frames for each trial - double % params.matrixFactorizationMethod - matrix factorization algorithm - 'PCA' or 'NMF' -% params.synergyExtrapolationCategorization - variability of synergy vector weights across trials for +% params.synergyExtrapolationCategorization - variability of synergy vector weights across trials for % SynX reconstruction -% params.residualCategorization - variability of synergy vector weights across trials for +% params.residualCategorization - variability of synergy vector weights across trials for % residual excitation reconstruction % TrialIndex - trial index for each task - cell array (nTask cells) % nOtherDesignVariables - number of other design variables - double @@ -56,7 +56,7 @@ function [muscleExcitations, muscleExcitationsNoTDelay] = ... calcMuscleExcitationsSynX(experimentalData, timeDelay, ... - emgScalingFactor, synergyExtrapolationVariables, params) + emgScalingFactor, synergyExtrapolationVariables, params) experimentalData.emgDataExpanded = round(experimentalData.emgDataExpanded, 4); @@ -79,10 +79,10 @@ % Create EMG splines emgSplines = createEmgSignals(muscleExcitationsNoTDelay, ... experimentalData.emgTime, timeDelay); -% Interpolate Emg +% Interpolate Emg muscleExcitations = evaluateEmgSplines(experimentalData.emgTime, ... emgSplines, timeDelay); -% muscleExcitations = permute(muscleExcitations,[1 3 2]); +% muscleExcitations = permute(muscleExcitations,[1 3 2]); end function UnmeasuredExcitations = getUnmeasuredMuscleExcitations(params, ... @@ -93,33 +93,33 @@ size(params.missingEmgChannelGroups,2), ... length(params.synergyCategorizationOfTrials)); for j = 1:size(params.missingEmgChannelGroups,2) -if strcmpi(params.matrixFactorizationMethod,'PCA') - SynxVars = reshape(extrapolationWeights',(params.numberOfSynergies + 1) * ... - length(params.synergyCategorizationOfTrials), ... - size(params.missingEmgChannelGroups, 2))'; - for i = 1 : length(params.synergyCategorizationOfTrials) - UnmeasuredExcitations(:, j, i) = extrapolationCommands{i} * SynxVars(j, ... - (i - 1) * (params.numberOfSynergies + 1) + 1 : (i - 1) * ... - (params.numberOfSynergies + 1) + params.numberOfSynergies)' + ... - repmat(SynxVars(j, i * (params.numberOfSynergies + 1)), ... - size(emgData, 3) * ... - size(params.synergyCategorizationOfTrials{i}, 2), 1); - end -elseif strcmp(params.matrixFactorizationMethod,'NMF') - SynxVars = reshape(synxVars', params.numberOfSynergies * ... - length(params.synergyCategorizationOfTrials), ... - size(params.missingEmgChannelGroups, 2))'; - for i = 1 : length(params.synergyCategorizationOfTrials) - UnmeasuredExcitations = zeros(emgData * size( ... - params.synergyCategorizationOfTrials{i}, 2), ... - size(params.missingEmgChannelGroups, ... - 2), length(params.synergyCategorizationOfTrials)); - UnmeasuredExcitations(:, j, i) = extrapolationCommands{i} * ... - SynxVars(j, (i - 1) * params.numberOfSynergies + 1 : ... - i * params.numberOfSynergies, :)'; + if strcmpi(params.matrixFactorizationMethod,'PCA') + SynxVars = reshape(extrapolationWeights',(params.numberOfSynergies + 1) * ... + length(params.synergyCategorizationOfTrials), ... + size(params.missingEmgChannelGroups, 2))'; + for i = 1 : length(params.synergyCategorizationOfTrials) + UnmeasuredExcitations(:, j, i) = extrapolationCommands{i} * SynxVars(j, ... + (i - 1) * (params.numberOfSynergies + 1) + 1 : (i - 1) * ... + (params.numberOfSynergies + 1) + params.numberOfSynergies)' + ... + repmat(SynxVars(j, i * (params.numberOfSynergies + 1)), ... + size(emgData, 3) * ... + size(params.synergyCategorizationOfTrials{i}, 2), 1); + end + elseif strcmp(params.matrixFactorizationMethod,'NMF') + SynxVars = reshape(synxVars', params.numberOfSynergies * ... + length(params.synergyCategorizationOfTrials), ... + size(params.missingEmgChannelGroups, 2))'; + for i = 1 : length(params.synergyCategorizationOfTrials) + UnmeasuredExcitations = zeros(emgData * size( ... + params.synergyCategorizationOfTrials{i}, 2), ... + size(params.missingEmgChannelGroups, ... + 2), length(params.synergyCategorizationOfTrials)); + UnmeasuredExcitations(:, j, i) = extrapolationCommands{i} * ... + SynxVars(j, (i - 1) * params.numberOfSynergies + 1 : ... + i * params.numberOfSynergies, :)'; + end end end -end UnmeasuredExcitations = permute(UnmeasuredExcitations, [3 2 1]); end @@ -146,7 +146,7 @@ end concatResidualExcitations = []; for i = 1:size(params.residualCategorizationOfTrials, 2) - concatResidualExcitations = [concatResidualExcitations; + concatResidualExcitations = [concatResidualExcitations; residualExcitations{i}]; end residualExcitations = concatResidualExcitations; @@ -155,11 +155,11 @@ function emgData = updateEmgSignals(missingEmgChannelGroups, emgData, ... unmeasuredEmgSignals) -for i = 1 : size(missingEmgChannelGroups, 2) -for j = 1 : size(missingEmgChannelGroups{i}, 2) - emgData(:, missingEmgChannelGroups{i}(1, j), :) = ... - unmeasuredEmgSignals(:, i, :); -end +for i = 1 : size(missingEmgChannelGroups, 2) + for j = 1 : size(missingEmgChannelGroups{i}, 2) + emgData(:, missingEmgChannelGroups{i}(1, j), :) = ... + unmeasuredEmgSignals(:, i, :); + end end end @@ -168,12 +168,12 @@ %Distribute residual muscle excitations to measured muscles residualExcitationsExpanded = zeros(size(emgData, 3) * ... -size(emgData, 1), size(emgData, 2)); + size(emgData, 1), size(emgData, 2)); for i = 1 : size(currentEmgChannelGroups, 2) -for j = 1 : size(currentEmgChannelGroups{i}, 2) - residualExcitationsExpanded(:, currentEmgChannelGroups{i}(1, j)) = ... - residualExcitations(:, i); -end + for j = 1 : size(currentEmgChannelGroups{i}, 2) + residualExcitationsExpanded(:, currentEmgChannelGroups{i}(1, j)) = ... + residualExcitations(:, i); + end end residualExcitationsExpanded = reshape(residualExcitationsExpanded, ... size(emgData, 3), size(emgData, 1), size(emgData, 2)); @@ -186,13 +186,13 @@ if size(timeDelay, 2) <= 2 for j = 1 : size(emgData, 1) emgSplines{j} = spline(emgTime(j, 1 : 4 : end), ... - emgData(j, :, 1 : 4 : end)'); + emgData(j, :, 1 : 4 : end)'); end else for i = 1 : size(emgData, 2) for j = 1:size(emgData, 1) emgSplines{j, i} = spline(emgTime(j, 1 : 4 : end), ... - emgData(j, i, 1 : 4 : end)); + emgData(j, i, 1 : 4 : end)); end end end diff --git a/src/MuscleTendonPersonalization/SynergyExtrapolation/getLinearInequalityConstraints.m b/src/MuscleTendonPersonalization/SynergyExtrapolation/getLinearInequalityConstraints.m index a73c6776..1814bb71 100644 --- a/src/MuscleTendonPersonalization/SynergyExtrapolation/getLinearInequalityConstraints.m +++ b/src/MuscleTendonPersonalization/SynergyExtrapolation/getLinearInequalityConstraints.m @@ -74,17 +74,32 @@ fullCommands = cat(1, fullCommands, synergyCommands{i}); end fullCommands(:, end+1) = ones(size(fullCommands, 1), 1); +fullEmg = reshape(emgData, size(emgData, 1) * size(emgData, 2), ... + size(emgData, 3)); -totalGroups = length(params.missingEmgChannelGroups) + ... - length(params.currentEmgChannelGroups); -A = zeros(totalGroups * size(fullCommands, 1), numOtherDesignVariables + totalGroups * size(fullCommands, 2)); -for i = 1:totalGroups - A(size(fullCommands, 1) * (i - 1) + 1 : size(fullCommands, 1) * i, ... - numOtherDesignVariables + size(fullCommands, 2) * (i - 1) + 1 : ... - numOtherDesignVariables + size(fullCommands, 2) * i) = fullCommands; +if strcmp(params.matrixFactorizationMethod, 'PCA') + totalGroups = length(params.missingEmgChannelGroups) + ... + length(params.currentEmgChannelGroups); + A = zeros(totalGroups * size(fullCommands, 1), ... + numOtherDesignVariables + totalGroups * size(fullCommands, 2)); + for i = 1:totalGroups + A(size(fullCommands, 1) * (i - 1) + 1 : size(fullCommands, 1) * i, ... + numOtherDesignVariables + size(fullCommands, 2) * (i - 1) + 1 : ... + numOtherDesignVariables + size(fullCommands, 2) * i) = fullCommands; + end + A = cat(1, A, -A); + + b = ones(length(params.missingEmgChannelGroups) * ... + size(fullCommands, 1), 1); + for i = 1:length(params.currentEmgChannelGroups) + b = [b; 1 - fullEmg(:, i)]; + end + b = [b; zeros(length(params.missingEmgChannelGroups) * ... + size(fullCommands, 1), 1)]; + for i = 1:length(params.currentEmgChannelGroups) + b = [b; fullEmg(:, i)]; + end end -A = cat(1, A, -A); -b = ones(totalGroups * size(fullCommands, 1) * 2, 1); end From e18e260ec4315c07e07bff789120db599688a90e Mon Sep 17 00:00:00 2001 From: stwilliams333 Date: Mon, 12 Jun 2023 17:43:56 -0500 Subject: [PATCH 2/5] Fix constraints for multiple trials --- .../SynergyExtrapolation/getLinearInequalityConstraints.m | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/MuscleTendonPersonalization/SynergyExtrapolation/getLinearInequalityConstraints.m b/src/MuscleTendonPersonalization/SynergyExtrapolation/getLinearInequalityConstraints.m index 1814bb71..9b8227b9 100644 --- a/src/MuscleTendonPersonalization/SynergyExtrapolation/getLinearInequalityConstraints.m +++ b/src/MuscleTendonPersonalization/SynergyExtrapolation/getLinearInequalityConstraints.m @@ -78,7 +78,7 @@ size(emgData, 3)); if strcmp(params.matrixFactorizationMethod, 'PCA') - totalGroups = length(params.missingEmgChannelGroups) + ... + totalGroups = length(params.missingEmgChannelGroups) * size(emgData, 2) + ... length(params.currentEmgChannelGroups); A = zeros(totalGroups * size(fullCommands, 1), ... numOtherDesignVariables + totalGroups * size(fullCommands, 2)); @@ -90,12 +90,12 @@ A = cat(1, A, -A); b = ones(length(params.missingEmgChannelGroups) * ... - size(fullCommands, 1), 1); + size(emgData, 2) * size(fullCommands, 1), 1); for i = 1:length(params.currentEmgChannelGroups) b = [b; 1 - fullEmg(:, i)]; end b = [b; zeros(length(params.missingEmgChannelGroups) * ... - size(fullCommands, 1), 1)]; + size(emgData, 2) * size(fullCommands, 1), 1)]; for i = 1:length(params.currentEmgChannelGroups) b = [b; fullEmg(:, i)]; end From dc16caa725b70f55d5b594542f11dc6d18eeb11c Mon Sep 17 00:00:00 2001 From: stwilliams333 Date: Tue, 13 Jun 2023 14:39:20 -0500 Subject: [PATCH 3/5] Fix saving MTP without precal --- .../MuscleTendonPersonalizationTool.m | 2 + .../getLinearInequalityConstraints.m | 98 +++++++++---------- 2 files changed, 51 insertions(+), 49 deletions(-) diff --git a/src/MuscleTendonPersonalization/MuscleTendonPersonalizationTool.m b/src/MuscleTendonPersonalization/MuscleTendonPersonalizationTool.m index 26bf1cf4..a17db3f4 100644 --- a/src/MuscleTendonPersonalization/MuscleTendonPersonalizationTool.m +++ b/src/MuscleTendonPersonalization/MuscleTendonPersonalizationTool.m @@ -38,6 +38,8 @@ function MuscleTendonPersonalizationTool(settingsFileName) optimizedInitialGuess = MuscleTendonLengthInitialization(precalInputs); inputs = updateMtpInitialGuess(inputs, precalInputs, ... optimizedInitialGuess); +else + precalInputs = struct('optimizeIsometricMaxForce', false); end optimizedParams = MuscleTendonPersonalization(inputs, params); diff --git a/src/MuscleTendonPersonalization/SynergyExtrapolation/getLinearInequalityConstraints.m b/src/MuscleTendonPersonalization/SynergyExtrapolation/getLinearInequalityConstraints.m index 9b8227b9..67f50956 100644 --- a/src/MuscleTendonPersonalization/SynergyExtrapolation/getLinearInequalityConstraints.m +++ b/src/MuscleTendonPersonalization/SynergyExtrapolation/getLinearInequalityConstraints.m @@ -49,57 +49,57 @@ function [A,b] = getLinearInequalityConstraints(params, ... numOtherDesignVariables, synergyCommands, emgData) -% if strcmp(params.matrixFactorizationMethod, 'PCA') -% matrixFactorizationFactor = 1; -% elseif strcmp(params.matrixFactorizationMethod, 'NMF') -% matrixFactorizationFactor = 0; -% end -% % Construct the linear constraint, SynX and Residual A Matrix -% aMatrixSynergy = getSynxAMatrix(params, synergyCommands, ... -% emgData, matrixFactorizationFactor); -% aMatrixResidual = getResidualAMatrix(emgData, params.numberOfSynergies, ... -% params.missingEmgChannelGroups, params.residualCategorizationOfTrials, ... -% matrixFactorizationFactor); -% % Concatenate non-SynX design variables with SynX and Residual Matrices -% A = [zeros(2 * size(params.missingEmgChannelGroups, 2) * size(emgData, 1) * ... -% size(emgData, 2), numOtherDesignVariables) [-aMatrixSynergy ... -% aMatrixResidual; aMatrixSynergy aMatrixResidual]]; -% % b (between 0 and 1) -% b = [zeros(size(emgData, 1) * size(emgData, 2) * size( ... -% params.missingEmgChannelGroups, 2), 1); ones(size(emgData, 1) * ... -% size(emgData, 2) * size(params.missingEmgChannelGroups, 2), 1)]; - -fullCommands = synergyCommands{1}; -for i = 2:length(synergyCommands) - fullCommands = cat(1, fullCommands, synergyCommands{i}); -end -fullCommands(:, end+1) = ones(size(fullCommands, 1), 1); -fullEmg = reshape(emgData, size(emgData, 1) * size(emgData, 2), ... - size(emgData, 3)); - if strcmp(params.matrixFactorizationMethod, 'PCA') - totalGroups = length(params.missingEmgChannelGroups) * size(emgData, 2) + ... - length(params.currentEmgChannelGroups); - A = zeros(totalGroups * size(fullCommands, 1), ... - numOtherDesignVariables + totalGroups * size(fullCommands, 2)); - for i = 1:totalGroups - A(size(fullCommands, 1) * (i - 1) + 1 : size(fullCommands, 1) * i, ... - numOtherDesignVariables + size(fullCommands, 2) * (i - 1) + 1 : ... - numOtherDesignVariables + size(fullCommands, 2) * i) = fullCommands; - end - A = cat(1, A, -A); - - b = ones(length(params.missingEmgChannelGroups) * ... - size(emgData, 2) * size(fullCommands, 1), 1); - for i = 1:length(params.currentEmgChannelGroups) - b = [b; 1 - fullEmg(:, i)]; - end - b = [b; zeros(length(params.missingEmgChannelGroups) * ... - size(emgData, 2) * size(fullCommands, 1), 1)]; - for i = 1:length(params.currentEmgChannelGroups) - b = [b; fullEmg(:, i)]; - end + matrixFactorizationFactor = 1; +elseif strcmp(params.matrixFactorizationMethod, 'NMF') + matrixFactorizationFactor = 0; end +% Construct the linear constraint, SynX and Residual A Matrix +aMatrixSynergy = getSynxAMatrix(params, synergyCommands, ... + emgData, matrixFactorizationFactor); +aMatrixResidual = getResidualAMatrix(emgData, params.numberOfSynergies, ... + params.missingEmgChannelGroups, params.residualCategorizationOfTrials, ... + matrixFactorizationFactor); +% Concatenate non-SynX design variables with SynX and Residual Matrices +A = [zeros(2 * size(params.missingEmgChannelGroups, 2) * size(emgData, 1) * ... + size(emgData, 2), numOtherDesignVariables) [-aMatrixSynergy ... + aMatrixResidual; aMatrixSynergy aMatrixResidual]]; +% b (between 0 and 1) +b = [zeros(size(emgData, 1) * size(emgData, 2) * size( ... + params.missingEmgChannelGroups, 2), 1); ones(size(emgData, 1) * ... + size(emgData, 2) * size(params.missingEmgChannelGroups, 2), 1)]; + +% fullCommands = synergyCommands{1}; +% for i = 2:length(synergyCommands) +% fullCommands = cat(1, fullCommands, synergyCommands{i}); +% end +% fullCommands(:, end+1) = ones(size(fullCommands, 1), 1); +% fullEmg = reshape(emgData, size(emgData, 1) * size(emgData, 2), ... +% size(emgData, 3)); +% +% if strcmp(params.matrixFactorizationMethod, 'PCA') +% totalGroups = length(params.missingEmgChannelGroups) * size(emgData, 2) + ... +% length(params.currentEmgChannelGroups); +% A = zeros(totalGroups * size(fullCommands, 1), ... +% numOtherDesignVariables + totalGroups * size(fullCommands, 2)); +% for i = 1:totalGroups +% A(size(fullCommands, 1) * (i - 1) + 1 : size(fullCommands, 1) * i, ... +% numOtherDesignVariables + size(fullCommands, 2) * (i - 1) + 1 : ... +% numOtherDesignVariables + size(fullCommands, 2) * i) = fullCommands; +% end +% A = cat(1, A, -A); +% +% b = ones(length(params.missingEmgChannelGroups) * ... +% size(emgData, 2) * size(fullCommands, 1), 1); +% for i = 1:length(params.currentEmgChannelGroups) +% b = [b; 1 - fullEmg(:, i)]; +% end +% b = [b; zeros(length(params.missingEmgChannelGroups) * ... +% size(emgData, 2) * size(fullCommands, 1), 1)]; +% for i = 1:length(params.currentEmgChannelGroups) +% b = [b; fullEmg(:, i)]; +% end +% end end From 89844cef532845e461d8348f60ea6af65550e90d Mon Sep 17 00:00:00 2001 From: stwilliams333 Date: Tue, 13 Jun 2023 14:50:53 -0500 Subject: [PATCH 4/5] Use excitation penalty --- src/MuscleTendonPersonalization/Optimization/calcMtpCost.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MuscleTendonPersonalization/Optimization/calcMtpCost.m b/src/MuscleTendonPersonalization/Optimization/calcMtpCost.m index 7ec645d8..a1783a61 100644 --- a/src/MuscleTendonPersonalization/Optimization/calcMtpCost.m +++ b/src/MuscleTendonPersonalization/Optimization/calcMtpCost.m @@ -76,7 +76,7 @@ synxModeledValues, modeledValues, experimentalData, costTerm); case "muscle_excitation_penalty" cost = calcMuscleExcitationPenaltyCost( ... - synxModeledValues, experimentalData, costTerm); + synxModeledValues, experimentalData);%, costTerm); otherwise throw(MException('', 'Cost term %s is not valid for MTP', ... costTerm.type)) From 5cac3876cd482683cd42c2fc1a27858ce580af00 Mon Sep 17 00:00:00 2001 From: stwilliams333 Date: Wed, 14 Jun 2023 16:45:15 -0500 Subject: [PATCH 5/5] Excitation cost allowable error and error center from xml --- .../CostTerms/calcMuscleExcitationPenaltyCost.m | 6 ++++-- src/MuscleTendonPersonalization/Optimization/calcMtpCost.m | 2 +- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/MuscleTendonPersonalization/Optimization/CostTerms/calcMuscleExcitationPenaltyCost.m b/src/MuscleTendonPersonalization/Optimization/CostTerms/calcMuscleExcitationPenaltyCost.m index 5de8a3e2..505251ae 100644 --- a/src/MuscleTendonPersonalization/Optimization/CostTerms/calcMuscleExcitationPenaltyCost.m +++ b/src/MuscleTendonPersonalization/Optimization/CostTerms/calcMuscleExcitationPenaltyCost.m @@ -27,13 +27,15 @@ function cost = calcMuscleExcitationPenaltyCost(modeledValues, ... experimentalData, costTerm) +errorCenter = valueOrAlternate(costTerm, "errorCenter", 0.5); +maximumAllowableError = valueOrAlternate(costTerm, "maxAllowableError", 0.25); muscleExcitationsConstraint = modeledValues.muscleExcitationsNoTDelay(: , ... setdiff(1 : size(modeledValues.muscleExcitationsNoTDelay, 2), ... [experimentalData.synergyExtrapolation.missingEmgChannelGroups{:}]), ... experimentalData.numPaddingFrames + 1 : ... size(modeledValues.muscleExcitationsNoTDelay, 3) - ... experimentalData.numPaddingFrames); -cost = 120 * (muscleExcitationsConstraint - 0.5) .^ 8; +cost = 30 / maximumAllowableError * (muscleExcitationsConstraint - errorCenter) .^ 8; cost(isnan(cost))=0; -cost = sum([sqrt(0.1) .* cost ].^ 2, 'all'); +cost = sum((sqrt(0.1) .* cost).^ 2, 'all'); end \ No newline at end of file diff --git a/src/MuscleTendonPersonalization/Optimization/calcMtpCost.m b/src/MuscleTendonPersonalization/Optimization/calcMtpCost.m index a1783a61..7ec645d8 100644 --- a/src/MuscleTendonPersonalization/Optimization/calcMtpCost.m +++ b/src/MuscleTendonPersonalization/Optimization/calcMtpCost.m @@ -76,7 +76,7 @@ synxModeledValues, modeledValues, experimentalData, costTerm); case "muscle_excitation_penalty" cost = calcMuscleExcitationPenaltyCost( ... - synxModeledValues, experimentalData);%, costTerm); + synxModeledValues, experimentalData, costTerm); otherwise throw(MException('', 'Cost term %s is not valid for MTP', ... costTerm.type))