Yearly Archives: 2008

禁网

向大家宣布我的space目前被禁止国内互联网访问了,我前面就奇怪,这段时间访问量一下子下来了。奥运真是很和谐,封网的同志们也辛苦了。刚才翻了一下以前的帖子,删除了有关杨佳和瓮安的两则评论。刚才跟朋友开玩笑,如果因此惹事,可要记得救我出来,哈哈。

无题

奥运两天后就开始,最高领导人也指示说当前的首要任务是搞好奥运。从个人讲我还是很高兴看到奥运在北京开,既然已经付出了这么多的代价,更是一定要搞好了,更好展示中国的形象。

然而,前不久的四川地震好像一下子被大家遗忘到后面了,或者是那些矛盾都已经安置好了,怕不是这么回事。

今天收到Oxfam America的捐款匹配计划,据说是dollar-for-dollar。想着既然有这匹配计划,不妨咱们也慈善一把,让美国佬也多少出点钱接济咱国家灾后需要钱的人。

其实咱也不宽裕,平时能省就省,咱又不是富豪阶层,甚至连中产都算不上,今天听到一朋友在Ohio买了1万4的二手车(才开了2000km的polo),都为他肉疼了一把。不过捐款这事从来跟有没有钱没太大关系。前面看到新闻说,原三峡某公司的工程师年纪轻轻就过世了,亲朋好友发现他在世的时候捐了好多钱帮过好多人,不由让人感动了一把。可怜的是现在的所谓门户网站上到处是王菲又怀孕了、成龙下车捡垃圾、奥运村布置的如何舒服,以及准备了多少多少避孕套方便运动员在中国搞女人之类的新闻,而这等真正感人的事一般都躲在某个角落,也许是真没几个人关心了,或者更准确的说记者编辑们不关心。

你看,一不小心又小小愤青了一把,就此打住。

使用matlab拟合Gamma分布

南卓铜 ([email protected])

任务描述: 有一批数据,需要将其histogram拟合成Gamma分布,histogram分20个bin,从0到100。结果如下图所示:

hist_dist_final

实现步骤:

1. 装载数据进来,数据文件是data.for.hist.txt,文件结构是一列无header的数据,共43个数据。

>> data=load(‘data.for.hist.txt’)

data是43×1的矢量。

2. 打开distribution fitting tool

>> dfittool

选择Display type为Density(PDF)。点击Data…,弹出Data对话框。

在此对话框内,设置Data为 ‘data’矢量,Censoring和Frequency为 none。设置data set name为dataset,然后点击 Create Data Set。

点击Set Bin Rules,设置Bin width为5。其余不变,OK确认关闭Set Bin Rules。这时,我们将data的histogram设置为每个bin长度为5。在不同的应用里,bin根据具体情况调整。关闭Data对话框。

3. Gamma拟合

点击New Fit…,在New Fit对话框里选择Data为dataset,Distribution为Gamma,点击Apply进行拟合。点击Close关闭New Fit对话框。

此时,看到下图的效果:

image

我们需要修改它的X域,使之匹配我们的histogram的X最大值100。点击Tools菜单下的Axes Limit Control,将X Lower Limit和X Upper Limit更改为 0和100。

image

4. 导出成figure

使用File菜单下的Print to Figure导出成figure。对figure进行必要修改。

>>figure(1)  (注,如果已经打开多个figure窗口,导出的figure序号可能不是1而是其它数字,注意figure窗口标题上的Figure n里的n)
>>legend hide
>>xlabel ‘1/32 degree, mm/day’
>>ylabel ‘Frequency’

此时效果如下,

image

5. 修改figure格式

设置figure为Tools-> edit plot状态,打开View->Property Editor。

将histogram的plot type从line改成Area,配置合适的face color和edge color。

将x轴的范围X Limits设置为0 to 100。

image

这时,histogram压着拟合曲线,选择拟合曲线,右键Cut后,再Paste进来,可以将拟合曲线带到最上面。

6. 修改数据

