Fractional Programming

First, the general form of fractional programming is that:

$$Minimize   λ = ƒ(x) = \frac{a(x)}{b(x)} (x∈S)$$ $$s.t.  ∀x∈S  b(x)>0$$

      so how do we solve this problem? a effective solution is parametric search: means we assume x is our answer and check whether x can Minimize λ; similar to binary search
      First,we must prove why we can use parametric search:

      suppose λ^*=ƒ(x^*) is the optimal solution,we have

$$λ^* = ƒ(x^*) = \frac{a(x^*)}{b(x^*)}$$ $$λ^* *b^*=a^*$$ $$a^*-λ^* *b^*=0$$

     then we can get a new function: $g(λ)=min(a(x)-λ*b(x))  x∈S$

      Let's talk something  about $g(λ)$

      1. monotonicity：$g(λ)$ is a decreasing monotone function, ∀x, if $λ_1  < λ_2$,then $g(λ_1) > g(λ_2)$
      Prov:  suppose when $λ= λ_1$, $min(a(x)-λ_1*b(x))=a(x_1)-λ_1*b(x_1)$    $$g(λ_1)=min(a(x)-λ_1*b(x))$$ $$=a(x_1)-λ_1*b(x_1)$$ $$>a(x_1)-λ_2*b(x_1)$$ $$≥min(a(x)-λ_2*b(x))=g(λ_2)$$

     the last step shows that the answer $x_1$make $g(λ_1)$Minimize, does not means that $g(λ_2)$Minimize too.
     monotonicity means that we can use binary search, but what is termination ?
     2. Dinkelbach Theorem: suppose $λ^*$ is the optimal solution, then $g(λ)=0$ if and only if $λ=λ^*$
    First we prove necessary condition: $λ^*$ is the optimal solution ,if $λ=λ^*$, \⇒(g(λ)=0\): ∀x∈S, we have

$$λ^*=\frac{a(x^*)}{b(x^*)}≤\frac{a(x)}{b(x)}  ⇒  a(x)-λ^**b(x)≥0$$

but to $x^*$:   $$λ^*=\frac{a(x^*)}{b(x^*)}  ⇒  a(x^*)-λ^**b(x^*)=0$$

so, the minimize value of $g(λ^*)$ is 0, when $x=x^*$
    now we move to sufficient condition:  if $λ^*$ is optimal solution, $g(λ)=0  ⇒  λ=λ^*$ if exist $x_1$  we have $g(λ_1)=min(a(x)-λ_1*b(x))=0$, we will make sure $λ_1=ƒ(x_1)$ is the optimal solution
    we use contradition to prove: if there exist other $λ_2=ƒ(x_2)$,$λ_2<λ_1$

$$λ_2=\frac{a(x_2)}{b(x_2)}<λ_1  ⇒  a(x_2)-λ_1*b(x_2)<0$$

    so we get a contradition .

    we can get conclusion now

    core idea of this solution is to transform fractional programming to non-fractional programming, 0-1 fractional programming is special case of fractional programming. 

    Two example about 0-1 fractional programming:
    maximum density subgraph,  the optimum ratio spanning tree

ripple simulation

        Simulate ripple of raindrop Using GLSL shader language to speedup, when ripple  calculation completed. compress it into 1-Dimension texture for rendering

Data structure:

sky:  sky box


wavefild:  water surface


wavesource: raindrop representation


waveEvaluator: height field. represent as 1-dimension texture, used wavesource as input, each pixel value is the superposition of each raindrop in this pixel
the equation of ripple as follow:

$W(t,d) = \frac{A}{1+λ_1*t^2+λ_2*d^2}cos(\frac{Ω}{1+λ_3*t}*(t+\frac{d}{v})+π)$

where :
         t:     time
         d:     distance
        λ_1:   coefficient of space
        λ_2:   coefficient of time
        λ_3:   coefficient of frequency
          A:    amplitude
          v:    velocity
          Ω:   initialization of frequency


wavefielddraw: for scene rendering, included sky, water surface, height fileld. using GPU speedup


GPU pipeline

result:

No raindrop

Raindrop number: 100

Raindrop number: 500

Appendix:

GLSL Shader language for computing height field

