GL终究还是来了

本文算是最近三周内的第三个“来了”系列(博主是不是整天无所事事?)…不过这次的确是因为工作需要的原因,GL的主题将会持续一段时间。至于为什么我习惯称之gl,那么请留意gl的发音/gl/,前面再加个open…总会有人真正体会到她的魅力所在的。不过这也确实不是我第一次尝试学习gl了:一年前我曾冒着烈日学习nehe怎么用gl画会转的三角形…两天后被迫放假回家,于是又练了个新号,直到前几天又被个文明贼盗了2w多g…这里似乎离题了!总之现在就算不上是gl新手了吧…

上述事实再次告诉我们,当你像我这么老的时候,好奇心就不能再是生活的第一目标了,而应优先满足需求…在满足各种实际需求的过程中,我们永远的好奇心也将逐步得到慰藉。

PS:Kinect for Windows SDK近期连续跳票…广大Kinect开发者表示情绪灰常稳定,是到了该OpenNI收拾残局的时候了。

PPS:本周读书去…

Comments

CUDA学习的现状和困扰

首先需要说明的是,CUDA并非GPGPU(GPU通用计算)的全部。但是我们仍然应该首先感谢Nvidia,因为CUDA使GPGPU走向大众化,使得我等终于也能从中一窥究竟了。不过随之而来的疑问却是,CUDA真能担负得起GPGPU的重担吗?至少在目前看来,答案并不太明显。

战前准备

当然,要使用到CUDA,最好还是拥有一块N系列的显卡,尽管CUDA支持emu模式(该模式将把cuda程序完全放在CPU中运行),但device加速计算才是我们真正的目的。拥有并装配任意一块GeForce 8、9、GTX或Quadro、Tesla系列显卡后,在NVIDIA的网站上下载最新适用的driver、toolkit和sdk,安装在相应的平台上就可以了。下图简示了cuda应用程序的软硬件体系结构:

[singlepic id=31 w=320 h=240 mode=watermark float=center]

图中我们可以看出GPU在现行体系结构的计算机中所扮演的角色,以及cuda应用程序驱动GPU设备的控制流。上述环境搭建好后,就可以进入程序设计阶段了(的确,这才刚刚开始)——实际上当前cuda应用中还普遍存在一个问题,即在当前ide中配置应用程序框架以便于开发,针对该问题我们将在今后的文章中做一些说明。

CUDA的extended C基本组成

1、基本结构和kernel function

这里给出一段简单易懂的基于cuda的C语言程序模板:

[c]

Main(){

//Allocate memory on GPU

float *Md;

cudaMalloc((void**)&Md, size);

//Copy data from CPU to GPU

cudaMemcpy(Md, M, size, cudaMemcpyHostToDevice);

//Call GPU kernel function

kernel<<<dimGrid, dimBlock>>> (arguments);

//Copy data from GPU back to CPU

CopyFromDeviceMatrix(M, Md);

//Free device matrices

FreeDeviceMatrix(Md);

}

[/c]

从这段程序中我们可以看出cuda编程的两大特点:一是在显存中的空间分配操作和CPU操作内存的思路基本一致,同时避免不了Host与Device之间的大规模数据交换;二是GPU计算的算法核心在于核函数kernel function的设计。通过核函数的调用方式kernel<<<dimGrid, dimBlock>>> (arguments)可以看出,kernel是自定义的核函数名,符号“<<<”和“>>>”中含有两个参数,分别是Grid维数和Block维数,在cuda中Grid表示一个block组,而Block表示一个thread组,显然参数dimGrid和dimBlock分别定义了Grid和Block的规模,其维数均可达到三维。下面是一个定义Grid和Block的例子:

[c]

dim3 dimGrid(2, 2);

dim3 dimBlock(4, 2, 2);

[/c]

应当注意,每一个thread block最多只能分配512个threads,其中每调用一次kernel函数就相当于在GPU上分配了一个Grid。此外,为了使线程内部处理各自的数据块,核函数内部还有为每个线程分配的threadIdx、blockIdx、threadDim和blockDim信息,示例程序如下:

[c]

//CUDA program

global void inc_gpu(float a, float b, int N)

{

for (intidx = 0; idx<N; idx++)

a[idx] = a[idx] + b;

}