注意到目前为止,Y轴仍是Density。我们需要将Y坐标修改为Frequency(或者count)。Count=Density x bin width x data count,这里bin width是5,data count是43。因此,我们需要对Ydata乘以5×43=215。

选择拟合曲线,

>>ydata=get(gco,’YData’)
>>ydata_215 = ydata*215
>>set(gco,’YDATA’,ydata_215)

选择histogram,

>>set(gco,’YDATA’,get(gco,’YDATA’)*215)

修改xlim和ylim

>>ylim([0 10])
>>xlim([0 100])

将figure还原为非编辑状态。Tools->Edit Plot,并关闭Property editor。

image

以上步骤中数次修改xlim和ylim,并不是每次都需要,这里仅是为了演示需要。

7. 保存图形

>>print -dpng ‘histogram_dist.png’

在figure对话框的File-> Save As…保存为MATLAB Figure(.fig),以后可以通过File->Open…打开。

结论

本文演示了如何使用Matlab对数据的histogram使用Gamma分布进行拟合。介绍了matlab的常用命令,Distribution Fitting Tool,figure的格式化编辑功能,对figure对像数据的修改等。也介绍了比如如何调整figure里对像的z-order等实用技巧。

data.for.hist.txt

围棋是外星人发明的

聂棋圣已经老了,风光不再。所以仅有的新闻都是在那故作姿态,标新立异以求出位。这不,才有美国宇航局(NASA)的前宇航员说美国政府保留有外星人尸体,这边老聂就说了,围棋是外星人发明的,原因是围棋太难了,太深奥了。我也会点点围棋,不过水平实在不怎么样,我也感觉围棋很深奥,易学难精。然而,这深奥跟发明却是完全两码子事,古人发明出十九乘十九路格子(有说古围棋是15×15路,这本身也是围棋发展的一种证据),拿黑白子来充填形成一种规则,我想不出来有多么的复杂,在漫长的人类历史中,发明这种有趣的规则看起来好像不是特别的困难,因为人的想像力有时候真是很伟大。数学上不也有这么多的假说未得到严谨的证明吗,自然科学里观察到现象却无法解释的例子也比比皆是,但不能证明或者无法解释并不代表假说不能被提出或者现象不能被观测到。围棋的魅力之一就在于简单的规则之下有着非凡的奥秘。

再接下,老聂又说了,“俞斌的确这样说过,但在我看来,恐怕再过3000年,电脑(下围棋)也胜不了人脑。”他太当围棋是一码事了,大概感觉围棋就是天底下最难的事,或者太过维护咱们古老国粹的面子了。俞大侠尽管拿过几个围棋冠军,他自己也承认从围棋水平上讲,不算最顶尖,他也好像编写过围棋程序,不过水平没有高过同类产品。而且我个人怀疑俞大侠的围棋程序恐是跟人合作的结果,他大概也没有时间和能力去独立完成一个商业化的围棋程序。俞大侠对计算机底层的理解怕也有限(无知无畏,才敢断言3000年后大家尸骨无存的世界)。至于老聂,根本不懂计算机,也不会理解计算机的发展和未来可能的潜力,他说的三千年电脑(下围棋)胜不了人脑的断言跟一个白痴发表类似声明一样的缺少权威性。

理由?围棋的计算机实现本质上讲可以通过穷举法完美实现,但这种算法目前来讲不现实,其原因是当前的计算机达不到如此巨大的计算能力。但3000年后呢,开玩笑,46年第一台计算机发明到现在才百年不到,已经发展到目前程度,凭老聂或者俞斌这些非专业人士可以想像3000年后计算能力会发展到什么程度,真是天方夜谈,事实上即使是世界上最顶尖的计算机专家都没有这个预言的能力吧。如果3000年后计算能力不成问题,围棋的实现方案就是小学生计算机课本里学习递归算法的练习题。这还在建立在没有其它优化算法的基础上,一旦人工智能发展到一定程序,计算机具备一定的知识学习和推理功能,围棋的思维方式并没有比其他多数人类智力活动更为特殊和困难——大部分的智力活动目前的计算机都无法模拟,而不仅只有围棋。除了人工智能,围棋算法实现上也是有可能取得突破。

