function linhet_demo_optimization

prob = coco_prob;
prob = coco_set(prob, 'ode', 'autonomous', false);
prob = coco_set(prob, 'cont', 'NAdapt', 1);

%% first run to find initial fold
prob1 = coco_set(prob, 'coll', 'NTST', 100);

% zero problems
[t0, x0]  = ode45(@(t,x) linode(t, x, 0.98), 1.56164+[0 2*pi/0.98], [1.01958; 0]);
coll_args = {@linode, @linode_dx, @linode_dp, @linode_dt, ...
  @linode_dxdx, @linode_dxdp, @linode_dpdp, @linode_dtdx, ...
  @linode_dtdp, @linode_dtdt, t0, x0, 'om', 0.98};
prob1 = ode_isol2coll(prob1, 'orb1', coll_args{:});

[data1, uidx1] = coco_get_func_data(prob1, 'orb1.coll', 'data', 'uidx');
maps1 = data1.coll_seg.maps;
bc_funcs = {@linode_bc, @linode_bc_du, @linode_bc_dudu};
prob1 = coco_add_func(prob1, 'po1', bc_funcs{:}, data1, 'zero', 'uidx', ...
  uidx1([maps1.x0_idx; maps1.x1_idx; maps1.T_idx; maps1.p_idx]));

[t0, x0]  = ode45(@(t,x) linode(t, x, 0.88), 1.49982+[0 2*pi/0.88], [1.10077; 0]);
coll_args = {@linode, @linode_dx, @linode_dp, @linode_dt, ...
  @linode_dxdx, @linode_dxdp, @linode_dpdp, @linode_dtdx, ...
  @linode_dtdp, @linode_dtdt, t0, x0, 0.88};
prob1 = ode_isol2coll(prob1, 'orb2', coll_args{:});

[data2, uidx2] = coco_get_func_data(prob1, 'orb2.coll', 'data', 'uidx');
maps2 = data2.coll_seg.maps;
prob1 = coco_add_func(prob1, 'po2', bc_funcs{:}, data2, 'zero', 'uidx', ...
  uidx2([maps2.x0_idx; maps2.x1_idx; maps2.T_idx; maps2.p_idx]));

prob1 = coco_add_glue(prob1, 'glue', uidx1(maps1.p_idx), uidx2(maps2.p_idx), -0.1);
prob1 = coco_add_glue(prob1, 'amp', uidx1(maps1.x0_idx(1)), uidx2(maps2.x0_idx(1)), 'ampdiff', 'inactive');

% adjoints
prob1 = adjt_isol2coll(prob1, 'orb1');

[data1, axidx1] = coco_get_adjt_data(prob1, 'orb1.coll', 'data', 'axidx');
opt1 = data1.coll_opt;
prob1 = coco_add_adjt(prob1, 'po1', 'aidx', ...
  axidx1([opt1.x0_idx; opt1.x1_idx; opt1.T_idx; opt1.p_idx]));

prob1 = adjt_isol2coll(prob1, 'orb2');

[data2, axidx2] = coco_get_adjt_data(prob1, 'orb2.coll', 'data', 'axidx');
opt2 = data2.coll_opt;
prob1 = coco_add_adjt(prob1, 'po2', 'aidx', ...
  axidx2([opt2.x0_idx; opt2.x1_idx; opt2.T_idx; opt2.p_idx]));

prob1 = coco_add_adjt(prob1, 'glue', 'aidx', [axidx1(opt1.p_idx); axidx2(opt2.p_idx)]);
prob1 = coco_add_adjt(prob1, 'amp', 'd.ampdiff', 'aidx', [axidx1(opt1.x0_idx(1)); axidx2(opt2.x0_idx(1))]);

cont_pars = {'ampdiff' 'om' 'd.ampdiff' };% 'd.om'};
coco(prob1, 'linode1', [], 1, cont_pars, {[],[0.4 1.4]});

return
%% branch switch from fold to grow nontrivial adjoint
bd1   = coco_bd_read('linode1');
BPlab = coco_bd_labs(bd1, 'BP');

% zero problems
prob2 = ode_BP2coll(prob, '', 'linode1', BPlab(1));

[data, uidx] = coco_get_func_data(prob2, 'coll', 'data', 'uidx');
maps = data.coll_seg.maps;
prob2 = coco_add_func(prob2, 'po', bc_funcs{:}, data, 'zero', 'uidx', ...
  uidx([maps.x0_idx; maps.x1_idx; maps.T0_idx; maps.T_idx]));

prob2 = coco_add_pars(prob2, 'vel', uidx(maps.x0_idx(2)), 'v');

% branch switch data
chart = coco_read_solution('linode1', BPlab(1), 'chart');
cdata = coco_get_chart_data(chart, 'lsol');

% adjoints
prob2 = adjt_BP2coll(prob2, '', 'linode1', BPlab(1));

[chart, lidx] = coco_read_adjoint('po', 'linode1', BPlab(1), 'chart', 'lidx');
[data, axidx] = coco_get_adjt_data(prob2, 'coll', 'data', 'axidx');
opt = data.coll_opt;
prob2   = coco_add_adjt(prob2, 'po', 'aidx', ...
  axidx([opt.x0_idx; opt.x1_idx; opt.T0_idx; opt.T_idx]), ...
  'l0', chart.x, 'tl0', cdata.v(lidx));

[chart, lidx] = coco_read_adjoint('vel', 'linode1', BPlab(1), 'chart', 'lidx');
prob2 = coco_add_adjt(prob2, 'vel', 'd.v', 'aidx', ...
  axidx(opt.x0_idx(2)), 'l0', chart.x, 'tl0', cdata.v(lidx));
