[37] | 1 | function dfdx = shadowjacobian(f,x) |
---|
| 2 | % See SDPVAR/jacobian |
---|
| 3 | |
---|
| 4 | % Author Johan Löfberg |
---|
| 5 | % $Id: shadowjacobian.m,v 1.6 2005/05/10 15:04:52 joloef Exp $ |
---|
| 6 | |
---|
| 7 | if nargin==1 |
---|
| 8 | if isa(f,'sdpvar') |
---|
| 9 | x = recover(depends(f)); |
---|
| 10 | else |
---|
| 11 | x = 0; |
---|
| 12 | end |
---|
| 13 | else |
---|
| 14 | if length(getvariables(x))<length(x) |
---|
| 15 | error('x should be a vector of scalar independant variables'); |
---|
| 16 | end |
---|
| 17 | end |
---|
| 18 | |
---|
| 19 | if isa(f,'double') |
---|
| 20 | dfdx = zeros(size(f,1),length(x)); |
---|
| 21 | end |
---|
| 22 | |
---|
| 23 | [n,m]=size(f); |
---|
| 24 | |
---|
| 25 | if m>1 |
---|
| 26 | error('Jacobian only defined for column vectors.') |
---|
| 27 | end |
---|
| 28 | |
---|
| 29 | if n*m==1 |
---|
| 30 | dfdx = scalar_jacobian(f,x); |
---|
| 31 | % Argh, fix this (sorts inside scalar_jacobian |
---|
| 32 | for i = 1:length(x) |
---|
| 33 | var(i,1)=getvariables(x(i)); |
---|
| 34 | end |
---|
| 35 | [i,j]=sort(var); |
---|
| 36 | dfdx = dfdx(1,j); |
---|
| 37 | return |
---|
| 38 | |
---|
| 39 | else |
---|
| 40 | dfdx = []; |
---|
| 41 | AllVars = recover(unique([depends(f) getvariables(x)])); |
---|
| 42 | |
---|
| 43 | for i = 1:length(f) |
---|
| 44 | dfdx = [dfdx;scalar_jacobian(f(i),x,AllVars)]; |
---|
| 45 | end |
---|
| 46 | % Argh, fix this (sorts inside scalar_jacobian |
---|
| 47 | for i = 1:length(x) |
---|
| 48 | var(i,1)=getvariables(x(i)); |
---|
| 49 | end |
---|
| 50 | [i,j]=sort(var); |
---|
| 51 | dfdx = dfdx(:,j); |
---|
| 52 | end |
---|
| 53 | |
---|
| 54 | function [dfdx,dummy] = scalar_jacobian(f,x,AllVars) |
---|
| 55 | |
---|
| 56 | % Aaaaaah, horrible code. FIX! |
---|
| 57 | |
---|
| 58 | if isa(f,'double') |
---|
| 59 | dfdx = zeros(1,length(x)); |
---|
| 60 | return |
---|
| 61 | end |
---|
| 62 | |
---|
| 63 | if nargin==2 |
---|
| 64 | AllVars = recover(uniquestripped([depends(f) getvariables(x)])); |
---|
| 65 | end |
---|
| 66 | |
---|
| 67 | exponent_p = exponents(f,AllVars); |
---|
| 68 | |
---|
| 69 | coefficients = getbase(f); |
---|
| 70 | coefficients = coefficients(2:end); |
---|
| 71 | coefficients = coefficients(:); |
---|
| 72 | if nnz(exponent_p(1,:))==0 |
---|
| 73 | exponent_p=exponent_p(2:end,:); |
---|
| 74 | end |
---|
| 75 | |
---|
| 76 | x_variables = getvariables(x); |
---|
| 77 | AllVars_variables = getvariables(AllVars); |
---|
| 78 | %AllDeriv = []; |
---|
| 79 | AllDeriv2 = []; |
---|
| 80 | for k = 1:length(x) |
---|
| 81 | wrt = find(ismembc(AllVars_variables,x_variables(k))); |
---|
| 82 | deriv = exponent_p; |
---|
| 83 | deriv(:,wrt) = deriv(:,wrt)-1; |
---|
| 84 | % AllDeriv = [AllDeriv;deriv]; |
---|
| 85 | keep{k} = find(deriv(:,wrt)~=-1); |
---|
| 86 | AllDeriv2 = [AllDeriv2;deriv(keep{k},:)]; |
---|
| 87 | end |
---|
| 88 | |
---|
| 89 | if size(AllDeriv2,1)==0 |
---|
| 90 | dummy = 1; |
---|
| 91 | else |
---|
| 92 | dummy = recovermonoms(AllDeriv2,AllVars); |
---|
| 93 | end |
---|
| 94 | |
---|
| 95 | top = 1; |
---|
| 96 | dfdx=[]; |
---|
| 97 | for k = 1:length(x) |
---|
| 98 | wrt = find(ismembc(AllVars_variables,x_variables(k))); |
---|
| 99 | m = length(keep{k}); |
---|
| 100 | if m>0 |
---|
| 101 | poly = sum((coefficients(keep{k}(:)).*exponent_p(keep{k}(:),wrt)).*dummy(top:top+m-1),1); |
---|
| 102 | top = top + m; |
---|
| 103 | else |
---|
| 104 | poly = 0; |
---|
| 105 | end |
---|
| 106 | % %else |
---|
| 107 | % % poly = 0; |
---|
| 108 | % %end |
---|
| 109 | % tip = tip + LengthUsed{k};%size(exponent_p,1); |
---|
| 110 | % end |
---|
| 111 | %j = find(AllDeriv(top:top+size(exponent_p,1)-1,wrt)~=-1); |
---|
| 112 | %if ~isempty(j) |
---|
| 113 | % poly = sum(coefficients(j).*exponent_p(j,wrt).*dummy(top+j-1)); |
---|
| 114 | %else |
---|
| 115 | % poly = 0; |
---|
| 116 | %end |
---|
| 117 | %top = top + size(exponent_p,1); |
---|
| 118 | |
---|
| 119 | dfdx = [dfdx ; poly]; |
---|
| 120 | end |
---|
| 121 | dfdx = dfdx'; |
---|