From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 988 invoked by alias); 3 Apr 2014 21:53:12 -0000 Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org Received: (qmail 978 invoked by uid 89); 3 Apr 2014 21:53:11 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=-0.1 required=5.0 tests=AWL,BAYES_00,FREEMAIL_ENVFROM_END_DIGIT,FREEMAIL_FROM,RCVD_IN_DNSWL_NONE,SPF_PASS autolearn=no version=3.3.2 X-HELO: mail-pd0-f172.google.com Received: from mail-pd0-f172.google.com (HELO mail-pd0-f172.google.com) (209.85.192.172) by sourceware.org (qpsmtpd/0.93/v0.84-503-g423c35a) with (AES128-SHA encrypted) ESMTPS; Thu, 03 Apr 2014 21:53:08 +0000 Received: by mail-pd0-f172.google.com with SMTP id p10so2395247pdj.3 for ; Thu, 03 Apr 2014 14:53:06 -0700 (PDT) X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20130820; h=content-type:mime-version:subject:from:in-reply-to:date:cc :content-transfer-encoding:message-id:references:to; bh=ri+DMLDuAZHVST/MiQxq7VILS0TuH+FBqv4FU01jhvg=; b=Lt09L+qtWyFlHy01M5nZM0gtVTTsbnf9svnaWmlq1Goiq7N6zGs/o66a0FOWy+moEi 3j0tYDNiBELs7uShWCuahKl1PiNZvVKBtPRchFbV2IkjEovlj4swp+ZIqSqESYsjrnKQ v7AuopR8NBSKyuqnBqbxwKyWN0uv5MW4P+fKaIBAsQqWpTytTzWnNPZwDS4Gfm+czpQv 6TbtZqtpsNQkb235OiToZCu8Z1eA7PM2sQe2Z91jsXubhAfs9UzhpJKt2r4ClWc8ZLwH 8JfSiDd7w6gvALOoyuZyU9sW6bgsd2ty+80uQPg5Ta8ASvc/Bw+k34jSm+sBtUD58sCv UscA== X-Received: by 10.68.217.234 with SMTP id pb10mr10425904pbc.142.1396561986817; Thu, 03 Apr 2014 14:53:06 -0700 (PDT) Received: from dhcp-128-189-73-105.ubcsecure.wireless.ubc.ca (dhcp-128-189-73-105.ubcsecure.wireless.ubc.ca. [128.189.73.105]) by mx.google.com with ESMTPSA id yz5sm30576611pac.9.2014.04.03.14.53.06 for (version=TLSv1 cipher=ECDHE-RSA-RC4-SHA bits=128/128); Thu, 03 Apr 2014 14:53:06 -0700 (PDT) Content-Type: text/plain; charset=windows-1252 Mime-Version: 1.0 (Mac OS X Mail 7.2 \(1874\)) Subject: Re: Steffen interpolation From: =?windows-1252?Q?Jean-Fran=E7ois_Caron?= In-Reply-To: <533DD7E5.7090900@colorado.edu> Date: Thu, 03 Apr 2014 21:53:00 -0000 Cc: "gsl-discuss@sourceware.org" Content-Transfer-Encoding: quoted-printable Message-Id: <6878829B-EB42-4C43-BC6D-4E7BB08FC833@phas.ubc.ca> References: <57ABFACA-CA7C-439D-9695-F136F0142156@phas.ubc.ca> <8EA5CA7B-E4C4-48A0-A9EB-BA77F3F1AC11@phas.ubc.ca> <532B23E3.7050700@colorado.edu> <532B2ACF.1060608@colorado.edu> <532B2D8E.8090902@colorado.edu> <532B33B2.6060903@colorado.edu> ,<39D96706-0826-4235-9634-18EE4B97EDB4@phas.ubc.ca> <93B0BFDB4CD56A40BBAE7FB8D8984BF701211E565EBB@EXC4.ad.colorado.edu> <0065F39F-DA69-4394-BA53-E9E9CE163583@phas.ubc.ca> <5339D20F.9030609@colorado.edu> <5339EB15.50803@colorado.edu> <958F2963-D3D1-4BC0-8081-EB56644E3AEE@phas.ubc.ca> <533DD7E5.7090900@colorado.edu> To: Patrick Alken X-SW-Source: 2014-q2/txt/msg00001.txt.bz2 Looks good to me. Jean-Fran=E7ois On Apr 3, 2014, at 14:51 , Patrick Alken wrote: > Ok I used your new text and modified it slightly to say that the method u= ses piecewise cubic polynomials in each interval: >=20 > ---- > Steffen's method guarantees the monotonicity of the interpolating function > between the given data points. Therefore, minima and maxima can only occur > exactly at the data points, and there can never be spurious oscillations > between data points. The interpolated function is piecewise cubic > in each interval. The resulting curve and its first derivative > are guaranteed to be continuous, but the second derivative may be > discontinuous. > ---- >=20 > Does this look ok? >=20 > I added your name to test.c >=20 > Patrick >=20 > On 04/03/2014 01:58 PM, Jean-Fran=E7ois Caron wrote: >> Hi Patrick, yes feel free to change the example dataset. I used it beca= use it=92s the same as I put into the test.c code, and other interpolation = methods used randomly-generated data. >>=20 >> For the description in the docs, I might recommend a different wording: >>=20 >> @deffn {Interpolation Type} gsl_interp_steffen >> Steffen=92s method guarantees the monoticity of the interpolating functi= on >> between the given data points. Thus minima and maxima can only occur >> exactly at the data points, and there can never be spurious oscillations= between data points. >> The interpolated function and its first derivative are guaranteed to be = continuous, >> but the second derivative may be discontinuous. >> @end deffn >>=20 >> Thanks for supporting my work! I=92m very excited to be officially cont= ributing to an open-source project. Could you check the copyright & attrib= ution parts of the code files that I modified? I=92m not sure what is corr= ect, but I see author=92s names and dates. I added mine to the steffen.c, = but should I add it also to test.c and the others? >>=20 >> Jean-Fran=E7ois >>=20 >> On Mar 31, 2014, at 15:24 , Patrick Alken w= rote: >>=20 >>> I couldn't reproduce the figure in Steffen's paper, so I found another = dataset which nicely illustrates oscillation issues with Akima: >>>=20 >>> J. M. Hyman, Accurate Monotonicity preserving cubic interpolation, >>> SIAM J. Sci. Stat. Comput. 4, 4, 1983. >>>=20 >>> The dataset is simpler than your randomly generated plot and I think it= s a little easier to compare the different methods. >>>=20 >>> I added an example program and a figure to the manual (in the steffen b= ranch). >>>=20 >>> I am hoping to finish everything up and merge into master by the end of= the week. >>>=20 >>> Thanks again, >>> Patrick >>>=20 >>>=20 >>> On 03/31/2014 02:37 PM, Patrick Alken wrote: >>>> Ok I made a new branch 'steffen' in the GSL repository with your latest >>>> changes, thanks for all your work on this. I still want to update the >>>> docs a little and do some more testing on my own before merging it into >>>> master. I made a blurb about gsl_interp_steffen in the docs: >>>>=20 >>>> ---- >>>> @deffn {Interpolation Type} gsl_interp_steffen >>>> Steffen's method for monotonic interpolation (not allowing minima or >>>> maxima >>>> to occur between adjacent data points). The resulting curve is >>>> piecewise cubic on each interval with the slope at each grid point >>>> chosen to ensure monotonicity and prevent undesired oscillations. The >>>> first-order derivative is everywhere continuous. >>>> @end deffn >>>> ---- >>>>=20 >>>> Can you read this and make sure I haven't said anything inaccurate? Or >>>> let me know any suggestions you think its important to add for the use= rs >>>> benefit to understand what this method does. >>>>=20 >>>> Thanks, >>>> Patrick >>>>=20 >>>> On 03/27/2014 11:17 AM, Jean-Fran=E7ois Caron wrote: >>>>> By the way, my the second test function in interpolation/test.c uses = randomly-generated data points, but actually serves to nicely illustrate th= e difference between major non-linear interpolation methods. See the linke= d graph for a comparison of the interpolation for those data using my imple= mentation of steffen, and the existing GSL akima and cubic spline methods. >>>>>=20 >>>>> https://github.com/jfcaron3/gsl-steffen-devel/blob/steffen/interpolat= ion/compare.pdf (I couldn=92t send a pdf to the mailing list, and I don=92t= know how to view a pdf on github=92s website, but I guess you can just get= the image when you clone the repo.) >>>>>=20 >>>>> While the cubic spline and akima methods preserve continuity of the s= econd derivatives, they are not monotonic and can have oscillations that ar= e often undesireable. The steffen method sacrifices continuity of the seco= nd derivative (but maintains it for the first) in order to maintain monotic= ity, which also eliminates weird oscillations. In Steffen=92s paper, there= is also an example graph where the akima method is unstable (a very small = change in one data point makes a large change in the interpolated function)= , while the steffen method is stable by construction. >>>>>=20 >>>>> Jean-Fran=E7ois >>>>>=20 >>>>>> On Mar 27, 2014, at 01:10 , Patrick Alken wrote: >>>>>>=20 >>>>>>> The code is looking very good - I will try to find time in the next= few days to do some tests and import it into GSL >>>>>>>=20 >>>>>>> Thanks >>>>>>> Patrick >>>>>>> ________________________________________ >>>>>>> From: gsl-discuss-owner@sourceware.org [gsl-discuss-owner@sourcewar= e.org] On Behalf Of Jean-Fran=E7ois Caron [jfcaron@phas.ubc.ca] >>>>>>> Sent: Wednesday, March 26, 2014 7:10 PM >>>>>>> To: gsl-discuss@sourceware.org >>>>>>> Subject: Re: Compiling & Testing New Interpolation Type >>>>>>>=20 >>>>>>> I have now fixed the problems with the tests and added a more robus= t test with lots of data points. I am effectively ready to give a pull req= uest from my github repo. Let me know what I need to do to facilitate this. >>>>>>>=20 >>>>>>> Jean-Fran=E7ois >>>>>>>=20 >>>>>>> On Mar 25, 2014, at 15:51 , Jean-Fran=E7ois Caron wrote: >>>>>>>=20 >>>>>>>> Git and Github weren=92t as intimidating as I expected. I have a = repo here with the =93steffen=94 branch including my changes: >>>>>>>>=20 >>>>>>>> https://github.com/jfcaron3/gsl-steffen-devel >>>>>>>>=20 >>>>>>>> The Savannah git repo didn=92t include a configure script, and I g= ot my modified GSL+Steffen code to compile by directly modifying interpolat= ion/Makefile AFTER running ./configure, so I=92m not sure how to compile th= e files cloned from my github repo. At least it=92s easier to see the chan= ges now. >>>>>>>>=20 >>>>>>>> Jean-Fran=E7ois >>>>>>>>=20 >>>>>>>> On Mar 25, 2014, at 14:56 , Jean-Fran=E7ois Caron wrote: >>>>>>>>=20 >>>>>>>>> I=92ve improved my initial code greatly. You can find it here: >>>>>>>>>=20 >>>>>>>>> http://bazaar.launchpad.net/~jfcaron/+junk/my_steffen/files >>>>>>>>>=20 >>>>>>>>> You can compile it into GSL by adding in the interpolation/Makefi= le references to =93steffen.c=94, =93steffen.lo=94, and =93steffen.Plo=94 e= xactly where there are currently references to =93akima.*=94. >>>>>>>>>=20 >>>>>>>>> I=92ve tried adding an =93integ=94 method, but I=92m afraid I don= =92t even understand the workings of the integ methods for the existing int= erpolation types. I couldn=92t just copy from the akima.c integ method bec= ause they use a build-in spline calculation function (which I also don=92t = understand). Reading uncommented C code is pretty hard. My test program s= hows that the integration method isn=92t obviously broken, but it fails the= tests I wrote in interpolation/test.c The actual interpolation and deriva= tives seem to work and pass the tests. >>>>>>>>>=20 >>>>>>>>> I=92ve not used github before, so I guess my next move should be = to learn the basics and start using that, since otherwise describing my add= itions & changes are hard to follow. In the meantime, is anyone able to ex= plain how the heck the =93integ=94 methods work? >>>>>>>>>=20 >>>>>>>>> Jean-Fran=E7ois >>>>>>>>>=20 >>>>>>>>> On Mar 20, 2014, at 11:30 , Patrick Alken wrote: >>>>>>>>>=20 >>>>>>>>>> Yes that green curve is rather strange and doesn't seem much bet= ter than the cubic spline. I like simplicity too so lets proceed with impor= ting the steffen code. >>>>>>>>>>=20 >>>>>>>>>> On 03/20/2014 12:18 PM, Jean-Fran=E7ois Caron wrote: >>>>>>>>>>> Definitely an advantage of a) is that it is conceptually simple= . b) is 44 pages while a) is only 7. Even if b) is somehow mathematically= superior, I like the idea of understanding the tools that I am using (and = being able to explain it to my academic supervisor/conference attendees). >>>>>>>>>>>=20 >>>>>>>>>>> The MESA astrophysics library (C++ unfortunately) actually incl= udes both types, and has a little graph to show differences: >>>>>>>>>>> http://mesa.sourceforge.net/interp_1D.html >>>>>>>>>>>=20 >>>>>>>>>>> Actually their graph is confusing, blue is supposed to be a), g= reen b), but the green curve isn=92t piece-wise monotonic between the data = points. I=92m starting to think maybe Stetten and Huynh mean different thi= ngs when they say =93monotonic=94. I=92ll try to read Huynh=92s paper to s= ee if they define what they are trying to do. Steffen is pretty clear abou= t his technique being a for an interpolating function that is monotonic bet= ween data points - i.e. the interpolating function doesn=92t change sign be= tween data points, and extrema can only occur at said data points. >>>>>>>>>>>=20 >>>>>>>>>>> Jean-Fran=E7ois >>>>>>>>>>>=20 >>>>>>>>>>> On Mar 20, 2014, at 11:03 , Patrick Alken wrote: >>>>>>>>>>>=20 >>>>>>>>>>>> I see question 1) is answered by section 4 of Steffen's paper = - the method works on all data sets, and preserves monotonicity in each int= erval, which is nice. They also state that method (c) has some serious draw= backs. >>>>>>>>>>>>=20 >>>>>>>>>>>> Unfortunately paper (b) doesn't reference (a) and so its diffi= cult to tell whether (b) offers any advantage over (a) >>>>>>>>>>>>=20 >>>>>>>>>>>> On 03/20/2014 11:52 AM, Patrick Alken wrote: >>>>>>>>>>>>> Hi, I'm moving this discussion over to gsl-discuss which is m= ore suited >>>>>>>>>>>>> for development issues. >>>>>>>>>>>>>=20 >>>>>>>>>>>>> I have 2 naive questions which you may be able to answer sinc= e you've >>>>>>>>>>>>> been working on this code. >>>>>>>>>>>>>=20 >>>>>>>>>>>>> 1) If the Steffen algorithm is applied to non-monotonic data,= will it >>>>>>>>>>>>> still provide a solution or does the method encounter an erro= r? >>>>>>>>>>>>>=20 >>>>>>>>>>>>> 2) Earlier on the GSL list it was mentioned that there are 3 = different >>>>>>>>>>>>> methods for interpolating monotonic data: >>>>>>>>>>>>>=20 >>>>>>>>>>>>> (a) M.Steffen, "A simple method for monotonic interpolation i= n one >>>>>>>>>>>>> dimension", Astron. Astrophys. 239, 443-450 (1990). >>>>>>>>>>>>>=20 >>>>>>>>>>>>> (b) H.T.Huynh, "Accurate Monotone Cubic Interpolation", SIAM = J. Numer. >>>>>>>>>>>>> Anal. 30, 57-100 (1993). >>>>>>>>>>>>>=20 >>>>>>>>>>>>> (c) Fritsch, F. N.; Carlson, R. E., "Monotone Piecewise Cubic >>>>>>>>>>>>> Interpolation", SIAM J. Numer. Anal. 17 (2), 238=96246 (1980). >>>>>>>>>>>>>=20 >>>>>>>>>>>>> I haven't looked at (c) but it seems that (a) and (b) both us= e piecewise >>>>>>>>>>>>> cubic polynomials and preserve monotonicity. Do you happen to= know if >>>>>>>>>>>>> one method is superior to the other? If one method is signifi= cantly >>>>>>>>>>>>> better than the other two it would make more sense to include= that one >>>>>>>>>>>>> in GSL. >>>>>>>>>>>>>=20 >>>>>>>>>>>>> Patrick >>>>>>>>>>>>>=20 >>>>>>>>>>>>> On 03/20/2014 11:37 AM, Jean-Fran=E7ois Caron wrote: >>>>>>>>>>>>>> Yes, I didn=92t bother doing the integration function at the= time because I was having trouble just compiling. I will add the integrat= ion function, and re-write the eval and deriv/deriv2 functions to use Horne= r=92s scheme for the polynomials. I can generate some comparison graphs us= ing fake data like in Steffen=92s paper, that sounds easy enough. >>>>>>>>>>>>>>=20 >>>>>>>>>>>>>> I=92ll look at the interpolation/test.c file and see if I ca= n come up with similar tests. >>>>>>>>>>>>>>=20 >>>>>>>>>>>>>> Thanks for offering to help with the integration into GSL it= self. I don=92t know a lot of the procedures (or even politics sometimes!)= involved. >>>>>>>>>>>>>>=20 >>>>>>>>>>>>>> Jean-Fran=E7ois >>>>>>>>>>>>>>=20 >>>>>>>>>>>>>> On Mar 20, 2014, at 10:22 , Patrick Alken wrote: >>>>>>>>>>>>>>=20 >>>>>>>>>>>>>>> I did notice you talking about 1.6 in your earlier messages= , but assumed it was a typo and you meant 1.16, oops. >>>>>>>>>>>>>>>=20 >>>>>>>>>>>>>>> On 03/20/2014 11:11 AM, Jean-Fran=E7ois Caron wrote: >>>>>>>>>>>>>>>> My original problem was that I wanted to add an interpolat= ion type to GSL. Specifically I want monotonic cubic-splines following the= description in Steffen (1990): http://adsabs.harvard.edu/full/1990A%26A...= 239..443S >>>>>>>>>>>>>>> I took a quick look at your code earlier and it looks prett= y nice. I noticed you commented out the _integ function - is this something= you could add to make it feature complete with the other interpolation typ= es? >>>>>>>>>>>>>>>=20 >>>>>>>>>>>>>>> It is important to add automated tests for this. Can you lo= ok at interpolation/test.c and design similar tests for your new method? Al= so I think it would be nice to add a figure to the manual illustrating the = differences between cubic, akima, and your new steffen method (similar to t= he figures in the Steffen paper). This would help users a lot when trying t= o decide what method to use. Do you happen to have a dataset which shows a = nice contrast like Figs 1, 3 and 8 from that paper? >>>>>>>>>>>>>>>=20 >>>>>>>>>>>>>>> When everything is ready I would be happy to add it to GSL,= as we are already planning to update the interpolation module for the next= release. When I find some time I want to import the 2D interpolation exten= sion discussed previously, and also add Hermite interpolation. >>>>>>>>>>>>>>>=20 >>>>>>>>>>>>>>> It would be easiest for us if you could clone the GSL git r= epository and make your changes there. You could make a new branch called '= steffen' or something and publish it to github, or just send a patch file t= o me, whichever is easiest. >>>>>>>>>>>>>>>=20 >>>>>>>>>>>>>>> Patrick >>>>>>>>>>>>>>>=20 >>>>>>>>>>>>>>> On Mar 19, 2014, at 18:40 , Dave Allured - NOAA Affiliate <= dave.allured@noaa.gov> wrote: >>>>>>>>>>>>>>>>> More data. I tried the same plain build recipe, GSL 1.16= on our test >>>>>>>>>>>>>>>>> machine which is at Mac OS 10.9.3. Got another perfect b= uild, no make >>>>>>>>>>>>>>>>> check errors, no PPC-related issues. Outputs on request,= please be >>>>>>>>>>>>>>>>> specific. >>>>>>>>>>>>>>>>>=20 >>>>>>>>>>>>>>>>> CC=3Dclang >>>>>>>>>>>>>>>>> CFLAGS=3D-g >>>>>>>>>>>>>>>>> ./configure --prefix /Users/dallured/Disk/3rd/gsl/1.16.os= 10.9 >>>>>>>>>>>>>>>>>=20 >>>>>>>>>>>>>>>>> mac27:~/Disk/3rd/gsl/1.16.os10.9 57> sw_vers >>>>>>>>>>>>>>>>> ProductName: Mac OS X >>>>>>>>>>>>>>>>> ProductVersion: 10.9.3 >>>>>>>>>>>>>>>>> BuildVersion: 13D17 >>>>>>>>>>>>>>>>>=20 >>>>>>>>>>>>>>>>> mac27:~/Disk/3rd/gsl/1.16.os10.9/src 36> \ >>>>>>>>>>>>>>>>> ? grep -i '# [a-z]' ../logfiles/make-check.0319a.log | so= rt | uniq -c >>>>>>>>>>>>>>>>> 45 # ERROR: 0 >>>>>>>>>>>>>>>>> 45 # FAIL: 0 >>>>>>>>>>>>>>>>> 42 # PASS: 1 >>>>>>>>>>>>>>>>> 3 # PASS: 2 >>>>>>>>>>>>>>>>> 45 # SKIP: 0 >>>>>>>>>>>>>>>>> 42 # TOTAL: 1 >>>>>>>>>>>>>>>>> 3 # TOTAL: 2 >>>>>>>>>>>>>>>>> 45 # XFAIL: 0 >>>>>>>>>>>>>>>>> 45 # XPASS: 0 >>>>>>>>>>>>>>>>>=20 >>>>>>>>>>>>>>>>> mac27:~/Disk/3rd/gsl/1.16.os10.9 62> \ >>>>>>>>>>>>>>>>> ? grep -c -i ppc logfiles/*319a*log >>>>>>>>>>>>>>>>> logfiles/configure.0319a.os10.9.log:0 >>>>>>>>>>>>>>>>> logfiles/install.0319a.log:0 >>>>>>>>>>>>>>>>> logfiles/make-check.0319a.log:0 >>>>>>>>>>>>>>>>> logfiles/make.0319a.log:0 >>>>>>>>>>>>>>>>>=20 >>>>>>>>>>>>>>>>> mac27:~/Disk/3rd/gsl/1.16.os10.9 65> \ >>>>>>>>>>>>>>>>> ? grep -i ppc src/config.h src/config.log src/config.stat= us >>>>>>>>>>>>>>>>> src/config.h:/* #undef HAVE_GNUPPC_IEEE_INTERFACE */ >>>>>>>>>>>>>>>>> src/config.log:HAVE_GNUPPC_IEEE_INTERFACE=3D'' >>>>>>>>>>>>>>>>> src/config.status:S["HAVE_GNUPPC_IEEE_INTERFACE"]=3D"" >>>>>>>>>>>>>>>>>=20 >>>>>>>>>>>>>>>>> --Dave >>>>>>>>>>>>>>>>>=20 >>>>>>>>>>>>>>>>> On Wed, Mar 19, 2014 at 5:27 PM, Jean-Francois Caron >>>>>>>>>>>>>>>>> wrote: >>>>>>>>>>>>>>>>>> Dave is correct, I am using an "i686" 64-bit x86 mac. F= or some reason >>>>>>>>>>>>>>>>>> it is still looking for the PPC mac header file. The ./= configure >>>>>>>>>>>>>>>>>> stage correctly identifies my system, so it's a bit stra= nge. Also GSL >>>>>>>>>>>>>>>>>> installs without errors when I do it from MacPorts, and = MacPorts >>>>>>>>>>>>>>>>>> doesn't seem to do anything other than ./configure && ma= ke, from my >>>>>>>>>>>>>>>>>> reading of the portfile. >>>>>>>>>>>>>>>>>>=20 >>>>>>>>>>>>>>>>>> When I get back to my Mac, I will look at the NOTES file= to see if >>>>>>>>>>>>>>>>>> anything needs to be done for 10.9. >>>>>>>>>>>>>>>>>>=20 >>>>>>>>>>>>>>>>>> Jean-Fran=E7ois >=20