int idx =blockIdx.x* blockDim.x+ threadIdx.x;

if (idx < N)

a[idx] = a[idx] + b;

}

void main()

{

dim3 dimBlock (blocksize);

dim3 dimGrid( ceil( N / (float)blocksize) );

inc_gpu<<<dimGrid, dimBlock>>>(a, b, N);

}

[/c]

2、类型关键字

为了处理Grid、Block、Thread层级内的内存共享和线程通信问题,cuda特别给出了自有的声明关键字:

device

该空间储存于GPU上,属于全局可读写空间,和应用程序具有相同的生命期,且可被grid中所有线程存取,CPU代码通过runtime函数存取。

constant

该空间储存于GPU上,属于全局只读空间,和应用程序具有相同的生命期,且可被grid中所有线程存取,CPU代码通过runtime函数存取。

shared

该空间存储于GPU上Block内的共享存储器,和Block具有相同的生命期,只能被Block内的全部线程存取。

Local变量

储存于SM内的寄存器和local memory中,和thread具有相同的生命期,Thread私有。

同时CUDA还引入了几个新的内置数据类型,除了前面已经给出的dim3外,还针对C语言的内置数据类型进行了进一步扩展,如int4四维向量,其四个维数分别存储于int4.w、int4.x、int4.y、int4.z中。此外还有一个独特的Texture type,其一般声明表达式为texture<Type, Dim, ReadMode> texRef,其中type为C的内置数据类型,维数最多为3,ReadMode包括了cudaReadModeNormalizedFloat和cudaReadModeElementType两种形式。

CUDA在函数声明中也引入了新的关键字:

device float DeviceFunc()

该类型函数只能在GPU内部被调用,且只能在GPU内部运行。device函数不能使用&取存储器地址、不支持递归、不支持静态变量,也不支持变长参数的声明方式。

global void KernelFunc()

一般为核函数,在CPU内调用,但在GPU内运行,且返回值只能为void。前文已经详细说明了核函数的一般使用方法,应注意核函数还拥有一个配置参数,其调用方式一般如下:

[c]

……

size_t SharedMemBytes = 64; // 64 bytes of shared memory

KernelFunc<<< DimGrid, DimBlock, SharedMemBytes >>>(…);

[/c]

host float HostFunc()

只在CPU中调用和运行的程序。

还有一种声明方式,即devicehost 组合使用,那么函数分别在CPU和GPU中被编译保存。

CUDA还有自己一套非常实用的数学函数库,其功能包括pow, sqrt, cbrt, hypot, exp, exp2, expm1, log, log2, log10,log1p, sin, cos, tan,asin, acos, atan, atan2, sinh, cosh,tanh, asinh, acosh, atanh, ceil, floor, trunc, round, etc。但是上述函数都不能像matlab中一样进行向量运算(但是有了GPU和Kernel…呵呵),此外在函数名前方加入__前缀就可以达到更快的速度,但会损失一定精度。

CUDA Libraries

如果只有CUDA drivers API,那么我们对CUDA的研究可能真的要止步于此了,因为基于GPGPU的程序设计目前还缺乏方法论支持——当然这是科学家的义务(如果我们要想在此做一些贡献,可能不会先发在blog上:)。不过正是因为更多有趣的Libraries,使我们不得不对CUDA继续感起兴趣,尽管这种编程方式的难度与matlab相比远高了许多数量级。笔者使用的CUDA Toolkit 4.0 RC2(发布于2011年4月)大体上包含了四个CUDA应用库,包括GPU-accelerated BLAS library基础线性代数子程函数库、GPU-accelerated FFT library快速傅里叶变换函数库、GPU-accelerated Sparse Matrix library稀疏矩阵函数库、GPU-accelerated RNG library随机数生成函数库。

最常用的函数库是CUBLAS基础线性代数子程函数库, 其基本功能是进行矩阵运算。有经验的读者可能会想到直接用CUDA API编写矩阵运算程序可能并不复杂,但是鉴于面对的是建立在本身就不怎么友好的标准C之上的CUDA,CUBLAS还是吸引了众多关注。这里给出基于CUDA 4.0 CUBLAS(CUBLAS v2.0)程序的基本结构:

[c]

cublasHandle_t handle ;

stat = cublasCreate(&handle ) ;

stat = cublasSetMatrix ( M , N , s i z e o f (* a ) , a , M , devPtrA , M ) ;

modify ( handle , devPtrA , M , N , 2 , 3 , 1 6.0 f , 12.0 f ) ;

stat = cublasGetMatrix ( M , N , s i z e o f (* a ) , devPtrA , M , a , M ) ;

cudaFree ( devPtrA ) ;

cublasDestroy ( handle ) ;

[/c]

这段代码摘自CUBLAS Library手册中的Example code,其功能是将当前矩阵改为以1为起始索引。handle是CUBLAS新版本中加入的线程安全句柄,这保证了CUBLAS库函数的可重入性,从而允许其在CPU的不同线程中被安全调用。需要注意的是,CUBLAS采用了以列为主序、1为起始索引的存储结构,以便更大程度上兼容Fortran编程环境,然而数据组织这种方式与C/C++完全不同,这就需要在向存储器中填入数据时考虑选择列为主序的问题,其中官方手册上也做了如下说明:

当应用程序的形式是在C中嵌入的Fortran代码,那么为了适应索引由1起始的习惯和以列为主序的存储方式,可使用以下宏调用进行数组操作:

[c]

define IDX2F( i , j , ld ) ( ( ( ( j )-1) *( ld ) )+(( i )-1) )

[/c]

如果使用的是本地C/C++代码,那么可使用以下宏调用进行数组操作:

[c]

define IDX2C( i , j , ld ) ( (( j) * (ld))+( i) )

[/c]

上述两段宏定义的“id”是指矩阵的主维数,如当矩阵是以行为主序存储时,主维数id即为列数,反之亦然,主维数在CUBLAS库函数中被大量使用。

后记

CUDA入门时问题往往集中在基本的extended C以及库函数的调用,其次就是编译、链接和调试了。但是当我们真正熟悉它的基本使用方法,才发现CUDA的难度集中在算法和程序设计乃至优化等方面,要知道,在传统CPU上许多数值计算方法都尚未得到高效的解决,更何况要让众多的开发者把精力投入到这种新的思考模式——也许这正是AMD至今把GPGPU捧在科学研究之巅的原因。目前,CUDA在GPGPU上的主要应用还是基本的数值计算,我们知道在纯粹的浮点运算上GPU已经把CPU远远抛在了后面,但处于冯诺依曼体系下的GPGPU仍然需要唯CPU马首是瞻,这就引出Host和Device间数据传输瓶颈的问题。即使是在PCI-E 2.0理论传输峰值达到8.0g/s的情况下,仍然与高端GPU动辄上百g/s的片内存储带宽相距甚远。更何况CUDA的数据传输机制为了避免Host对Memory的影响,往往先在Memory中申请一块独占空间再进行DMA传输——这也是往往使用cudaMallocHost的page locked内存空间申请模式的原因。

因此,我认为,如果作为一名普通爱好者,那么学习CUDA、偶尔加速并优化自己的实际工作——如快速视频编解转码,也小有一番风味。但如果热衷于高性能计算特别是GPGPU的研究,那么是否采用CUDA——这个让GPU走下神坛的“利器”并借助它去接近理想,恐怕还是有待商榷的。

Comments

啃奶来了—Kinect for Windows抢先体验

如果说微软将“被迫”开放Kinect for Windows SDK,无疑这次是PrimeSense立了首功。然而要不是我冒着生命危险在鸟不拉屎的地方蹲守了近十天,终于在一个学院林立的院落里拿下Kinect for Taiwan Only,并跋涉1200公里侥幸逃生的话,本文的出炉可能还要推迟一个月左右。总之无论如何,感谢MS,感谢PS,感谢万恶的铁道部,我终于活着且独立完成了这篇文章。

Kinect到底好不好玩,我不知道。不过我知道为了测试一个hand tracker手臂举得酸痛,然而结果就像是我们数年前的理想如今突然变成了现实,而且这种久逢甘霖的冲动仍在持续,套用某牛逼CEO创业时的一句话:“说实话,我们并不清楚它未来将会变成什么。”由于微软将于2011年5月16日正式开放Kinect for Windows SDK beta,尽管是non-commercial use,但对我们来说已经足矣。本文在某种程度上是基于已有hacker的经验谈,基本框架一般是OpenNI+SensorKinect+NITE这样的半官方组合,而要想全面入门体感开发还是需要先研读一下OpenNI UserGuide这种小部头。