dv_int = [chart.x(1) 1];

cont_pars = {'d.v', 'v', 'p'};
coco(prob2, 'linode2', [], cont_pars, {dv_int, [.9 2]});

end

function y = linode(t, x, p)

x1 = x(1,:);
x2 = x(2,:);
om = p(1,:);

y(1,:) = x2;
y(2,:) = -x2-x1+cos(om.*t);

end

function J = linode_dx(t, x, p)  %#ok<INUSD>

J = zeros(2,2,numel(t));
J(1,2,:) = 1;
J(2,1,:) = -1;
J(2,2,:) = -1;

end

function J = linode_dp(t, x, p)  %#ok<INUSL>

om = p(1,:);

J = zeros(2,1,numel(t));
J(2,1,:) = -t.*sin(om.*t);

end

function J = linode_dt(t, x, p)  %#ok<INUSL>

om = p(1,:);

J = zeros(2,numel(t));
J(2,:) = -om.*sin(om.*t);

end

function dJ = linode_dxdx(t, x, p) %#ok<INUSD>
dJ = zeros(2,2,2,numel(t));
end

function dJ = linode_dxdp(t, x, p) %#ok<INUSD>
dJ = zeros(2,2,1,numel(t));
end

function dJ = linode_dpdp(t, x, p)  %#ok<INUSL>

om = p(1,:);

dJ = zeros(2,1,1,numel(t));
dJ(2,1,1,:) = -t.^2.*cos(om.*t);

end

function dJ = linode_dtdx(t, x, p) %#ok<INUSD>
dJ = zeros(2,2,numel(t));
end

function dJ = linode_dtdp(t, x, p)  %#ok<INUSL>

om = p(1,:);

dJ = zeros(2,1,numel(t));
dJ(2,1,:) = -om.*t.*cos(om.*t)-sin(om.*t);

end

function dJ = linode_dtdt(t, x, p)  %#ok<INUSL>

om = p(1,:);

dJ = zeros(2,numel(t));
dJ(2,:) = -om.^2.*cos(om.*t);

end

function [data, y] = linode_bc(prob, data, u) %#ok<INUSL>

x0 = u(1:2);
x1 = u(3:4);
T  = u(5);
om = u(6);

y = [x1(1:2)-x0(1:2); x0(2); om*T-2*pi];% periodicity in R2 x S1 on Poincare section

end

function [data, J] = linode_bc_du(prob, data, u) %#ok<INUSL>

T  = u(5);
om = u(6);

J = [-1 0 1 0 0 0; 0 -1 0 1 0 0; 0 0 0 1 0 0; 0 0 0 0 om T];

end

function [data, dJ] = linode_bc_dudu(prob, data, u) %#ok<INUSD,INUSL>

dJ = zeros(4,6,6);
dJ(4,5,6) = 1;
dJ(4,6,5) = 1;

end

% a different vector field to play with

function y = ode(t, x, p)

x1 = x(1,:);
x2 = x(2,:);
p1 = p(1,:);
p2 = p(2,:);
p3 = p(3,:);

y(1,:) = x1.*x2.*t;
y(2,:) = -p2.^2.*x2-p1.*x1.^2+cos(p3.*t);

end

function J = ode_dx(t, x, p)

x1 = x(1,:);
x2 = x(2,:);
p1 = p(1,:);
p2 = p(2,:);

J = zeros(2,2,numel(t));
J(1,1,:) = x2.*t;
J(1,2,:) = x1.*t;
J(2,1,:) = -2*p1.*x1;
J(2,2,:) = -p2.^2;

end

function J = ode_dp(t, x, p)

x1 = x(1,:);
x2 = x(2,:);
p2 = p(2,:);
p3 = p(3,:);

J = zeros(2,3,numel(t));
J(2,1,:) = -x1.^2;
J(2,2,:) = -2*p2.*x2;
J(2,3,:) = -t.*sin(p3.*t);

end

function J = ode_dt(t, x, p)

x1 = x(1,:);
x2 = x(2,:);
p3 = p(3,:);

J = zeros(2,numel(t));
J(1,:) = x1.*x2;
J(2,:) = -p3.*sin(p3.*t);

end

function dJ = ode_dxdx(t, x, p)  %#ok<INUSL>

p1 = p(1,:);

dJ = zeros(2,2,2,numel(t));
dJ(1,1,2,:) = t;
dJ(1,2,1,:) = t;
dJ(2,1,1,:) = -2*p1;

end

function dJ = ode_dxdp(t, x, p)

x1 = x(1,:);
p2 = p(2,:);

dJ = zeros(2,2,3,numel(t));
dJ(2,1,1,:) = -2*x1;
dJ(2,2,2,:) = -2*p2;

end

function dJ = ode_dpdp(t, x, p)

x2 = x(2,:);
p3 = p(3,:);

dJ = zeros(2,3,3,numel(t));
dJ(2,2,2,:) = -2*x2;
dJ(2,3,3,:) = -t.^2.*cos(p3.*t);

end

function dJ = ode_dtdx(t, x, p) %#ok<INUSD>

x1 = x(1,:);
x2 = x(2,:);

dJ = zeros(2,2,numel(t));
dJ(1,1,:) = x2;
dJ(1,2,:) = x1;

end

function dJ = ode_dtdp(t, x, p)  %#ok<INUSL>

p3 = p(3,:);

dJ = zeros(2,3,numel(t));
dJ(2,3,:) = -sin(p3.*t)-p3.*t.*cos(p3.*t);

end

function dJ = ode_dtdt(t, x, p)  %#ok<INUSL>

p3 = p(3,:);

dJ = zeros(2,numel(t));
dJ(2,:) = -p3.^2.*cos(p3.*t);

end