目前看,科技发展迅猛,30年尚不敢斗胆预言。相比之下,人类在漫长进化时间里,几千年基本可以忽略,所以以我学过的进化论知识,我可以预言3000年,其时人类最顶尖的围棋选手水平怕不能让现在的李昌镐二子,如果还是一样竞技规则的话。我感觉这个预言比老聂的劳什子预言要可靠的多。

纵观老聂以前的一些言论,多是没经过脑袋的,我更愿意相信他的出发点是因为对围棋理解日深,从而对围棋产生了一种敬畏心理,夸大其辞,做感叹之言论,未料被无良记者断章取义,哗众取宠。不过老聂能让人开心一乐,还是要谢谢了。

Excel 2007 Chart导出为image

南卓铜 ([email protected])

总结如下,

1) copy chart in Excel 2007,paste到paint(画图)或其它图像处理软件,然后保存为图像,如jpeg。

2) 写vba代码,ActiveChart.Export FileName:=’image.gif’, FilterName:=’GIF’,可以保存为JPG, GIF和PNG。需要知道如果在Excel里写VBA代码,并执行,有可能要求打开宏安全控制选项。但无法控制图像质量。MATLAB Handle Graphics

3) copy chart,打开PowerPoint,将chart使用special paste…到Powerpoint,可供选项有emf, png, gif, jpeg等。在powerpoint中右键点击粘贴进来的图像,可以save as picture…,保存为需要的格式。

4) 将包括chart的worksheet保存为了html,将同时生成.gif的chart图像。

5)第三方add-in,比如PUP 7,不过需要掏钱。

你还有什么tip?

Land cover数据Goode格式在ArcGIS中的访问

南卓铜 ([email protected])

一些有名的Land cover数据集并没有使用通用的GIS格式,比如在以下链接,

http://www.geog.umd.edu/landcover/tree-cover-poster/download.html

下载broadleaf Layer,gzip解压缩后为 694,417,920 bytes,扩展名为.img,但不同于Erdas Imagine .img格式。

据介绍Goode binary file是plain binary,即相当于BSQ (Band sequential),其metadata有以下信息:

40031 pixels by 17347 lines,且是8bit unsigned integer,应当总大小为 694,417,757bytes,而解压缩后的大小大了163bytes,目前我还不知道原因。

根据metadata信息,可以添加头文件,broadleaf.hdr,内容如下,

nrows      17347
ncols     40031
nbands     1
nbits     8
layout     bsq
skipbytes     0
ulxmap     -20015500.000
ulymap     8673500.000
xdim     1000
ydim     1000

将broadleaf.img改名为 broadleaf.bsq,在ArcGIS中才可以正常访问,也可以导出为其它格式,比如Grid。不过注意其空间参考系仍为unknown。在ArcCatalog中打开其属性页,在spatial reference中选projected coordinates/World/Goode Homolosine (Land).prj。可以将bsq设置为 Goode投影。

如果希望转换成其它常用格式,

image

上图将bsq转成ESRI binary grid格式。图形化结果如下图。

image

在Workstation下好像没有找到Goode投影,但在Desktop下可以找到.prj。因为如果bsq是unknown时,导成esri grid,在ArcCatalog里对grid进行定义投影(属性->define spatial reference),没法子定义到Goode投影,很奇怪。所以这里,我是先对bsq定义投影,再转成grid。此外,通过ArcToolbox的Data Management下的projection define来做,也是可以。

笑 (转)

林徽因

笑的是她的眼睛、口唇,
和唇边浑园的酒窝。
艳丽如同露珠,
朵朵的笑向贝齿的闪光里躲。
那是笑——神的笑,美的笑;
水的映影,风的轻歌。

笑的是她惺松的头发
散乱地挨着她的耳朵。
轻歌如同花影,
痒痒的甜蜜
涌进了你的心窝。
那是笑——诗的笑,画的笑;
云的留痕,浪的柔波。