下面是啃奶来自NiViewer的第一张深度对比图:

[singlepic id=30 w=320 h=240 mode=watermark float=center]

根据PS的说法,自然交互(Natural Interaction)就是一种基于声音和视觉等人类体感信息的人机交互方式。目前看来,NI设备完全可以取代如遥控器、鼠标等远程控制外设。当前日常生活中常见的自然交互方式包括:

演说以及口令识别,NI设备通过声音接收装置接收语音指令;

手势识别,通过预设计一定的交互范型,人们可以通过简单的手势控制一般设备,如卧室的光照强度;

躯体运动跟踪,NI设备通过识别当前图像中的完整人类躯体,快速分析躯体骨架信息,并将实际躯体运动姿态转换为骨架运动信息。

值得注意的是,前面的图中展示了一张Kinect的深度对比图,也就是说大多数NI设备通过技术手段可以感应物体的空间深度信息,目前已有利用Kinect进行三维场景重建和空间交互的hacker例子。下图是利用NITE的StickFigure追踪人体骨架姿态信息的例子:

[singlepic id=28 w=320 h=240 mode=watermark float=center]

简单修改一下NITE的PointViewer例子,我们可以把手势跟踪进一步改进为手势鼠标控制(已做消除抖动处理),如下图:

[singlepic id=29 w=320 h=240 mode=watermark float=center]

由于目前官方SDK尚未推出,基于Kinect的开发框架方案已有许多,但缺乏权威性。较令人信服的组合除了我们上述建议的三件套,再加上OpenCV算法包和OpenGL库就足以做出十分不错的应用了。即日起本站新增加Kinect项目板块,将重点关注基于Kinect的应用开发。值得一提的是,目前在kinect开发领域拥有一个人气不错的中文社区http://www.cnkinect.com%EF%BC%8C%E8%BF%99%E5%9C%A8%E5%9B%BD%E5%86%85%E5%AD%A6%E9%99%A2%E6%B4%BE%E6%97%A5%E7%9B%8A%E8%A1%B0%E5%BE%AE%E7%9A%84%E4%BB%8A%E5%A4%A9%E7%BB%9D%E5%AF%B9%E6%98%AF%E9%9A%BE%E8%83%BD%E5%8F%AF%E8%B4%B5%E7%9A%84%EF%BC%8C%E4%BB%8A%E5%90%8E%E6%88%91%E4%BB%AC%E4%B9%9F%E5%B0%86%E6%8C%81%E7%BB%AD%E5%85%B3%E6%B3%A8%E8%BF%99%E4%B8%80%E7%A4%BE%E5%8C%BA%E7%9A%84%E5%8A%A8%E6%80%81%E3%80%82

Comments

真累啊

本来就懒得跑,没消停几天就又得出差…人生真真是个杯具。

Comments

