From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 19294 invoked by alias); 27 Mar 2014 17:17:58 -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 19274 invoked by uid 89); 27 Mar 2014 17:17:57 -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-f174.google.com Received: from mail-pd0-f174.google.com (HELO mail-pd0-f174.google.com) (209.85.192.174) by sourceware.org (qpsmtpd/0.93/v0.84-503-g423c35a) with (AES128-SHA encrypted) ESMTPS; Thu, 27 Mar 2014 17:17:55 +0000 Received: by mail-pd0-f174.google.com with SMTP id y13so3611246pdi.5 for ; Thu, 27 Mar 2014 10:17:53 -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 :content-transfer-encoding:message-id:references:to; bh=ciLOzUzdA7OSjjl7QSnecXiQMLuOIO+y/kEvt0WdMsQ=; b=do5zFehOLdbjjlHnuA2H/q2YY0CRZ1dIOW2Yo/mgthOS2n/3Dfubs8rjWOFfI0ug+t DZ4v3fTZ8VjiFgAldWjhe8jSbHs+F1HfAJKy/rriQtxru8alDiITSrYDeCLwsI95FVGL j2afHtwU1DqADlbJWi38ADlWxUCTglZ2Gu1z27by7GTzGuSOObNfv+fFTeVy4ch35RLj ZB0uD1MNj0rPQolboGIVH91Hrt4OE7lodl+ndbIMeVkPs7NNdOYb5PvMe8C6x+V6U44/ 3Jz0DlrXR31w9vBjb2PDxVERspIfY1z3XB2reUw/VKDJka2JC15VCaFp1gblRzn3oaz7 /Wkw== X-Received: by 10.66.243.131 with SMTP id wy3mr3242317pac.32.1395940673464; Thu, 27 Mar 2014 10:17:53 -0700 (PDT) Received: from dhcp-128-189-72-21.ubcsecure.wireless.ubc.ca (dhcp-128-189-72-21.ubcsecure.wireless.ubc.ca. [128.189.72.21]) by mx.google.com with ESMTPSA id gj9sm11510277pbc.7.2014.03.27.10.17.52 for (version=TLSv1 cipher=ECDHE-RSA-RC4-SHA bits=128/128); Thu, 27 Mar 2014 10:17:52 -0700 (PDT) Content-Type: text/plain; charset=windows-1252 Mime-Version: 1.0 (Mac OS X Mail 7.2 \(1874\)) Subject: Re: Compiling & Testing New Interpolation Type From: =?windows-1252?Q?Jean-Fran=E7ois_Caron?= In-Reply-To: Date: Thu, 27 Mar 2014 17:17:00 -0000 Content-Transfer-Encoding: quoted-printable Message-Id: <0065F39F-DA69-4394-BA53-E9E9CE163583@phas.ubc.ca> References: <57ABFACA-CA7C-439D-9695-F136F0142156@phas.ubc.ca> <5328EF5C.5010300@colorado.edu> <0665FACA-1454-4C7C-80B5-30D08E71A1B7@phas.ubc.ca> <5329F67E.7080203@colorado.edu> <532A124B.4090608@colorado.edu> <532A13D0.2050609@colorado.edu> <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> To: "gsl-discuss@sourceware.org" X-SW-Source: 2014-q1/txt/msg00034.txt.bz2 By the way, my the second test function in interpolation/test.c uses random= ly-generated data points, but actually serves to nicely illustrate the diff= erence between major non-linear interpolation methods. See the linked grap= h for a comparison of the interpolation for those data using my implementat= ion of steffen, and the existing GSL akima and cubic spline methods.=20=20 https://github.com/jfcaron3/gsl-steffen-devel/blob/steffen/interpolation/co= mpare.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 i= mage when you clone the repo.) While the cubic spline and akima methods preserve continuity of the second = derivatives, they are not monotonic and can have oscillations that are ofte= n undesireable. The steffen method sacrifices continuity of the second der= ivative (but maintains it for the first) in order to maintain monoticity, w= hich also eliminates weird oscillations. In Steffen=92s paper, there is al= so 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), whil= e the steffen method is stable by construction.=20=20 Jean-Fran=E7ois >=20 > On Mar 27, 2014, at 01:10 , Patrick Alken wr= ote: >=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@sourceware.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 robust tes= t with lots of data points. I am effectively ready to give a pull request = 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 got my= modified GSL+Steffen code to compile by directly modifying interpolation/M= akefile AFTER running ./configure, so I=92m not sure how to compile the fil= es cloned from my github repo. At least it=92s easier to see the changes n= ow. >>>=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/Makefile re= ferences to =93steffen.c=94, =93steffen.lo=94, and =93steffen.Plo=94 exactl= y 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 interpol= ation types. I couldn=92t just copy from the akima.c integ method because = they use a build-in spline calculation function (which I also don=92t under= stand). Reading uncommented C code is pretty hard. My test program shows = that the integration method isn=92t obviously broken, but it fails the test= s I wrote in interpolation/test.c The actual interpolation and derivatives= seem to work and pass the tests. >>>>=20 >>>> I=92ve not used github before, so I guess my next move should be to le= arn the basics and start using that, since otherwise describing my addition= s & changes are hard to follow. In the meantime, is anyone able to explain= 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 better t= han the cubic spline. I like simplicity too so lets proceed with importing = 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 supe= rior, 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 includes = 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), green = b), but the green curve isn=92t piece-wise monotonic between the data point= s. I=92m starting to think maybe Stetten and Huynh mean different things w= hen they say =93monotonic=94. I=92ll try to read Huynh=92s paper to see if= they define what they are trying to do. Steffen is pretty clear about his= technique being a for an interpolating function that is monotonic between = data points - i.e. the interpolating function doesn=92t change sign between= 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 interval= , which is nice. They also state that method (c) has some serious drawbacks. >>>>>>>=20 >>>>>>> Unfortunately paper (b) doesn't reference (a) and so its difficult = 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 more s= uited >>>>>>>> for development issues. >>>>>>>>=20 >>>>>>>> I have 2 naive questions which you may be able to answer since 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 error? >>>>>>>>=20 >>>>>>>> 2) Earlier on the GSL list it was mentioned that there are 3 diffe= rent >>>>>>>> methods for interpolating monotonic data: >>>>>>>>=20 >>>>>>>> (a) M.Steffen, "A simple method for monotonic interpolation in one >>>>>>>> dimension", Astron. Astrophys. 239, 443-450 (1990). >>>>>>>>=20 >>>>>>>> (b) H.T.Huynh, "Accurate Monotone Cubic Interpolation", SIAM J. Nu= mer. >>>>>>>> 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 use pie= cewise >>>>>>>> cubic polynomials and preserve monotonicity. Do you happen to know= if >>>>>>>> one method is superior to the other? If one method is significantly >>>>>>>> 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 integration f= unction, and re-write the eval and deriv/deriv2 functions to use Horner=92s= scheme for the polynomials. I can generate some comparison graphs using f= ake data like in Steffen=92s paper, that sounds easy enough. >>>>>>>>>=20 >>>>>>>>> I=92ll look at the interpolation/test.c file and see if I can com= e up with similar tests. >>>>>>>>>=20 >>>>>>>>> Thanks for offering to help with the integration into GSL itself.= I don=92t know a lot of the procedures (or even politics sometimes!) invo= lved. >>>>>>>>>=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 interpolation t= ype to GSL. Specifically I want monotonic cubic-splines following the desc= ription in Steffen (1990): http://adsabs.harvard.edu/full/1990A%26A...239..= 443S >>>>>>>>>> I took a quick look at your code earlier and it looks pretty nic= e. I noticed you commented out the _integ function - is this something you = could add to make it feature complete with the other interpolation types? >>>>>>>>>>=20 >>>>>>>>>> It is important to add automated tests for this. Can you look at= interpolation/test.c and design similar tests for your new method? Also I = think it would be nice to add a figure to the manual illustrating the diffe= rences between cubic, akima, and your new steffen method (similar to the fi= gures in the Steffen paper). This would help users a lot when trying to dec= ide 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 w= e are already planning to update the interpolation module for the next rele= ase. When I find some time I want to import the 2D interpolation extension = discussed previously, and also add Hermite interpolation. >>>>>>>>>>=20 >>>>>>>>>> It would be easiest for us if you could clone the GSL git reposi= tory and make your changes there. You could make a new branch called 'steff= en' or something and publish it to github, or just send a patch file to me,= whichever is easiest. >>>>>>>>>>=20 >>>>>>>>>> Patrick >>>>>>>>>>=20 >>>>>>>>>> On Mar 19, 2014, at 18:40 , Dave Allured - NOAA Affiliate wrote: >>>>>>>>>>>> More data. I tried the same plain build recipe, GSL 1.16 on o= ur test >>>>>>>>>>>> machine which is at Mac OS 10.9.3. Got another perfect build,= no make >>>>>>>>>>>> check errors, no PPC-related issues. Outputs on request, plea= se be >>>>>>>>>>>> specific. >>>>>>>>>>>>=20 >>>>>>>>>>>> CC=3Dclang >>>>>>>>>>>> CFLAGS=3D-g >>>>>>>>>>>> ./configure --prefix /Users/dallured/Disk/3rd/gsl/1.16.os10.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 | sort | = 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.status >>>>>>>>>>>> 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. For so= me reason >>>>>>>>>>>>> it is still looking for the PPC mac header file. The ./confi= gure >>>>>>>>>>>>> stage correctly identifies my system, so it's a bit strange. = Also GSL >>>>>>>>>>>>> installs without errors when I do it from MacPorts, and MacPo= rts >>>>>>>>>>>>> doesn't seem to do anything other than ./configure && make, f= rom my >>>>>>>>>>>>> reading of the portfile. >>>>>>>>>>>>>=20 >>>>>>>>>>>>> When I get back to my Mac, I will look at the NOTES file to s= ee if >>>>>>>>>>>>> anything needs to be done for 10.9. >>>>>>>>>>>>>=20 >>>>>>>>>>>>> Jean-Fran=E7ois >>>>>=20 >>>>=20 >>>=20 >>=20 >=20