Climate Audit

Mark T.

UC,

I have not looked in-depth at the code that resulted in this: reg_ts 1996×2592×19 786392064 double,
but, if it is possible, break it into 19 separate 1996×2592 blocks. MATLAB is grumpy when it comes to memory allocation. In particular, it requires contiguous memory for a single array. Even on my machine with 4 GB of RAM (well, 3.2 GB usable because of Winders), the largest array I’ve ever seen it create was on the order of 500 MB. In fact, you can run a check:

[uV sV] = memory

and it will return a element in uV called

MaxPossibleArrayBytes: 958464000

Of course, I could play all day and never get a 958 MB array, but it will put you in the correct order of magnitude. 🙂

Mark

Sinan Unur

I have a dumb question. It’s been quite a while since I last used MatLab. Back then, MatLab source code would be contained in basic text files. I have been looking at the *.mat files in http://www.meteo.psu.edu/~mann/supplements/MultiproxyMeans07/data/reconstructions/cps/

Those files seem to be in some kind of binary format. What am I missing? Where can I find the actual source code?

And, what do cps and eiv stand for?

I am not sure if I can be of any help and I am sorry I got stumped at such an elementary level. Any pointers would be much appreciated.

Sinan

Mark T

MAT files are just data files that hold MATLAB variables upon executing a standard save command. They are not really binary (not like a bunch of packed 32-bit floats, for example), but MATLAB’s proprietary format.

64-bit OS may help, but it really depends upon Winders memory management. By going to 64-bit OS, and more memory, there is certainly more available for MATLAB to use. The problem regards contiguous space, which is not often granted by Winders (not that it tries to break memory up, I’ve just never been able to get large arrays even when I have the memory available).

From a coding standpoint, it doesn’t really make sense to store things in large 3D arrays in MATLAB anyway. First, that array needs to be generated prior to running the loop in which it gets used – reg_ts = zeros(1996,2592,19); – else MATLAB must redefine the array (and space it occupies) every time you iterate and add to it. This can immediately be the source of “out of memory” errors. It also slows MATLAB down considerably (keep in mind, the first run of a MATLAB script is interpreted, not compiled, and it does not optimize array growth even on consecutive runs). Also, there’s simply no benefit other than a reduction in the number of arrays that are used unless larger array operations are being performed (I can’t think of a reason for such an idea).

Mark

Steve McIntyre

I’ve uploaded an R script to calculate the two next nearest gridcells for each proxy using great-circle distances as scripts/mann.2008/nextnearest.txt and added two columns next1, next2 to the details data frame recording the jones cell #s of the two closest gridcell with Mann08 instrumental values.

The script scripts/mann.2008/temperature.correlation.txt calculates correlations to Mann instrumental data for a variety of cases: two by two for 1) actual gridcell and 2) two nearest with Mann instrumental values; a) infilled; b)”original”, and collects the rtable and SI versions, for the periods 1850-1995, 1896-1995 and 1850-1949.

The scripts should be turnkey or close to turnkey. I’ve “replicated” the rtable correlations up to .01 for all but 27 series (and up to .001 for most of these). This doesn’t mean that the underlying calculations are “right” in any sense, merely that I’m pretty close to replicating this step. The 27 oddballs are:

id next_infill rtable si
325 chil016 -0.077 0.062 NA
326 chil017 -0.151 -0.037 NA
339 co522 0.084 0.105 NA
395 fisher_1994_agassiz 0.165 0.258 0.5238
396 fisher_1995_agassiz 0.065 -0.065 NA
400 fisher_2002_devon 0.121 -0.116 NA
405 fisher_2002_plateau -0.203 -0.283 NA
492 lutannt12 0.746 0.864 0.8637
512 lutannt30 0.478 0.561 0.5614
526 lutannt43 0.835 0.846 0.8460
549 lutannt64 0.677 0.623 0.6229
576 meese_1994_gisp2 -0.080 0.062 NA
799 nv037 -0.078 0.061 NA
857 or026 0.202 -0.134 NA
910 schweingruber_mxdabd_grid100 0.222 0.363 0.3626
911 schweingruber_mxdabd_grid101 0.465 0.439 0.4387
936 schweingruber_mxdabd_grid20 0.374 0.302 0.3021
946 schweingruber_mxdabd_grid3 0.280 0.215 0.2146
952 schweingruber_mxdabd_grid35 0.233 0.306 0.3059
964 schweingruber_mxdabd_grid48 0.312 0.282 0.2820
982 schweingruber_mxdabd_grid65 0.363 0.428 0.4284
984 schweingruber_mxdabd_grid67 0.528 0.559 0.5593
992 schweingruber_mxdabd_grid78 0.160 0.289 0.2886
1000 schweingruber_mxdabd_grid87 0.310 0.374 0.3742
1004 schweingruber_mxdabd_grid90 0.450 0.397 0.3973
1017 smith_2004_thicknessmm 0.111 0.151 0.1513
1041 stenni_1999_talos 0.161 0.142 0.1424