Web开发的若干安全问题小结

 万恶的项目终于要结题了…这时无论如何确该总结一下。实际上这次项目最大的收获有两部分,本文仅介绍其中之一—Web开发中的一些安全问题及其解决策略。

 项目的安全问题起源于首次测试版本提交。测试方使用了IBM Rational Appscan扫描提交的两部分应用程序,每个应用程序返回了将近二十个安全问题,具体列表这里先省略。Appscan是IBM旗下知名的web安全自动化测试工具,同类产品中较知名的还有HP的Webinspect。测试方回复的意思是说所有安全问题在项目验收前必须全部解决,而事实上本项目从设计阶段起就没有考虑过应用程序级的安全问题,对方也从未提出过此类要求。不过我们还是尝试先在逻辑上处理部分问题,自测后发现所有问题依然存在,包括众所周知且已有很多解决方案的Injection。

 纠结几日后大家终于决定重构,考虑在重构过程中再应用安全策略,于是又花了几周时间问题才得以基本解决。然而在解决web安全问题的过程中我们查阅了较多资料,由于缺乏经验,网上的许多解决方案基本都是无效的。当然也不能排除测试手段随着新的漏洞出现而加强的可能。这里我们仅就使用Appscan 7.8.0.2 Rules-Update 947的测试报告针对一些常见问题给出相应的解决办法。

 应当注意,目前Web安全领域拥有自己的国际性非盈利组织OWASP(Open Web Application Security Project),大部分Web安全问题的讨论都可以在https://www.owasp.org/index.php/Main_Page上找到相应的文章。如果读者对Web安全感兴趣,那么可以直接转向官网,本文接下来只考虑如何让应用程序通过Appscan 7.8的安全测试:)(如果仅仅是服务器或tomcat配置的问题,那么本文将不涉及)

 Injection

 SQL盲注问题十分普遍,一般考虑前端+后台验证。我们最初在后台对每个表单进行了关键词检测,然而测试依然不通过,而且仔细查阅资料后发现我们的过滤规则要比一般规则复杂得多,唯一问题在于处理策略不同:我们是直接给出错误提示,而大部分策略是将关键词截断后交付处理。郁闷的是改进后依然不能通过测试。于是开始求助,发现许多专业的JEE开发人员基本都直接使用一种servlet过滤器,项目地址http://securityfilter.sourceforge.net/,更神奇的是,配置该过滤器后问题的确不存在了,但是从表单直接进行关键词提交后台依然能收到完整语句,比较奇怪。此外,另一个问题Cross Site Scripting跨站点脚本攻击由于此过滤器的出现也不复存在了。

 Inadequate Account lockout

 不充分的账户封锁,这个问题看起来很简单,就是为了防止黑客暴力破解用户密码。尽管目前大多数站点使用的验证码技术已经足以克服此类问题,但Appscan测试对于验证码确实是无能为力的,因此测试方要求取消应用程序的验证码功能。后来使用了错误次数限制方案,如果某用户登录失败超过若干次,则直接禁止用户再尝试登录,直到申请管理员解封为止。

 Session Identifier Not Updated

 会话标识未更新,Appscan给出的描述是建议用户每次登录时需使用新的会话标识。网上给出了一种很简单的解决方案:

[c] session.invalidate();

 Cookie cookie = request.getCookies()[0];

 cookie.setMaxAge(0);[/c]

 上述代码是解决问题的核心,但不会使用照样无法解决任何问题。我们最终的解决方法是限制每次停留在登录页上的时间,如果超时提交登录则返回页面更新session后要求重新登录。另外,对于任何登录失败、用户注销、意外退出等情况都要彻底清除一遍sessionid。最后为了彻底解决这一问题,每次用户进入login页面后先检查sessionid是否更新,否则跳至过渡页清空sessionid后再返回login,问题解决。

 Unencrypted Login Request

 已解密的登录请求,要求就是数据加密传输。开始想当然地编写了一种简单的前端加密算法,再从后台进行解密,然而Appscan依然报此错(现在感觉是测试软件本身的问题,最后再谈)。于是决定转https,JEE平台转https的通行步骤如下:使用java keytool生成RSA密钥,打开tomcat的8443端口,并配置相应密钥即可,具体步骤网上很多,GOOGLE即可。更麻烦的问题随之而来,测试方认为尽管提供了8443的https访问,但仍然还可以通过8080使用普通的http访问,并要求禁止普通模式登录。

 但是由于要把8080端口提供给项目的WebService组件,不太可能完全禁掉。这里在项目web.xml里添加了一些配置解决:

[c]<security-constraint>

<!– Authorization setting for SSL –>

<web-resource-collection >

<web-resource-name >SSL</web-resource-name>

<url-pattern>*.jsp</url-pattern>

<url-pattern>*.action</url-pattern>

</web-resource-collection>

<user-data-constraint>

<transport-guarantee>CONFIDENTIAL</transport-guarantee>

</user-data-constraint>

</security-constraint>

<login-config>

<!– Authorization setting for SSL –>

<auth-method>CLIENT-CERT</auth-method>

<realm-name>Client Cert Users-only Area</realm-name>