朵朵和奶奶回家了

因为媳妇到呼和浩特开差,我妈不大会普通话,又带着朵朵,第一次坐飞机,我是比较担心。媳妇离开前,她跟我妈讲了一遍坐飞机的流程。小吴送她们去机场,又跟她们讲了一遍。我也很详细的讲了一遍,听到最后我妈说她都可以背下来了。

美国东部时间今天早上我给妈打电话,说已经顺得到家。还夸了一下朵朵,说找登机口的时候,是朵朵帮找到的,飞机上系安全带,我妈一下子还不会,是朵朵帮她系上的,还得意洋洋问奶奶说,你知道怎么解开吗?结果这么一掰就解开,我妈说她还真不知道。

我让朵朵接电话,本想夸她几句,结果接上问,是妈妈吗?我说是爸爸。她居然不接了,将电话又还她奶奶了。不过听说回家已经跟我弟的儿子玩到一起了,玩的很高兴。这就好。

原来我们是吃兴奋剂长大的

看到腾讯的新闻,

“奥运期间被禁止销售常用药名单”——奥运将至,根据国家食品药品监督管理局的规定,在全国范围内规范含兴奋剂生物和中药制剂的生产、经营和销售,在零售药店除胰岛素以外的蛋白同化制剂和肽类激素,一律停止销售。这意味着,皮炎平、派瑞松、咳喘宁、大活络丸等百姓常用药被禁止销售。……

“药监局规定奥运期间停售部分药品”——包括了平常用于止咳的急支糖浆、复方川贝止咳糖浆,都被列入被禁售的药品名录里了,就连几毛钱一板用于感冒的感冒胶囊,这次也被停止销售,而恰恰这些药里都含有能够兴奋中枢神经的麻黄碱。

怪不得国人都自信满满的,不少年轻同志还立志解救美国人民于万恶资本主义的水深火热之中,原来是因为我们都是吃兴奋剂长大的。我比较担心的一点是,奥运过后,这些兴奋剂药就回到药店正常销售了。

北京朋友都在感叹这段时间北京变化真大,变得漂亮、安全或者说“和谐”,这是真心话,不是在讽刺,尽管奥运的严格管制带来一些不便,但明显的总体上讲对北京看到的好处更多。我感叹于原来一些事情不是说就没有法子做好,而是平时根本没认真去做。中国很多事情都一样道理。

大量时间序列快照文件到时间序列文件相互转化的一个思路

南卓铜([email protected])

想像这样一种情境,我们有一个空间网格文件,分辨率为30×30 m,空间范围为30x30km^2,每个格子一个数据的话,也就是说有1,000,000个数据(假设是double类型)。每个网格文件可被认为是研究区域时间序列的一个快照。整个时间序列为11年(1997-2007)逐小时,即多达96,408个网格文件。我们要做的是,将全部的时间序列快照文件转换成每网格上的时间序列文件(将形成1M个文件)。类似的,我们还要做的是全部网格的时间序列文件时间序列快照文件的转换。

全部保持文件打开是不可能的,因为文件句柄是很宝贵的资源,每个进程都只允许有限个打开的文件(比如512个)。通常的迭代思路是,打开一个时间序列快照文件(网格文件),取其中一个网格上的数据,打开此网格的时间序列文件,写入时间和对应的数据,关闭两文件。这意味着需要打开和关闭文件1Mx96,408次。大家知道磁盘I/O的开支是十分昂贵的。此迭代方法将需要极长的时间来完成(一周甚至几周,取决于硬盘读取数据)。而内存却没有得到有效的利用。

对于内存十分大(比如64G内存)的服务器,也许可以将全部的待转换文件一次性装载到内存,每个在内存内分析,组合成输出格式,再一次性格式化输出。对于一般的个人机器此方法不通。

为了优化速度,我们采用这样的解决方案。针对时间序列文件时间序列快照文件(网格文件)的转换,