fragment shader file:
void main()
{
vec2 cpos = gl_TexCoord[0].xy * poolSize;
int i;
float Amplitude = 0.0;
for (i=0; i<n-1; i++)
{
int base = i * SourceSize;
vec4 v1 = texture1D(source, (base * InvN)*0.9998+0.0001);
vec4 v2 = texture1D(source, ((base + 1) * InvN)*0.9998+0.0001);
float InitT = v1.x;
float PosX = v1.y;
float PosY = v1.z;
float A = v1.w;
float Lambda = v2.x;
float Omega = v2.y;
float Velocity = v2.z;
float MaxPhase  = v2.w;
float dTime = time - InitT;
float dist = distance(cpos, v1.yz);
float x = dTime - dist / Velocity;
float phase = Omega/(1.0+Lambda2*dTime) * x - PI/2;
Amplitude += A /(1+Lambda*dist*dist+Lambda1*dTime*dTime) *
                             (cos(phase)) * (x>0?1:0)*(phase<MaxPhase?1:0);    
}
gl_FragData[0] = vec4(Amplitude,Amplitude,Amplitude,1.0) ;
}

vertex shader file:
void main()
{
gl_Position = gl_ModelViewProjectionMatrix * gl_Vertex;
gl_TexCoord[0] = gl_MultiTexCoord0;
}

GLSL Shader language for Rendering

fragment shader file:

void main()
{
vec3 lp = lightPosition;
        // vec3 eyedir=(eye-pos).xyz;
vec3 eyedir = (pos-eye).xyz;
vec3 n = normalize(normal.xyz);
vec3 I = normalize(lp - pos.xyz);
vec3 R = normalize(-reflect(-normalize(eyedir), n));
vec3 color = LightColor * max(0.0,dot(I, n)) +
LightSpecularColor * pow(max(0.0,dot(R,I)),4.0);
gl_FragColor = vec4(color,1.0);
}


vertex shader file:

void main()
{
vec4 h = texture2D(height, gl_MultiTexCoord0.xy);
vec4 h1 = texture2D(height, gl_MultiTexCoord0.xy+vec2(dx,0));
vec4 h2 = texture2D(height, gl_MultiTexCoord0.xy+vec2(-dx,0));
vec4 h3 = texture2D(height, gl_MultiTexCoord0.xy+vec2(0,dy));
vec4 h4 = texture2D(height, gl_MultiTexCoord0.xy+vec2(0,-dy));
pos = vec4(gl_Vertex.x, h.y, gl_Vertex.z, 1);
vec3 pos1 = vec3(gl_Vertex.x+dpx, h1.y, gl_Vertex.z);
vec3 pos2 = vec3(gl_Vertex.x-dpx, h2.y, gl_Vertex.z);
vec3 pos3 = vec3(gl_Vertex.x, h3.y, gl_Vertex.z+dpy);
vec3 pos4 = vec3(gl_Vertex.x, h4.y, gl_Vertex.z-dpy);

vec3 n1 = cross(pos1-pos.xyz, pos3-pos.xyz);
vec3 n3 = cross(pos2-pos.xyz, pos4-pos.xyz);

normal = -vec4(0.5*(normalize(n1)+normalize(n3)),1);
eye = gl_ModelViewMatrixInverse * vec4(0,0,0,1);
pos = gl_ModelViewMatrix * pos;
gl_Position = gl_ProjectionMatrix * pos;
gl_TexCoord[0] = gl_MultiTexCoord0;
}

Chebyshev distance and Manhattan distance

        Today ，when I solve the problem of Meeting Point, I learned new knowledge about Chebyshev distance ,we can check its definition from wiki http://en.wikipedia.org/wiki/Chebyshev_distance 
     
        Chebyshev distance is a metric defined on a vector space where the distance between two vector is the greatest of their differences along any coordinate dimension. for example, in two dimension, point(x,y) to its 8 adjacent points are 1(see images below, left is represent chebyshev distance, right is Manhattan distance)

chebyshev distance

Manhattan distance

   
       The relationship between Chebyshev distance (also know as Linf) and Manhattan distance(also know as L1) is subtle, in one dimension , they are equivalent; in two dimension, all points whose Chebyshev distance equal to r with point(x,y) become a circles in the form of square, with side of length 2r, very interesting, all points whose Manhattan distance equal  to r with point(x,y) also become a circles in the form of square, with side of length sqrt(2)*r , we can rotate an scale the axis to transform Chebyshev distance to manhattan distance.
       But  what's more, the equivalence is not generalize to high dimension

       so what is the advantage of Manhattan distance?  independent, more precisely, Manhattan distance calculation is independent between dimensions. 
    
      now you can finish the problem meeting point , of course , if you want to pass system test, you must implement tricky.
   
    

Python编程技巧小结#1

        学习Python有一段时间了，发觉Python的代码风格很飘逸，与C++的严谨风格不同，Python是一种更灵活简练的语言。



1.  循环式用xrange代替range
           我发觉在很多有关Python的书中都在循环中使用range，但很少涉及到xrange，其实
     如果循环的次数很大的话，xrange的效率要要明显高于range,毕竟xrange的实现不需要维
     护一个列表