</login-config>[/c]

 应注意,由于项目的一些组件无法通过https,因此url-pattern字段只对.jsp和.action进行了限制,如果不做特定限制,则系统默认是全部使用https传输。而且上述设置一旦在某个工程中出现,那么当前tomcat将全局采用这一配置。(这是实验结果,目前不评价其合理性)

 Cross-site request forgery

 跨站点请求伪造,网上研究文章一大堆,同时也是OWASP 2010 TOP 10级别。问题是很久都无法根本解决,经过两周研究后终于得以通过Appscan的CSRF检测。其实基本思路很简单,即在用户的提交的表单中加入一个隐藏域,存储一个唯一验证标识码:如sessionid。用户提交表单后由后台首先验证sessionid与提交标识是否一致,否则报错。

 然而要过Appscan的检测还是不太容易。首先是针对所有action的sessionid验证标识码检查,例如login,如果在session中始终保留一个标识,那么可能会报Session Identifier Not Updated错,解决方法是在登录成功后修改一下验证标识码。其实许多无法彻底解决的原因在于报错处理,起初我们将此类错误全部统一转向了一个安全报错页面,后来始终都无法解决这一问题。后来发现如果将错误情况转至首页则Appscan就不会报错,这可能是一个软件特征匹配的bug。

 除此之外,包括许多tomcat的配置问题这里就不再一一说明了。本次项目所引出的web安全问题确实需要得到重视,然而我们也应看到,过渡迷信流行的自动化测试工具可能会导致在某些问题上缘木求鱼,尤其是对于Appscan这种判断防御行为的复杂软件,仅靠有限的规则设置就当做是web安全的唯一标准——这也并非是聪明人的选择。

web
Comments

即日起开始同步分发blog数据

 本博即将迎来第300篇文章了。重建三年多来我做过多次较大的改动,这一版本无疑是维持较久的(已超过一年时间)。目前觉得对大家有所帮助的还是有关于技术方面的文章,但显然技术需要“交流”才能进步。不过有趣的是本博自诩为“交流频道”,却几乎没有相关交流的内容——这其实与我始终坚守这个独立的小站点有关,更何况百度根本不买我等穷弱儒的账。

 然而独立建站是我十余年的信念,岂能轻易废止。在300篇纪念之际我决定争取彻底改变多年来封闭写作的习惯,将原创文章第一时间推送至renren.com和qqzone,同时进一步向google提交新摘要数据并加强seo。另一方面,由于时间的原因,我能拥有的亲自撰写文章的时间越来越少,但潜在的阅读量仍十分可观。因此决定今后在三个不同的blog上提交侧重不同的数据:

 http://www.hanyi.name mp77的技术交流频道:主站点,注重原创文章,实际上目前已很少有转载内容了。

 http://8621316.qzone.qq.com QQ空间日志:重点转载文章,内容将多样化。

 http://blog.renren.com/micropotato 人人网日志:同步hanyi.name数据,并记录转载+评论的短篇文章。

 采取以上策略的原因在于,由于长时间逐渐减少,在短时间的条件下通过压缩原创篇幅来扩大信息量。如果有朝一日非主站点的数据量足够吸引人的话,我仍会选择将其中数据反导入到主站点中。