在内存生成n个时间序列快照文件,每个网格上用Missing data来填充 (内存主要消耗在这里)
打开一个时间序列文件
    从时间序列文件读n行(每行包括时间和对应的double值)
    将读取的n个数据,写入对应的n个快照文件的对应位置
关闭此时间序列文件,迭代
将内存内的n个快照文件写出
从n+1的位置开始再迭代以上过程,一直到结果。(关键)

最后一次可能不正好等于n,需要程序作相应控制。n的取值需根据运行的内存情况进行高速。对2G内存的工作站,n取为5000-10000。内存越大,n值可以取越大,可以有更好的执行效果。

当然,也可以同时打开多个时间序列文件,以最大化优化性能,但带来的是迭代控制上的复杂。而且据我的有限测试,同时打开多个时间序列文件,性能并没有得到明显改善(可以理解,因为磁盘I/O的存取本质上讲是磁头的顺序读取,由同一个磁头臂来控制)。

其中需要注意的是地方,是如何控制下一次准确快速定位到n+1的位置上。时间序列文件是文本文件,顺序读取在性能上很受影响。比如在最后一个循环时,将先遍历前面的全部行,然后才到达需要的起始行。我们需要以二进制形式打开,并自行控制每个时间序列文件的起始读取位置(各文件位置可能不一样,由于每行数据长度不等)。在c#里,以StreamReader打开,无法通过base stream取得当前准确的位置(position)。我们构造了TimeSeriesDataFile类。初始化需要给定文件名和起始读取的位置。ReadLines函数可以返回给定数目的数据行,通过CurrentPosition属性取得下一次读取的起始位置。

using System;
using System.Collections.Generic;
using System.Text;

namespace nzt.TimeSeries2Spatial
{
    /// <summary>
    /// Access time series data text file as binary.
    /// </summary>
    class TimeSeriesDataFile
    {

        private string _filepath;
        private System.IO.FileStream _fs;
        private long _lastpos;
        private long _startpos;
        private const int MAXLINELENGTH=50; //bytes, ensure it larger than max length of each line.

        public TimeSeriesDataFile(string filename, long startposition)
        {
            _filepath = filename;
            _startpos = startposition;
            Open();

        }

        private void Open()
        {
            _fs = new System.IO.FileStream(_filepath, System.IO.FileMode.Open, System.IO.FileAccess.Read);
            _lastpos=_fs.Seek(_startpos, System.IO.SeekOrigin.Begin);
        }

        /// <summary>
        /// Read a number of lines from stream beginning at startposition.
        /// </summary>
        /// <param name="count">Number of lines to be expected to return</param>
        /// <returns></returns>
        public string[] ReadLines(int count)
        {
            if (_fs == null) return null;

            byte[] buffer = new byte[MAXLINELENGTH * count];

            int bytesRead=_fs.Read(buffer, 0, buffer.Length);

            if (bytesRead == 0) return null;

            //we have data in buffer now.
            List<String> sb_list = new List<String>();
            int c = count;
            StringBuilder sb = new StringBuilder();
            int i;
            for (i = 0; i < bytesRead; i++)
            {
                if (buffer[i] != ‘r’ && buffer[i] != ‘n’)
                    sb.Append((char)buffer[i]);
                else if (buffer[i] == ‘n’)
                {
                    sb_list.Add(sb.ToString());
                    if (–c<=0) break;
                    sb = new StringBuilder();
                }
            }
            if (c>0 && sb.Length>0) sb_list.Add(sb.ToString());

            _lastpos += i+1;

            return sb_list.ToArray();
        }

        public long CurrentPosition
        {
            get { return _lastpos; }
        }

        public void Close()
        {
            _fs.Close();
        }

        ~TimeSeriesDataFile()
        {
            Close();
        }

    }
}

对于时间序列快照文件(网格文件)时间序列文件的转换,应用同样的思路。但由于时间序列快照文件(网格文件)一般较小,比如几百KB(相对,11年的逐小时时间序列文件则到2MB以上),则无须对StreamReader进行改造,可以一次性load到内存,在内存进行定位分析。

如果大家还有好的解决方案,也请分享。