* Bug 15299 fixed: plot() for polynomials & rationals
[scilab.git] / scilab / modules / graphics / macros / plot.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2004-2006 - INRIA - Fabrice Leray
3 // Copyright (C) 2008 - INRIA - Jean-Baptiste Silvy
4 // Copyright (C) 2012 - 2016 - Scilab Enterprises
5 // Copyright (C) 2018 - 2020 - Samuel GOUGEON
6 //
7 // This file is hereby licensed under the terms of the GNU GPL v2.0,
8 // pursuant to article 5.3.4 of the CeCILL v.2.1.
9 // This file was originally licensed under the terms of the CeCILL v2.1,
10 // and continues to be available under such terms.
11 // For more information, see the COPYING file which you should have received
12 // along with this program.
13
14 function varargout = plot(varargin)
15     // Try to build a new better parser that could manage things like:
16     // plot(x,y,'X',1:10); // where X stands for Xdata (Matlab recognizes
17     //it and treats it well...)
18     // HISTORY:
19     // 2018:
20     //   plot(x, fun) : call to fun() vectorized
21     //   plot(x, list(fun, params)) implemented
22     // 2019:
23     //   logflag implemented.
24     //   y: support to integers added
25
26     [lhs,rhs]=argn(0);
27
28     if ~rhs
29         //LineSpec and PropertySpec examples:
30         t = 0:%pi/20:2*%pi;
31         tt = t';
32         clf("reset");
33         drawlater();
34         subplot(211);
35         plot(tt, sin(tt), "ro-.", tt, cos(tt), "cya+", tt, abs(sin(tt)), "--mo");
36         subplot(212);
37         plot([t ;t],[sin(t); cos(t)],"xdat",[1:2]);
38         drawnow();
39         return;
40     end
41
42
43     CurColor = 0; // current color used if no color specified via LineSpec
44     // nor PropertyName
45
46
47     ListArg = varargin;
48     nextArgin = 1;
49
50     //detect and set the current axes now:
51     // -----------------------------------
52     if type(ListArg(1)) == 9
53         hdle = ListArg(1);
54         if (hdle.type == "Axes")
55             sca(ListArg(1));
56             ListArg(1) = null(); // remove this parameter from the list
57         else
58             msg = _("%s: Argument #%d: Graphic handle(s) of type ""%s"" expected.\n")
59             warning(msprintf(msg, "plot", 1, "Axes"))
60             return;
61         end
62         nextArgin = 2
63     end;
64
65     // Possible log flags
66     // ------------------
67     tmp = ListArg(1);
68     if type(tmp)==10 then
69         if size(tmp,"*")>1
70             msg = _("%s: Argument #%d: log flags: scalar expected.\n")
71             warning(msprintf(msg, "plot", nextArgin))
72             return
73         end
74         tmp = convstr(tmp);
75         tmp2 = strsubst(tmp, "/\s|l|n/", "", "r");
76         if tmp2<>""
77             msg = _("%s: Argument #%d: log flags: wrong value. Ignored.\n")
78             warning(msprintf(msg, "plot", nextArgin))
79             tmp = "nnn"
80         end
81         logflags = part(tmp+"n",1:2);
82         ListArg(1) = null()
83         nextArgin = nextArgin + 1;
84     else
85         logflags = "nn";
86     end
87     nv = size(ListArg)
88
89     argTypes=[];
90     couple=[];
91
92     typeOfPlot = "plot";
93     provided_data = 2;
94
95     for curArgIndex=1:nv
96         tmp = type(ListArg(curArgIndex))
97         if tmp == 16
98             if typeof(ListArg(curArgIndex))=="rational"
99                 tmp = 18
100             end
101         end
102         argTypes(curArgIndex,1) = tmp
103     end
104
105     Ttmp=argTypes;
106
107     for i=1:nv-1
108         acceptedTypes = [];
109         // double, integers, polynomials, rationals (18), macro function or primitive,
110         //    or list(macro|primitive, params) accepted as second argument
111         acceptedTypes = find(Ttmp(i,1)==1 & (or(Ttmp(i+1,1)==[1,2,8,13,130,15,18])))
112         if (acceptedTypes<>[]) then
113             couple=[couple i];
114             Ttmp(i,1)  = 99; // Replace a known type by 99 (no meaning) to count it once only!
115             Ttmp(i+1,1)= 99; // to avoid having (x1,y1,x2,y2) ->couple=[1,2,3]
116             // With this trick, couple=[1,3];
117         end
118     end
119
120     if (couple==[]) // No data couple found
121         // Search for at least a single data , i.e.: plot(y)
122
123         if ((argTypes(1,1)==1 | argTypes(1,1)==8) & ListArg(1)<>[]) then // case plot(SINGLE y,...)
124             couple = 1;
125             provided_data = 1;
126
127             if (modulo(nv-couple,2)<>0) then
128                 P1 = couple+2 // Position of the first PropertyName field
129             else
130                 P1 = couple+1
131             end
132
133         else
134             warning("Error inside input argument : no data");
135             return;
136         end
137
138     else
139
140         // Some test to check wrong inputs
141         //
142         // 1. Test if 2 data couples (first : type==1, second : type=[1,13,130])
143         // are at least separated by 2 indices
144         if (size(couple, "*") > 1 && couple(2:$)-couple(1:$-1)<2)
145             warning("Error inside input argument !");
146             return;
147         end
148
149         // 2. Test if no string couples happen before P1 (see below for P1 definition)
150         for index=1:couple($)
151             acceptedTypes=[];
152             acceptedTypes=find(Ttmp(index,1)==10 & Ttmp(index+1,1)==10)
153
154             if (acceptedTypes<>[]) then
155                 warning("Error inside input argument : String argument is an unknown option.");
156                 return;
157             end
158         end
159
160         if (modulo(nv-(couple($)+1),2)<>0) then
161             P1 = couple($)+3 // Position of the first PropertyName field
162         else
163             P1 = couple($)+2
164         end
165
166     end
167
168     numplot = size(couple,"*");
169
170     xyIndexLineSpec = zeros(numplot,3);
171     // xyIndexLineSpec is a matrix storing the index of x, y and linespec
172     // if one of these indices is 0 => it does not exist
173     // (which is possible for x and linepsec, not for y)
174
175     if (provided_data == 2) then
176
177         for curCouple=1:size(couple,"*")
178             xyIndexLineSpec(curCouple,1:2) = couple(curCouple) +[0,1] // x,y index storage
179
180             if (couple(curCouple)+2 < P1)
181                 if (argTypes(couple(curCouple)+2,1)==10) then // LineSpec treatment
182                     xyIndexLineSpec(curCouple,3) = couple(curCouple)+2;
183                 end
184             end
185         end
186     else
187         // we are in the case where: plot(SINGLE y,... x not specified
188         // or plot(handle,SINGLE y,...
189         xyIndexLineSpec(1,1) = 0; // no x specified
190         xyIndexLineSpec(1,2) = couple;
191
192         if (couple+1 < P1)
193             if (argTypes(couple+1,1)==10) then // LineSpec treatment
194                 xyIndexLineSpec(1,3) = couple+1;
195             end
196         end
197     end
198
199
200
201     // delay the drawing commands
202     // smart drawlater
203     current_figure=gcf();
204     cur_draw_mode = current_figure.immediate_drawing;
205     current_figure.immediate_drawing = "off";
206
207     // check whether this is the first plot for the axes in which we will draw
208     curAxes = gca();
209     // save auto_clear state.
210     OldAutoClear = curAxes.auto_clear;
211
212     isFirstPlot = (curAxes.children == [])
213
214     //Now, we plot the decomposed plots one by one with their own linespec
215     // provided_data = 2 : x and y are provided
216
217     FinalAgreg=[]; // Final Compound containing all the new created plots.
218
219     for i = 1:numplot
220         // Set off auto_clear for allowing multiple graphics entity
221         // will be restored behond
222         if i>1 then
223             curAxes.auto_clear="off";
224         end
225
226         //default values
227         Marker=[];
228         MarkerSize=1;
229         Color=[];
230         LineStyle=1;
231         Line = %F;
232         Marker = %F;
233
234         if (provided_data == 2) then
235
236             // A function (macro or primitive) is given:
237             if (type(ListArg(xyIndexLineSpec(i,2))) == 13 | ..
238                 type(ListArg(xyIndexLineSpec(i,2))) == 130| ..
239                 type(ListArg(xyIndexLineSpec(i,2))) == 15)
240                 //   We need to build the vector or matrix.
241                 firstarg = ListArg(xyIndexLineSpec(i,1));
242                 sizefirstarg = size(firstarg);
243                 secondarg = ListArg(xyIndexLineSpec(i,2));
244                 params = list();
245                 withParams = type(secondarg)==15
246                 if withParams
247                     if size(secondarg)~=2 | and(type(secondarg(1))~=[13 130])
248                         ResetFigureDDM(current_figure, cur_draw_mode);
249                         msg = _("%s: wrong list() specification for the curve #%d.\n")
250                         error(msprintf(msg, "plot", i))
251                     end
252                     buildFunc = secondarg(1)
253                     secondarg(1) = null()
254                     params = secondarg
255                 else
256                     buildFunc = secondarg
257                 end
258
259                 // We test if the function is vectorized:
260                 isvectorized = %t;
261                 try
262                     s1 = min(3,sizefirstarg(1,1))
263                     s2 = min(3,sizefirstarg(1,2))
264                     tmp = buildFunc(firstarg(1:s1,1:s2), params(:))
265                     isvectorized = and(size(tmp)==[s1 s2])  | size(tmp,1)==s1*s2;;
266                 catch
267                     isvectorized = %f;
268                 end
269
270                 // We evaluate ordinates accordingly:
271                 try
272                     if isvectorized
273                         tmp = buildFunc(firstarg, params(:));
274                     else
275                         tmp = [];
276                         for ii = 1:sizefirstarg(1,2)
277                             for jj = 1:sizefirstarg(1,1)
278                                 tmp(jj,ii) = buildFunc(firstarg(jj,ii), params(:));
279                             end
280                         end
281                     end
282                 catch // An error has occurred:
283                     // reset data
284                     ResetFigureDDM(current_figure, cur_draw_mode);
285
286                     // get error info
287                     [err_message, err_number, err_line, err_func] = lasterror(%t);
288
289                     // yield it
290                     if err_func~="", err_func = """"+err_func+"""", end
291                     msg1 = gettext("%s: Error : unable to evaluate input function %s.")
292                     msg2 = gettext("Error %d at line %d of the function: ''%s''")
293                     error(msprintf(msg1 + ascii(10) + msg2, "plot", ..
294                     err_func, err_number, err_line, err_message));
295                 end
296                 // All right: go on plotting:
297                 ListArg(xyIndexLineSpec(i,2)) = tmp;
298                 // if there is another iteration, we will have error message redefining function.
299                 // we need to clear here and not before, because user must see the warning if needed.
300                 clear buildFunc secondarg;
301
302             // Polynomials or rationals
303             // ------------------------
304             elseif type(ListArg(xyIndexLineSpec(i,2)))==2 | ..
305                    typeof(ListArg(xyIndexLineSpec(i,2)))=="rational"
306                 x = ListArg(xyIndexLineSpec(i,1))
307                 p = ListArg(xyIndexLineSpec(i,2))
308                 if isvector(x)
309                     ListArg(xyIndexLineSpec(i,1)) = x(:) // without warning..
310                     ListArg(xyIndexLineSpec(i,2)) = horner(matrix(p,1,-1),x(:))
311                 else
312                     if size(x,2) <> length(p)
313                         ResetFigureDDM(current_figure, cur_draw_mode);
314                         msg = _("%s: Plot #%d: Numbers of x columns and of %s must match.\n")
315                         if type(p)==2
316                             tmp = _("polynomials")
317                         else
318                             tmp = _("rationals")
319                         end
320                         error(msprintf(msg, "plot", i, tmp))
321                     end
322                     tmp = []
323                     for c = 1:size(x,2)
324                         tmp = [tmp horner(p(c), x(:,c))]
325                     end
326                     ListArg(xyIndexLineSpec(i,2)) = tmp
327                     clear tmp
328                 end
329             end
330             [X,Y] = checkXYPair(typeOfPlot,ListArg(xyIndexLineSpec(i,1)),ListArg(xyIndexLineSpec(i,2)),current_figure,cur_draw_mode)
331         else
332             if or(size(ListArg(xyIndexLineSpec(1,2)))==1)  // If this is a vector
333                 if size(ListArg(xyIndexLineSpec(1,2)), "r") == 1 then
334                     X=1:length(ListArg(xyIndexLineSpec(1,2))); // insert a column abcsissa vector of same length,
335                 else
336                     X=(1:length(ListArg(xyIndexLineSpec(1,2))))'; // insert a row abcsissa vector of same length,
337                 end
338             else                                  // if this is a matrix,
339                 X=(1:size(ListArg(xyIndexLineSpec(1,2)),1))'; // insert a row abcsissa vector with same size
340             end
341             // In both cases (matrix/vector), transpose it now so no warning is issued in checkXYPair().
342             [X,Y] = checkXYPair(typeOfPlot,X,ListArg(xyIndexLineSpec(1,2)),current_figure,cur_draw_mode)
343         end
344
345         // Case if 'Xdata', 'Ydata' or 'Zdata' have been set in (PropertyName,Propertyvalue) couples
346         // must be taken into account now
347
348         // P1 is the position of the first PropertyName field.
349         Property = P1;
350
351         while (Property <= nv-1)
352             PropertyName  = ListArg(Property);
353             PropertyValue = ListArg(Property+1);
354
355             // Xdata can ONLY be a vector (cf. Matlab help)
356             PName = getPlotPropertyName(PropertyName,current_figure,cur_draw_mode);
357             if (PName == "xdata")
358
359                 if (type(PropertyValue)<>1 | and(size(PropertyValue)<>1))
360                     warning("Xdata value must be a column or row vector.");
361                     ResetFigureDDM(current_figure, cur_draw_mode);
362                     return;
363                 else
364                     PropertyValue = PropertyValue(:); // force
365                     if or(size(X))==1  // If X is a vector (inevitably a column vector because checkXYPair always returns a column vector)
366                         X = PropertyValue; // X is replaced by PropertyValue
367                         [X,Y] = checkXYPair(typeOfPlot,X,Y,current_figure,cur_draw_mode)
368                     else // X is a matrix
369                         if size(PropertyValue,"*") == size(X,1)
370                             for j=1:size(PropertyValue,"*")
371                                 X(j,:) = PropertyValue(j,1);
372                             end
373                         else
374                             str="plot : incompatible dimensions in input arguments";
375                             warning(str);
376                             ResetFigureDDM(current_figure, cur_draw_mode);
377                         end
378                     end
379                 end
380
381                 // Ydata ONLY be a vector (contrary to what is said by the Matlab help)
382             elseif (PName == "ydata")
383
384                 if (type(PropertyValue)<>1 | and(size(PropertyValue)<>1))
385                     warning("Ydata value must be a column or row vector.");
386                     ResetFigureDDM(current_figure, cur_draw_mode);
387                     return;
388                 else
389                     PropertyValue = PropertyValue(:); // force
390                     if or(size(Y))==1  // If Y is a vector (inevitably a column vector because checkXYPair always returns a column vector)
391                         Y = PropertyValue; // Y is replaced by PropertyValue
392                         [X,Y] = checkXYPair(typeOfPlot,X,Y,current_figure,cur_draw_mode)
393                     else // Y is a matrix
394                         if size(PropertyValue,"*") == size(Y,1)
395                             for j=1:size(PropertyValue,"*")
396                                 Y(j,:) = PropertyValue(j);
397                             end
398                         else
399                             str="plot : incompatible dimensions in input arguments";
400                             warning(str);
401                             ResetFigureDDM(current_figure, cur_draw_mode);
402                         end
403                     end
404
405                 end
406
407                 // Zdata will be treated after plot building
408             end
409
410             Property = Property+2;
411         end
412
413         //Now we have an array xyIndexLineSpec [numplot x 3] containing indices pointing on T for :
414         // - x (<>0 if existing)
415         // - y
416         // - linespec (<>0 if existing)
417         // for each plot passed in argument
418         //       x | y | linespec
419         //       ----------------
420         //plot1   0|i1 |0    <=> plot(y)
421         //plot2  i2|i3 |0    <=> plot(x,y)
422         //plot3  i4|i5 |i6   <=> plot(x,y,LINESPEC)
423         //...
424
425         if (xyIndexLineSpec(i,3)<>0) then // if we have a line spec <=> index <> 0
426             [Color,Line,LineStyle,Marker,MarkerStyle,MarkerSize,fail] = getLineSpec(ListArg(xyIndexLineSpec(i,3)),current_figure,cur_draw_mode);
427         end
428
429         // The plot is made now :
430         err = execstr("plot2d(logflags,X,double(Y))","errcatch","m");
431
432         if err <> 0
433             mprintf("Error %d : in plot2d called by plot",err);
434             ResetFigureDDM(current_figure, cur_draw_mode);
435             return;
436         end
437
438         agreg=gce();  // when using plot2d, we always have an Compound as the current entity
439
440         FinalAgreg = [agreg FinalAgreg];
441
442         if Color==[]
443             DefaultColor = %T;
444         else
445             DefaultColor = %F;
446         end
447
448         for ii=size(agreg.children,"*"):-1:1
449             curPolyline=agreg.children(ii); // we apply linespec to the lines
450
451             // Color treatment : if no color specified by LineSpec nor PropertyName
452             // Set the default color to the curve
453             if DefaultColor == %T
454                 [Color,CurColor] = setDefaultColor(CurColor);
455             end
456
457             if (Marker == %T)
458                 curPolyline.mark_style=MarkerStyle;
459                 curPolyline.mark_mode ="on";
460                 curPolyline.mark_foreground = Color;
461                 curPolyline.mark_style=MarkerStyle;
462                 curPolyline.mark_size=MarkerSize;
463             else
464                 curPolyline.mark_mode ="off"
465             end
466
467             if (Line == %T)
468                 curPolyline.line_mode="on";
469                 curPolyline.foreground = Color;
470                 curPolyline.line_style = LineStyle;
471             else
472                 curPolyline.line_mode="off"
473             end
474
475             if (Line == %F & Marker ==%F) // no linespec nor PropertyName set
476                 curPolyline.line_mode="on";
477                 curPolyline.foreground = Color;
478                 curPolyline.line_style = LineStyle;
479             end
480
481         end
482     end
483
484     //Reset auto_clear Property
485     curAxes.auto_clear = OldAutoClear;
486
487     ///////////////////////////////////
488     //Global Property treatment      //
489     //PropertyName and PropertyValue //
490     ///////////////////////////////////
491
492
493
494     // Those properties will be applied to Agreg children
495     Agreg = glue(FinalAgreg(1:$))
496
497     nbCompound = find(Agreg.children.type=="Compound")
498
499     while (nbCompound<>[])
500         nbCompound=nbCompound(1);
501         unglue(Agreg.children(nbCompound));
502         nbCompound=find(Agreg.children.type=="Compound")
503     end
504
505
506
507     // P1 is the position of the first PropertyName field.
508     Property = P1;
509
510     Curves = Agreg.children
511     //Curves(:,1) = Curves(:,$:-1:1);
512
513     // set mark_size_unit to 'point' for all the curves
514     Curves.mark_size_unit="point";
515
516     while (Property <= nv-1)
517         setPlotProperty(ListArg(Property),ListArg(Property+1),Curves,current_figure,cur_draw_mode)
518
519         Property = Property+2;
520     end
521
522     // force drawing of box like in matlab
523     // for a first plot
524     // unless we are using centered axes
525     // to keep compatibility with Scilab 4
526     if  isFirstPlot & curAxes.x_location <> "origin" & curAxes.y_location <> "origin" then
527         curAxes.box = "on";
528     end
529
530
531     //postponed drawings are done now !
532     // smart drawnow
533     ResetFigureDDM(current_figure, cur_draw_mode)
534
535     if lhs
536         varargout(1) = Curves;
537     end
538 endfunction