可爱的TSP

 有意思的是,我接触到的第一个算法问题就是TSP(旅行商问题),自此以后我开始热衷于计算机科学,遗憾的是至今也未有什么收获。但是很显然,TSP及其衍生问题始终是横亘面前的一座高峰,对每一个人—包括所有计算机科学家,都是那么的不可逾越。也许真要等超越阿兰图灵的大师出现—我等凡人恐怕是看不到这一天了。

 不过这里仍要说明的是,我们学习和体会TSP,并不是为了吃饱撑的自虐,也非为了搬弄别人整出的玩意来给自己凑文章。TSP的解法之多令人惊讶,各种方法—即使是一名高中生和一位大学教授所掌握的,在TSP面前似乎都一律平等。但它为我们巩固和发掘基本算法提供了确切的衡量准则,我想这才应该是凡人研究TSP的真正目的。

 TSP问题最初的形式来源于爱尔兰数学家W. R. Hamilton发明的智力迷绳游戏,即在一张无向图中找到一条Hamilton回路。也就是从无向图中任意起点处找到一条路径,该路径唯一地穿过图中的每个顶点,最终回到起点。通常Hamilton回路也被称为Hamilton圈。经过若干年演进,普林斯顿大学的Hassler Whitney首次将该系列问题命名为Travelling Salesman Problem,即TSP,此时问题不仅限于找到一条回路了—而是要求找到花费尽可能少的解。1972年,伯克利的Richard M. Karp证明了Hamilton回路问题是NP完备的,那么作为Hamilton回路的扩展—TSP问题也至少具有NP难的特性。截止目前,应用传统方法解决TSP问题的有动态规划法、贪心法、回溯法、分支限界法、线性规划方法以及采用混合策略;而对于在有限时间内寻找尽可能优解的大规模TSP问题,目前更多的倾向于采用随机组合优化搜索方法,如模拟退火、禁忌搜索、神经网络、蚁群算法、遗传算法以及混合策略等。

 应当注意的是,一个实用的TSP解决方案,通常包括了两个阶段的改进算法,即第一步的初始路径生成和第二步的优化算法。当然,由于本文并非是研究某种特定的算法及其改进,因此并不会突出介绍其中某一个阶段。下面我们就常用的解决方法进行介绍并做简单评价。

 动态规划

 DP思路清晰,易于实现。显然很方便就能找到TSP的最优子结构,从而给出基本的状态转移方程Mk(i,V) = min(Mk-1(j,V{j}),dij)。需要指出的是,DP解决这一类图论问题都具有较高的时间和空间复杂度,一般可用于求解问题规模较小的TSP问题。

 贪心思想

 Greedy为我们贡献了许多著名的图论算法,如Dijkstra的SSP、Prim和Kruskal的MST等等。然而贪心多用于求近似解,而非定格于最优值。目前贪心思想常用于辅助随机搜索算法,类似一种启发式机制,但初始解在许多情况下与最优解相距较远。

 回溯法

 一般的TSP问题可以简单归纳为排列树问题,那么应用回溯法的确可以得出最优解。问题在于,回溯法的解空间非常庞大,一般的剪枝方法并不十分有效,其时间复杂度上界甚至达到n!。然而通过DP思想的启发,我们可以尝试设计一种自顶向下的动态规划算法解决TSP问题,即备忘录方法。备忘录方法相当于搜索和DP的一种折衷,但也丝毫没有改善较高空间复杂度的问题。

 分支限界法

 BB也是一种解空间搜索方法,但和回溯法有明显区别。回溯法通常是一种DFS搜索方法,而BB实际上是一种BFS方法。其中采用队列式的BB方法本质上就是BFS,只不过添加了一定程度上的剪枝策略。而采用优先队列的BB倾向于寻找一个最优解,在许多最优化问题中,如果限界函数和优先队列设计得当,那么BB的性能要优秀许多。然而不可否认,BB的时间复杂度上界实际上并不比回溯法低。

 整数线性规划

 线性规划方法本身就要比TSP问题复杂得多,这里不得不提LP的原因在于,诞生以来计算机科学本应是崇尚应用的科学,尽管学科本身目前还比不上数学、物理学的历史悠久和博大精深,但后者如今基本上需要计算机来完成从理论到实践的跨越。我认为如果一位科学家只热衷于研究理论,那么他应被称作是数学家或者物理学家。

 这里所说的整数线性规划ILP解决TSP问题,实际上是由美国兰德公司的科学家在上世纪50年代提出的。如今ILP被用于许多学科,最直观的—如现代物流专业,TSP是其必须解决的基本问题之一,目前ILP已成为运筹学中的专门课题。ILP问题描述如下:

 试求min(cx),使其满足Ax=b,其中A为m*n矩阵,x为n维列向量,x>=0。

 根据数据特征的不同,ILP问题又可分为0-1ILP问题、混合ILP问题和纯ILP问题。一般的TSP问题可描述为典型的0-1ILP问题,故可用各种基于线性规划单纯形算法的割平面方法尝试解决。

 实际上,针对大规模TSP问题,现有的确定性方法已经很难在规定时间内得出较好的解—当然,如果不计耗费,那么上述算法恐怕还是得到最优解的截至目前的绝对保证,其花费甚至达到了数百CPU年。那么,许多智能算法的大行其道也就有其合理性了。唯一的问题在于,大部分智能算法来源于生物、物理学仿真和人工智能,学习和理解起来并不容易,更何况要用于解决实际问题。况且当前智能算法的“遍地开花”也表明,此类方法的稳定性和准确性还需得到进一步检验和完善。

 由于时间紧张(万恶的项目居然又延期…),我们这里无法再就智能算法做进一步展开了。下个月的文章里,我们将尝试介绍并实现一种应用确定性方法和智能算法解决类TSP(难度超过非对称TSP)问题的例子。