2.  Python的三元表达式

            Python在诞生的时候的很长时间里，都没有加入类似于C++的三元操作符（?:）,以至
     于很多人都还认为Python中没有实现三元操作符，其实在Python 2.5之后，Python之父
     Guido van Rossum加入了这个新的特性：正确的表达是：V1 if C else V2  （详见python核
     心编程）

       其实还有很多方法可以模拟，我就很喜欢通过列表来模拟：

       C++ :  cout<<(n>=0)?"Positive":"Negative"<<endl;

       print "Positive" if n>=0 else "Negative"   # 标准的Python三元符

       print ["Positive","Negative"][n>=0]   # 列表模拟

3.  交换两个元素

            在C++中交换两个元素我们往往采用下面的代码

                   Type tmp=x ; x=y ; y=tmp;

             Python当然也可以这样实现，但这样显然没有很好利用Python的特性利用Python的 
      元祖赋值，我们可以直接写成下面更简洁的代码

                    x,y = y,x (够简洁吧)

4.   Python的变量是没有类型的，它可以与任意类型的内存对象绑定，这也是Python与大多
      数语言的一个重要区别

      如 ：a=1          #这时a与一个数值型绑定

             a="Test"     #a变为与一个字符型绑定

   

5.   Python的可变（mutable）类型与不可变（immutable）类型

           这个恐怕是Python初学者感到最头痛的问题之一了，Python把所有的类型划分为两 
      种：可变与不可变：从字面意思很好理解,可变类型就是指值可变的类型，不可变类型是
      指值不可变的类型（就像C++的常量），可变数据类型有列表，字典，不可变的数据类型
      有元祖，数值型和字符串，但慢慢深究，就会发现，其实里面大有奥秘

      首先看一看这个例子：

              a=1

              a=2

              这个运行不会出错，若按照前面所讲，a是一个数值型，应该是不可变的，那么这
       里为什么没有错呢？其实上面的a与下面的a已经不是同一个对象了，可以通过id查看内
       存对比更准确的说法是a抛弃旧的对象（1），而与新的对象（“Test”）绑定，原对象
      （1）并没有改变

            赋值：不管是可变还是不可变类型，变量赋值都是把两个变量与同一个内存绑定

             a=1

             b=a

             a is b   #结果为True

             a=[1,2,3]

             b=a

             a is b   #结果也为True 

            参数传递：与C++的参数传递分为值传递和引用传递不同，Python没有区分，统一进
            行赋值操作，也就是实参和形参都绑定到统一个对象中

            下面来看两个例子：
            def fun(a):

                    a=2

            data=1

             fun(data)

             print data  

             # 输出结果为1，也就是形参的改变并没有对实参产生影响，我们再看

             def fun(a):

                      a[0]=2

             data=[1,2,3]

             fun(data)

              print data  

              # 输出结果为[2,2,3]
               这里又该变了，那是否矛盾呢？当然不矛盾，对可变类型，形参与实参都绑定到同
         一类型中，修改形参值当然会改变实参，但对不可变类型来说，正如前面说的：改变
         形参其实是形参抛弃旧的对象，而与新的对象绑定
              
      因此，要将可变类型进行参数传递，尤其要小心防止被修改，若不想修改应该传递副本


6.  lambda的谨慎使用
             lambda创建一个匿名函数，通常与map，reduce ，filter等联合使用，但只有当函
      数很短（例如只有一个表达式时）才选择，不能有循环等语句，否则，应直接定义一个
      函数

        a=range(1,11)
        b = filter( lambda x:x%3==0 , a)
        b=[3,6,9]   #返回能被3整除的数
        总之，使用lambda只是使代码变得简洁，而不会带来性能上的提高

      总之，Python是一种灵活高效的语言，选择Python，就像以前曾经看过的一篇博客所说的，就要使Python有一种“一字值千金”的感觉

challenge solution #1

   这段时间写论文实在无聊，再做了两题：分别是Grid walking和Permutation game
 
   Grid walking题意是指在N维空间中，每一维都有一个范围，分别是[1,di]，给定一个在这N维空间中的初始位置（x1,x2,x3....xn），每一步可以选择任意一维，然后再这维上向前或向后移动一步，问走M步并且不越界（即1<=xi<=di）有多少种方式？很明显的DP问题，只要注意到每一维相互独立的，每一层单独考虑，那么转移方程不难列出，但我在实现时，却出了一个极为NC的错误，优先级没有考虑（囧），害我交了四次WA。

    Permutation game是一个博弈论问题，由于N不是很大，全部状态只有2^15,可以通过位运算，列出每一点的SG函数值，就可以知道必胜态和必败态了。

    

    争取这个星期做完20题。