[862] | 1 | function H=sidp_transformace(sidp_parameters, rsss_parameters,system,apriori)
|
---|
| 2 | %transformace
|
---|
| 3 | eps=10^-3;
|
---|
| 4 | system.dim=2;
|
---|
| 5 | apriori.eta_tilda0=eps;
|
---|
| 6 | apriori.eta_tilda0_range=apriori.y0_range/(system.sigma+apriori.y0_range);
|
---|
| 7 | apriori.beta_tilda0=eps;
|
---|
| 8 | apriori.beta_tilda0_range=apriori.b0_range^2/(apriori.P0+apriori.b0_range^2);
|
---|
| 9 |
|
---|
| 10 | %pomocne promenne a konstanty
|
---|
| 11 | step=[1 sidp_parameters.n_grid];
|
---|
| 12 |
|
---|
| 13 | zero=zeros(sidp_parameters.horizont,2); %zero(1,:) nema smysl, ale takhle je to pekne
|
---|
| 14 | range=zeros(sidp_parameters.horizont,2);
|
---|
| 15 | eigvec=zeros(sidp_parameters.horizont,2,2);
|
---|
| 16 |
|
---|
| 17 | [H,mi]=init_hyperstate(sidp_parameters,apriori); %vytvori H a mi, nastavi H0 a mi*
|
---|
| 18 |
|
---|
| 19 | for i=1:sidp_parameters.n_pass
|
---|
| 20 |
|
---|
| 21 | for j=1:sidp_parameters.n_iter
|
---|
| 22 | [i j]
|
---|
| 23 | search_region=sidp_parameters.gama^(j-1)*sidp_parameters.lambda^(i-1)*sidp_parameters.search_region_init;
|
---|
| 24 | [H zero range eigvec mi]=update_hyperstate(H,zero, range, eigvec,mi,sidp_parameters.n_grid,1+(i+j>2)); %vygeneruje trajektorie, v zasazenem regionu rovnomerne rozmisti body a prekopiruje nalezena mi*
|
---|
| 25 | %>1
|
---|
| 26 | for k=sidp_parameters.horizont:-1:1 %pozor na meze
|
---|
| 27 | level=sidp_parameters.horizont+1-k;
|
---|
| 28 | k
|
---|
| 29 |
|
---|
| 30 | for l=1:size(H,2)
|
---|
| 31 |
|
---|
| 32 | %if mod(l,1000)==0
|
---|
| 33 | % l
|
---|
| 34 | %end
|
---|
| 35 |
|
---|
| 36 | candidates=H(k,l,2)/H(k,l,1)*generate_candidates(H(k,l,1)/H(k,l,2)*mi(k,l), search_region, sidp_parameters.num_of_candidates,sidp_parameters.generate_candidates_mode);
|
---|
| 37 | %_H(k,l,2)/H(k,l,1)
|
---|
| 38 | %compare candidates
|
---|
| 39 |
|
---|
| 40 | %simple comparing
|
---|
| 41 | realization1=zeros(size(candidates,1),rsss_parameters.n0);
|
---|
| 42 | state(1,:)=H(k,l,:);
|
---|
| 43 |
|
---|
| 44 | for m=1:sidp_parameters.num_of_candidates
|
---|
| 45 |
|
---|
| 46 | best_control=candidates(m);
|
---|
| 47 | for n=1:rsss_parameters.n0
|
---|
| 48 |
|
---|
| 49 | %generate realization
|
---|
| 50 | for o=1:level
|
---|
| 51 |
|
---|
| 52 | eta_0=state(1)/(1+eps-state(1));
|
---|
| 53 | beta_0=sqrt(state(2)/(1+eps-state(2)));
|
---|
| 54 | best_control=-best_control *eta_0/beta_0;
|
---|
| 55 |
|
---|
| 56 | s=randn;
|
---|
| 57 | pom=sqrt(1+best_control^2);
|
---|
| 58 | eta_1=abs(eta_0+beta_0*best_control+pom*s);
|
---|
| 59 | next_state(1)=eta_1/(1+eta_1);
|
---|
| 60 |
|
---|
| 61 | realization1(m,n)=realization1(m,n)+(eta_1)^2;
|
---|
| 62 |
|
---|
| 63 | if o<level
|
---|
| 64 |
|
---|
| 65 | %kalman
|
---|
| 66 | beta_1=pom*beta_0+(best_control*s);
|
---|
| 67 | next_state(2)=beta_1^2/(1+beta_1^2);
|
---|
| 68 |
|
---|
| 69 | state=next_state;
|
---|
| 70 |
|
---|
| 71 | %find in H
|
---|
| 72 | stav=state-zero(k+o,:);
|
---|
| 73 | %for p=1:system.dim
|
---|
| 74 | %stav1(p)=(stav(1)*eigvec(k+o,1,p)+stav(2)*eigvec(k+o,2,p))*range(k+o,p)*(sidp_parameters.n_grid-1)+1;
|
---|
| 75 | %end
|
---|
| 76 | stav1=stav.*range(k+o,:)*(sidp_parameters.n_grid-1)+1;
|
---|
| 77 | index=min(max(round(stav1),1),sidp_parameters.n_grid)-[0 1];
|
---|
| 78 | % plot(H(k+o,:,1),H(k+o,:,2), 'r.' ); hold on; plot(state(1,1),state(1,2), 'b+' );
|
---|
| 79 | index=sum(index.*step);
|
---|
| 80 |
|
---|
| 81 | %find_in_hyperstate5(state, H(k+o,:,:))
|
---|
| 82 |
|
---|
| 83 | %plot(H(k+o,index,1),H(k+o,index,2), 'g+' ); hold off
|
---|
| 84 | %plot(H(k+o,index,1),H(k+o,index,2), 'g+' )
|
---|
| 85 |
|
---|
| 86 | % if imag(index)~=0
|
---|
| 87 | % ;
|
---|
| 88 | % end
|
---|
| 89 | best_control=mi(k+o,index);
|
---|
| 90 | end
|
---|
| 91 |
|
---|
| 92 | end
|
---|
| 93 |
|
---|
| 94 | end
|
---|
| 95 | end
|
---|
| 96 | mean_values=mean(realization1,2);
|
---|
| 97 | [min_val min_index]=min(mean_values);
|
---|
| 98 | best_control=candidates(min_index);
|
---|
| 99 | mi(k,l)=best_control;
|
---|
| 100 | end
|
---|
| 101 | end
|
---|
| 102 |
|
---|
| 103 | vypis(:,:)=H(:,:,1);
|
---|
| 104 | save 'eta_tilda.txt' vypis -ASCII;
|
---|
| 105 | vypis(:,:)=H(:,:,2);
|
---|
| 106 | save 'beta_tilda.txt' vypis -ASCII;
|
---|
| 107 | save 'mi.txt' mi -ASCII;
|
---|
| 108 | end
|
---|
| 109 |
|
---|
| 110 | end
|
---|
| 111 |
|
---|
| 112 | end
|
---|