Comments

忙II

 果然是忙I的升级版,半夜偷偷跑回家逛逛,天亮就得回实验室了。下周的文章主题是古老与永恒—关于tsp问题。

Comments

样条插值的一般方法(下)

 前一篇文章中,我们通过基本的三次Bezier曲线引入了样条曲线插值的一般方法—catmull-rom及其改进算法Bessel—Overhauser。这里首先需要澄清的是,本文题名为“样条插值的一般方法”,实际上仅限于在计算机图形学中作基础应用性的讨论,而非是对数值计算方法中有关插值或曲线拟合的一般性总结,事实上在现阶段研究大量的理论方法对实际工作而言并无多大意义—当然,如果读者来自数学相关专业,本文的内容可能没有任何帮助,甚至起到误导的作用。

 同时,在继续这次的讨论之前,我们需要对上节遗留的问题做一些补充。因为我们愈加注意到,当前的基本造型技术所采用的样条插值方法实际上是改进后的catmull-rom,许多文献介绍了各种各样的改进算法,如果读者确实需要更加灵活可控的样条插值,除了上文文末的参考文献外,可能还需要阅读但不仅限于Kockanek和Bartels等人的论文。

 有了上述样条曲线插值的基础,就可以尝试进行曲面插值了。在这之前,我们应先回顾一下Bezier曲面的构造方法。

 三次Bezier曲面的控制点有4*4个,从基本公式来看,单独固定u值或v值,所得到的截面曲线实际上是普通的三次Bezier曲线。下面的影片演示了三次Bezier曲面的边界曲线构造过程,完整空间的情况比较复杂,待webGL成熟后我们会改进演示……

 

 类似基于分段Bezier曲线的样条曲线插值,同样可以使用局部Bezier曲面拼接的方法构建样条曲面插值。关键问题在于如何尽可能确保曲面拼接的连续性。由已知,两个Bezier曲面在边界处C1连续,实际上即其在边界处的偏导数非零且相等。现在给定m*n个控制点,要求构造通过每个控制点的C1连续的Bezier曲面。分析如下:

 1、按照catmull-rom的思路,易知可对相邻两个控制点插值构建Bezier曲线。

 2、又由前文的结论得知,三次Bezier曲面在边界处完全是一条只受边界点控制的三次Bezier曲线。

 由1、2,将mn个控制点划分为相邻的22角点集,每个角点集内部的4个角点即局部三次Bezier曲面片的顶点。通过Catmull-rom或Bessel-Overhauser算法在每条边界插值2个额外控制点,从而构建4条边界曲线。为了保证相邻曲面片C1连续,需要考虑局部面片中间的4个控制点。

 这里引入三次Bezier曲面在角点处的扭向量的概念:对于三次Bezier曲面,扭向量是指曲面在角点处的二阶混合偏导数。为了保证构建的曲面在边界处偏导数相等,可以将连续条件公式做简单变换,并使用扭向量表示。那么只需给定扭向量即可简单算出剩余的4个控制点值了。

 最简单的扭向量计算方法是直接令其为0,那么其等价于Ferguson法(1964)构造的曲面。其缺陷是曲面在插值点处较“平”,这种情况在真实感图形绘制时会变得尤为明显。还有其它一些扭向量的估计方法,如标准公式,但其仅限于在控制点等间距的情况下使用。当控制点间距不相等时,我们建议参考Buss的《3D computer graphics》中提到的方法。当然,作者推荐的《Curves and Surfaces for Computer-Aided Geometric Design: A Practical Guide》by Gerald Farin也绝对是几何造型领域的必备导论书了。

Comments

假如《宫》被改编成RPG…

 这周休息,看完了最后十集的《宫》。突然想到,如果本剧能改编为3D RPG for pc就更好了。

Comments