* Bug #14582: gettext or _ were applied to broken literal strings
[scilab.git] / scilab / modules / scicos / macros / scicos_scicos / Compute_cic.sci
1 //  Scicos
2 //
3 //  Copyright (C) INRIA - METALAU Project <scicos@inria.fr>
4 //
5 // This program is free software; you can redistribute it and/or modify
6 // it under the terms of the GNU General Public License as published by
7 // the Free Software Foundation; either version 2 of the License, or
8 // (at your option) any later version.
9 //
10 // This program is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 // GNU General Public License for more details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with this program; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 //
19 // See the file ../license.txt
20 //
21 function ok=Compute_cic(method,Nunknowns)
22     global icpr
23     ok=%f
24
25     if icpr==list() then
26         return
27     end
28
29     tolerances=scs_m.props.tol
30     solver=tolerances(6)
31     if (and (solver <> [100 101 102])) then solver=100; end
32     tolerances(6)=solver
33     Atol=tolerances(1)
34     Rtol=tolerances(2)
35     %tcur=0;
36     tf=tolerances(3)
37     state=icpr.state;
38     nx=size(state.x,"r");
39     nx2=round(nx/2);
40     nxModelica=evstr(Nunknowns);
41
42     if nxModelica~=nx2 then
43         msg = _("Your model contains states defined in standard Scicos blocks.\nCurrent initialization interface does not support mixed models.")
44         messagebox(msprintf(msg),"error","modal")
45         return
46     end
47
48     //------------------------------
49     ierr=0;
50     try
51         ierr=execstr("[state,t]=scicosim(state,%tcur,tf,icpr.sim,''start'',tolerances)","errcatch")
52         if ierr<>0 then
53             messagebox(msprintf(_("Initialisation problem in %s "),"Scicosim-start"),"error","modal");
54             return
55         end
56     catch
57         messagebox(msprintf(_("Initialisation problem in %s "),"Scicosim-start"),"error","modal");
58         return
59     end
60     //--------------------------------------------------------------
61     Ida_ierr=0;
62     if method=="Ida (init)" then
63         try
64             // scicos menu, not any more in xcos
65             // setmenu(curwin,'stop')
66             Ida_ierr=execstr("[state,t]=scicosim(state,%tcur,tf,icpr.sim,''run'',tolerances)","errcatch")
67             // unsetmenu(curwin,'stop')
68             if Ida_ierr<>0 then,
69                 messagebox(msprintf(_("Initialisation problem in %s "),"Sundials"),"error","modal");
70             end
71         catch
72             messagebox(msprintf(_("Initialisation problem in %s "),"Sundials"),"error","modal");
73         end
74     end
75     //--------------------------------------------------------------
76     if method=="Kinsol" then
77         try
78             // scicos menu, not any more in xcos
79             // setmenu(curwin,'stop')
80             ierr=execstr("[state2,t]=scicosim(state,%tcur,tf,icpr.sim,''Kinsol'',tolerances)","errcatch")
81             if ierr==0 & (or(isnan(state2.x)) | or(isinf(state2.x))) then
82                 ierr=-1;
83             end
84             if ierr==0 then
85                 state=state2;
86             end
87             // unsetmenu(curwin,'stop')
88             if ierr<>0 then,
89                 messagebox(msprintf(_("Initialisation problem in %s "),"Kinsol"),"error","modal");
90             end
91         catch
92             messagebox(msprintf(_("Initialisation problem in %s "),"Kinsol"),"error","modal");
93         end
94     end
95     //--------------------------------------------------------------
96     if method=="Fsolve" then
97         try
98             x0=state.x(1:nx2);
99             [xres]=fsolve(x0,fsim);
100             ierr=0
101             if or(isnan(xres)) | or(isinf(xres)) then
102                 ierr=-1;
103             end
104             if ierr==0 then
105                 for i=1:nx2, state.x(i)=xres(i);end
106                 fsim(xres);// just to perform an idoit to update outputs in mixed_models
107             end
108         catch
109             messagebox(msprintf(_("Initialisation problem in %s "),"Fsolve"), "error","modal");
110         end
111     end
112     //--------------------------------------------------------------
113     if method=="Optim" then
114         try
115             x0=state.x(1:nx2);
116             [f,xres]=optim(fsumsquare,x0);
117             ierr=0
118             if or(isnan(xres)) | or(isinf(xres)) then
119                 ierr=-1;
120             end
121             if ierr==0 then
122                 for i=1:nx2, state.x(i)=xres(i);end
123                 fsim(xres);// just to perform an idoit to update outputs in mixed_models
124             end
125         catch
126             messagebox(msprintf(_("Initialisation problem in the %s method"),"Optim"), "error","modal");
127         end
128     end
129     //--------------------------------------------------------------
130     if method=="Nelder_Mead" then
131         try
132             x0=state.x(1:nx2);
133             [xmin,fmin,epsilo,ls,fs] = neldermead(rand(nx2,nx2+1),Atol,1,0.5,1.5)
134
135             for i=1:nx2, state.x(i)=xmin(i);end
136         catch
137             messagebox(msprintf(_("Initialisation problem in the %s method")," Nelder_Mead"), "error","modal");
138         end
139     end
140     //--------------------------------------------------------------
141     if method=="Hompack77" then
142         try
143             ierr=execstr("[state,t]=scicosim(state,%tcur,tf,icpr.sim,''hompack77'',tolerances)","errcatch")
144             if ierr<>0 then,
145                 messagebox(msprintf(_("Initialisation problem in the %s method"),"hompack77"), "error","modal");
146             end
147         catch
148             messagebox(msprintf(_("Initialisation problem in the %s method"),"hompack77"), "error","modal");
149         end
150
151     end
152     //--------------------------------------------------------------
153     if method=="Fsolve_Stepping" then
154         x0=state.x(1:nx2);
155         Res0=fsim(x0);
156         Lambda=0;
157         Steps=1000;
158         xres=x0;
159         for i=1:Steps
160             Lambda=i/Steps;
161             [xres]=fsolve(xres,fsim_step);
162             if modulo(i,10)==0 then,
163                 disp("Source-stepping: progress="+string(i/Steps*100)+"%, Error="+string(norm(fsim_step(xres))));
164             end
165         end
166         for i=1:nx2, state.x(i)=xres(i);end
167     end
168     //--------------------------------------------------------------
169     if method=="Sundials_Stepping" then
170     end
171     //--------------------------------------------------------------
172     Err="?"
173     if Ida_ierr==0 & ierr==0 then
174         ss=fsim(state.x(1:nx2));if ss~=[] then Err=string(max(abs(ss)));else Err="0";end//using inf norm
175     end
176
177     try
178         if  Ida_ierr==0 then //cossimdaskr is followed by a cosend in case of error
179             ierr=execstr("[state,t]=scicosim(state,%tcur,tf,icpr.sim,''finish'',tolerances)","errcatch")
180             if ierr<>0 then
181                 messagebox(_("Initialisation problem in the finish phase"), "error","modal");
182                 return
183             end
184         end
185     catch
186         messagebox(_("Initialisation problem in the finish phase"), "error","modal");
187         return
188     end
189     ok=%t;
190 endfunction
191
192 //------------------------------------------------------------
193 function  res=fsim(xin)
194     nx=size(xin,"r");
195     if nx==0 then res=[];return ;end
196     state1=state
197     for i=1:nx, state1.x(i)=xin(i);end
198     ierr=execstr("[state2,t]=scicosim(state1,%tcur,tf,icpr.sim,''linear'',tolerances)","errcatch")
199     res=state2.x(1:nx);
200 endfunction
201
202 function  res=fsim_step(xin)
203     res=fsim(xin)-(1-Lambda)*Res0;
204 endfunction
205
206
207 function  [sumsq,grad,ind]=fsumsquare(xin,ind)
208     nx=size(xin,"r");
209     tolerances=scs_m.props.tol;
210     atol=tolerances(1);
211     rtol=tolerances(2);
212
213     grad=[]
214     if ind==2 | ind==4 | ind==3 then
215         res=fsim(xin);
216         sumsq=0;  for i=1:nx,sumsq=sumsq+res(i)*res(i);end
217     end
218
219     if ind==3 | ind==4 then
220         for j=1:nx
221             xin_p=xin;
222             ewt_j=1/(abs(xin_p(j)*rtol+atol+%eps));
223             delta_j=max(abs(xin_p(j))*%eps,1/ewt_j);
224             xin_p(j)=xin_p(j)+delta_j;
225             res_p=fsim(xin_p);
226             sumsq_p=0;  for i=1:nx,sumsq_p=sumsq_p+(res_p(i)^2-res(i)^2)/delta_j;end
227             grad(j)=sumsq_p;
228         end
229     end
230
231 endfunction
232
233
234 function [xmin,fmin,epsilo,ls,fs] = neldermead(s,epsil,alpha,beta,gama)
235     // saved in the last versions
236 endfunction
237
238
239 function y=fsim2(x)
240     x0=x(1);x1=x(2);x2=x(3);x3=x(4);
241     x4=x(5);x5=x(6);x6=x(7);x7=x(8);
242
243     v0 = -x3;
244     v1 = -x2;
245
246     y(1) = x4+v0;
247     y(2) = x6+v1;
248     y(3) = x7-x6;
249     y(4) = x5-x4;
250     y(5) = 1e-14-x3*x2;
251     y(6) = v0+x2;
252     y(7) = v0+abs(x1);
253     y(8) = v1+abs(x0);
254
255     disp(y');
256 endfunction
257