• Wolfgang Flamme
  Posted Sep 25, 2008 at 5:09 AM | Permalink

  Re: Steve McIntyre (#13),

  http://www.meteo.psu.edu/~mann/supplements/MultiproxyMeans07/readme.txt

  **UPDATED: 23 September 2008 to correct errors with the longitude value in the “schweingruber mxd” grids.

  Steve: I notice that they didn’t acknowledge us. It seems to be a concerted Team policy not to mention me or CA. The PNAS SI isn’t corrected yet. http://www.pnas.org/content/suppl/2008/09/02/0805721105.DCSupplemental/SD1.xls

Steve McIntyre

I’ve sort of located the problem with series 910 and others may be similar. The gridded series coincides exactly with a gridcell (42.5N 107.5W). So there’s a tie for next nearest between 112.5W and 102.5W. The decision rule in propick.m according to Jean S is that tie goes to the Mann series that occurs first in his numbering system, which would be 112.5W. However the correlation to 102.5W is higher and this is the one that is used. Puzzling.

mannid mann lat long jones
1886 1886 1 42.5 -112.5 662
1887 1887 1 42.5 -107.5 663
1888 1888 1 42.5 -102.5 664

• Nick Moon
  Posted Sep 25, 2008 at 4:05 AM | Permalink

  Re: Steve McIntyre (#14),

  The decision rule in propick.m according to Jean S is that tie goes to the Mann series that occurs first in his numbering system, which would be 112.5W. However the correlation to 102.5W is higher and this is the one that is used. Puzzling.

  propick.m looks at the two closest gridcells – by distance and then as pointed out, by order in the data. It works out the correlation to each cell. If correlation is good enough it picks the cell that has the best correlation.

  Here is the code – I’ve snipped out other stuff – and one caveat – this is not a language I actually know. The correlations are done in a loop – variable name of k. It always loops twice, igrid=2, so once with k=1 and once with k=2. That’s the closest two cells or first in the list if tied on distance.

  for k=1:igrid % %numbers of gridpoints near proxy site used to cal. corrcoef
  n2=0; d2=[];

  … snip …

  rr=corrcoef(d2);
  aa(k,1)=(rr(1,2));

  … snip …
  end

  So, at this point, if I follow the code aa(1,1) has the closest cell’s correlation and aa(2,1) has the next.

  if sum(d3)>=igr % at least ‘igr’ correlation of ‘igrid’ pass the CL
  n4=n4+1;

  bb=abs(aa);
  ra=max(bb);

  Take the absolute values of aa, and find the largest value, into ra

  [iii,jjj]=find(bb==ra);
  if sum(iii)>2
  rra(n4)=aa(1,1); %07/20/06
  else
  rra(n4)=aa(iii,jjj);
  end

  d4(n4)=i+1;
  end

  I don’t understand the syntax of that bit. I’d expect iii and jjj to be the index of the element in bb that matches ra. So I can’t see why you would then sum it. However rra(n4) is going to end up set to the larger (positive or negative) in aa().

  There are several different ways in which this whole thing screws up. mxd series aren’t really proxies – so the lat,long are exactly equal to the grid cells. And because the algorithm only considers two close cells, it will use the matching one – and then whichever one it hits next. There being 4 adjacent cells that are exactly equally close. setting igrid to 5 would have been more logical – although frankly why not just see if the proxy matches any temperature anywhere on the planet 🙂

  The other issue, is what happens if there is no data for the instrumental grid cell, i.e the cell is named in the toofew list. This happens mainly close to the poles, so I’d guess it would affect ice core proxies. If there is no instrumental data, it’ll cast around for two cells with data. These might be some way apart and maybe 20 deg of lat away from the proxy location.

  • UC
    Posted Sep 25, 2008 at 6:36 AM | Permalink

    Re: Nick Moon (#17),

    although frankly why not just see if the proxy matches any temperature anywhere on the planet 🙂

    Cut the smiley. That’s what MBH98 did.

    • Nick Moon
      Posted Sep 25, 2008 at 7:23 AM | Permalink

      Re: UC (#21),

      Cut the smiley. That’s what MBH98 did.

      Ahh, but I thought we had moved on.

  • Spence_UK
    Posted Sep 27, 2008 at 3:36 AM | Permalink

    Re: Nick Moon (#17),

    I’ve been too busy to download and try the latest MBH stuff (hopefully get more time soon and have a play with it) but I have been browsing this thread with interest. This caught my eye. I would stress I haven’t seen the full code or context so I may be missing something.

    As noted by Mark T. and others, sum(iii) is a terrible way to test the length of the matrix. The behaviour of this test is going to be dependent on a number of different things; if bb is a row vector, iii will simply be a vector of ones and sum(iii) will indeed yield length (just very inefficiently with lots of extra computation compared to the length function). If bb is a column vector or matrix, iii will contain indices and this test will pass even for some non-vector results.

    Interestingly, in the former case, a two-element vector will fail the test, because the code uses greater than two instead of greater than or equal to two. This would cause the code to fall over so I’m assuming this never happens?

    But the curious thing is this. If the test passes (i.e. sum of iii greater than two), the code appears to just select the first result. But that isn’t what the code seems to do. Selection of the first result would be the code

    rra(n4)=aa(iii(1),jjj(1));

    but the code is actually

    rra(n4)=aa(1,1);

    which is presumably the first entry in the data, which may not be the maximum at all? Have I got this right? That means under certain conditions where you get multiple values through, a high value will be replaced with garbage. If that is most likely to happen for the Luterbacher series, the code may yield results with a correlation value closer to zero than should actually be the case.

    Repeated caveat: this is speculation from the code I’ve seen. I may have the context wrong.

    If I’ve got the right end of the stick, Dr Mann may be interested to know that the max function in MATLAB also provides the index into the array as a second output. This is a neater solution to the apparent intent:

    [maxValue, maxIndex] = max(abs(aa));
    rra(n4) = aa(maxIndex(1));

    The index into maxIndex will work whether maxIndex is a scalar, vector or matrix and is much simpler to follow. Also, Dr Mann, variables called “aa”, “bb” and “iii” really went out of programming fashion a few decades ago. MATLAB pre-processes the m-file into p-code before running it so there is no performance advantage to unreadable short names. And given the size of your arrays, I’m not convinced memory is an excuse, either.

    UC, you may be able to reduce some of the memory by converting those large arrays to single precision. If the code were well written, this shouldn’t have much affect on the results. Of course, I think we already know the answer to that question (as Jean S alludes to above). Unfortunately that means using machines with different maths co-processors (even different versions of the same chip, let alone entirely different architectures) will yield different results as well…

    • Spence_UK
      Posted Sep 27, 2008 at 4:32 AM | Permalink

      Re: Spence_UK (#52),

      Duh… penny drops that aa is already a restricted subset and aa(1,1) may have some special meaning… in which case aa(1,1) may do something not entirely unreasonable. I probably ought to find time to download the code and get some context before commenting 🙂

    • Nick Moon
      Posted Sep 28, 2008 at 5:20 AM | Permalink

      Re: Spence_UK (#53),

      aa contains the correlation coefficients, between the centre of the grid cell and the proxy. Though there are only ever 1 or 2 values in aa. The values maybe + or -. So Mann creates a copy with the absolute values – bb. Then he finds the largest value in that array. Then he uses that value to find the index. Then he uses the index to find the original coefficient in aa. I think there is then a little bit of code to cope with the special case that both coefficients are the same, if so he just takes the first one. It’s hard to believe there isn’t a much cleaner way of doing it. But code like this probably gets tweaked and modified and probably some of this stuff was doing useful work in a previous version.

    • Spence_UK
      Posted Sep 28, 2008 at 11:30 AM | Permalink

      Re: Nick Moon (#54),

      Thanks Nick, I finally downloaded some of the code to figure out more closely what was going on. Not the easiest code to follow.

      One neater and more generic way to do it is in the third blockquote of my #52. But having looked at the code more closely, Mann’s method will work, although for fairly non-obvious reasons. The vector aa is constructed as a 2-element column vector, which means iii will contain the indices into the vector and jjj will just be a vector of ones. iii can therefore only take on three possible values; [1], [2] or [1; 2]. jjj will simply be a vector of 1s as long as iii. Mann’s sum test will therefore always yield the correct answer, so long as igrid is 2. This test won’t work for igrid of 3 or more. Of course the value is hard wired, but the fact that it is parameterised gives it the feel of one of the many tweakable degrees of freedom available to the author – but the code doesn’t support it being used in this way.

      As for the code being tweaked and modified – yes, this is pretty evident in the comments. One convenience of MATLAB is the ability to quickly and easily try out little tweaks and changes. And it appears Mann’s passion for little tweaks and changes haven’t abated any in recent years, no matter how much they can get him into hot water.

      Here’s an interesting observation. On one pass around the larger outer loop, it could be set (and populated) to length 2. But say on a later pass, only one valid record may be found. This would still invoke the piece of code you list (since it requires a minimum of igr, 1 record, to be found). But then you get an inconvenient side effect – you’ve still got the correlation coefficient hanging around from the last time around the loop that you found 2 cases from.

      This may have no consequence whatsoever in the actual application, but nevertheless, this sort of coding is just asking for mistakes, and of course underscores the need for the original code to make any kind of replication. Anyhow, I’ve downloaded most of the files now, just need to generate clidata as Steve notes above and then I’ll be able to get propick running.

• Jean S
  Posted Sep 25, 2008 at 1:45 PM | Permalink

  Re: Steve McIntyre (#14),

  The decision rule in propick.m according to Jean S is that tie goes to the Mann series that occurs first in his numbering system, which would be 112.5W. However the correlation to 102.5W is higher and this is the one that is used. Puzzling.

  I think I figured this out … a sort of. It seems to be a rounding/floating point issue! Variable ‘d1’ has the distances calculated. I checked for the proxy #910, and the three smallest distances are:
  0
  410.1783
  410.1783
  and these correspond to:
  -107.5000 42.5000
  -102.5000 42.5000
  -112.5000 42.5000
  Now the difference 2nd smallest-3rd smallest is: -9.6634e-013 !!!
  So the rule is as I (and Nick Moon) explained, but it happens due to floating point issues that according
  to Mannian code 112.5W and 102.5W are not equally close to 107.5W …

  I’m 100% sure that there is an incredible amount of pointless calculation in this.

  Yes, that’s true. My favorite line (from ‘propick.m’) so far:

  itrdbmatrix=itrdbmatrix;

UC

OK, let’s test your memory. Type

reg_ts=NaN([1996,2592,19]);
saveregts=reg_ts(:,1:1728,:);

(from gridboxcps.m )

Who’s Matlab survives through these two commands ?

• Mark T.
  Posted Sep 25, 2008 at 9:22 AM | Permalink

  Re: UC (#20),

  Who’s Matlab survives through these two commands ?

  Without a doubt the two largest arrays I’ve ever had in a MATLAB session.

  Re: Nick Moon (#17),

  I don’t understand the syntax of that bit. I’d expect iii and jjj to be the index of the element in bb that matches ra. So I can’t see why you would then sum it.

  Yet another case of Mann not understanding MATLAB, IMO. If he wants to test whether there are more than two instances of bb==ra, all he needs to do is use length(iii) rather than sum(iii).

  Mark

Mark T.

Do I get some quatloos for getting it to run without crashing? I’m guessing, btw, that I can’t do a whole lot with MATLAB once I’ve scarfed up that much memory for two arrays. I probably can’t do much with Winders, either (I’m not working very hard at the moment, so it matters not, hehe). 🙂

Mark

Patrick M.

Based on the english in some of the readme files, I’m going to guess that they were written by somebody whose native language is not english.

I wonder how much of the code was written by Mann?

Steve: Which individual wrote the code is irrelevant for our purposes here, just as much as it’s irrelevant to a third party which Pfizer scientist developed a patent. From the directory labels, Zhang probably wrote some of the code; Rutherford probably wrote others, but unless you connect this to something in the output, forget about this line of analysis please.

Jean S

Oh, I forgot. It would be interesting to know if there are any proxies that are “significant/insignificant by chance”. That is, if there are any proxies whose selection depends from the floating point issues in Matlab…

Steve McIntyre

#33. Do you mean that they are insignificant with a “correct” pick two rule under Mannian order, but significant with an “incorrect” pick two?

Jean S

#35: Both ways. Those whose picking depends on how you order the equal distance cells. BTW, the ‘propick.m’ code is clearly designed to implement “pick m” rules. Would be interesting to know why they chose 2, not 1 or 3?

#34: Yes. The correlations are calculated in ‘propickl.m’. It inputs the “low frequency” splits ‘clidatal’ (instrumental) and ‘itrdbmatrixl’ (proxies). ‘clidatal’ is supposedly calculated in ‘clihybrid.m’ (which seems to be missing; at least I can’t locate it) and ‘itrdbmatrixl’ in calculated in ‘proxyhybrid.m’. Here’s the relevant filtering part (additionally there seems to be a lot of “normalization” going on, which is so typical to Mann):

[lowa,lowb]=butter(5,lowlim,’low’);

[snip]

%%%% Apply the zero-phase butterworth filter but only to the proxies
highf=zeros(nrows,ncols);
lowf=zeros(nrows,ncols);
lowf=lowf./lowf;
highf=highf./highf;
for k=1:ncols
clear y y2
data=stdkeepers(start(k):stop(k),k);
% lowf(start(k):stop(k),k)=lowpass(data,0.05,2,2);
frequency = 0.05;
[smoothed0,icb,ice,mse0] = lowpassmin(data,frequency);
lowf(start(k):stop(k),k)= smoothed0;
% lowf(start(k):stop(k),k)=filtfilt(lowa,lowb,data);
highf(start(k):stop(k),k)=data-lowf(start(k):stop(k),k);
end

So before choosing this 2nd latest Mannian filtering (choose the “best” “end condition”; ‘lowpassmin.m’) at least Mannian filtering version 2004 (with “minimum roughness”; ‘lowpass(data,0.05,2,2)’) and stardard Butterworth in Matlab (‘butter’ + ‘filtfilt’) were used. How come he needs to test so many things if his results are so robust?

• Rusty Scott
  Posted Sep 26, 2008 at 1:31 PM | Permalink

  Re: Jean S (#37),

  highf=zeros(nrows,ncols);
  lowf=zeros(nrows,ncols);
  lowf=lowf./lowf;
  highf=highf./highf;

  I found this bit of code confusing because in general I don’t like to divide numbers by zero, even if I’m dividing zero by zero.
  This is apparently Mannian long-hand for
  lowf = highf = NaN(nrows,ncols);
  which would be much easier to read. I can see why it is so hard to follow what Mann does in his code because he changes direction in the code so many times and leaves bits and pieces of old methods hanging around in the code.

  • RomanM
    Posted Sep 26, 2008 at 5:00 PM | Permalink

    Re: Rusty Scott (#48),

    This is apparently Mannian long-hand for lowf = highf = NaN(nrows,ncols);

    Curiously enough, I have an older version of matlab and the NaN function doesn’t seem to work on my system. Maybe that piece of code was written on an earlier version of Matlab.

Steve McIntyre

#18. Although they said it was updated on Sept 23 in this statement, elsewhere they say Sep 25 which seems more likely. Not that it matters particularly, but this may be helpful for those readers who get cross at me when I do not take everything at face value

Jeff Id

Correlation is ‘allegedly’ independent of the multiplier as well. I don’t recall the specific phrase but I believe Mann stated something to that effect. Therefore the offset and final result cannot be calculated without both.

Willis Eschenbach

Jeff Id, if RegEM is some flavor of linear fitting, the final result will likely be (or be able to be recast) in the form:

“temperature” = proxy1*m1 + proxy2*m2 + proxy3*m3 + …. + proxyn*mn +b

So we’d end up with 1210 coefficients, the majority of which will likely be ~ 0 …

w.

Mark T.

Can’t do the double equal thing in MATLAB like you can in C, btw. Mann just isn’t as familiar with the language as he could (or should) be. You’re right, it is confusing at times.

Mark

Spence_UK

OK, got doannsst.m and propick.m running okay. Just needed to redirect the loads and saves.

I think the coding flaw I outlined in #55 seems moot, since it seems to always pick up two locations (i.e. sum(d3) appears always 2, never 0 or 1, at least in the first half of propick). Buggy software along execution paths that are never tested… I probably shouldn’t be surprised!

Steve McIntyre

The starting point of Mann’s program is the infilling of proxies. This might be a good place to watch how the RegEM works as well.

• UC
  Posted Mar 16, 2009 at 2:53 PM | Permalink

  Re: Steve McIntyre (#59),

  This might be a good place to watch how the RegEM works as well.

  Got some code from Jean S.. The truth is out: EIV-Files go soon public 🙂

  • Jean S
    Posted Mar 16, 2009 at 3:47 PM | Permalink

    Re: UC (#61),
    here’s a small preview what to expect 😉 Some proxies weighted according to their weight in the (RegEM) EIV-reconstruction (nhnhscrihad). Blue indicates positive weight and red negative. The lower right corner is the actual (low split) reconstruction (for the step) scaled to the final scaling.

    It is rather interesting how some proxies respond with a different sign depending on the century they are used…