Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

a bug when use autodiff: [codegen_llvm.cpp:taichi::lang::CodeGenLLVM::visit@301] from != to #1569

Open
JJHu1993 opened this issue Jul 23, 2020 · 4 comments
Open

a bug when use autodiff: [codegen_llvm.cpp:taichi::lang::CodeGenLLVM::visit@301] from != to #1569

JJHu1993 opened this issue Jul 23, 2020 · 4 comments
Labels
autodiff This issue is related to automatic differentiation system potential bug Something that looks like a bug but not yet confirmed

Comments

@JJHu1993
Copy link

When I use autodiff, an error message appears:
[verify.cpp:taichi::lang::IRVerifier::basic_verify@46] IR broken: stmt 8233 cannot have operand 8226.

when I close advanced_optimization, the error becomes:
[codegen_llvm.cpp:taichi::lang::CodeGenLLVM::visit@301] from != to

the code:

import sys
import os
import taichi as ti
import time
import math
import numpy as np
ti.init(advanced_optimization=False)
# ti.init()
imgSize = 720
screenRes = ti.Vector([imgSize, imgSize])
img = ti.Vector(3, dt=ti.f32, shape=[imgSize,imgSize])
depth = ti.var(dt=ti.f32, shape=[imgSize,imgSize])
gui = ti.GUI('Cloth', res=(imgSize,imgSize))


clothWid  = 4.0
clothHgt  = 4.0
clothResX = 31
clothResY = 31
max_steps = 1024

# pos_pre    = ti.Vector(3, dt=ti.f32, shape=(clothResX+1, clothResY+1), needs_grad=True)
pos        = ti.Vector(3, dt=ti.f32, shape=(max_steps, clothResX+1, clothResY+1), needs_grad=True)
vel        = ti.Vector(3, dt=ti.f32, shape=(max_steps, clothResX+1, clothResY+1), needs_grad=True)
F          = ti.Vector(3, dt=ti.f32, shape=(max_steps, clothResX+1, clothResY+1), needs_grad=True)
J          = ti.Matrix(3, 3, dt=ti.f32, shape=(max_steps, clothResX+1, clothResY+1), needs_grad=True)
loss       = ti.var(dt=ti.f32, shape=(), needs_grad=True)

eye        = ti.Vector(3, dt=ti.f32, shape=())
target     = ti.Vector(3, dt=ti.f32, shape=())
up         = ti.Vector(3, dt=ti.f32, shape=())
gravity    = ti.Vector(3, dt=ti.f32, shape=())
collisionC = ti.Vector(3, dt=ti.f32, shape=())

mass       = ti.var(dt=ti.i32, shape=())
damping    = ti.var(dt=ti.i32, shape=())
pointSize  = ti.var(dt=ti.i32, shape=())
           
deltaT     = ti.var(dt=ti.f32, shape=())
KsStruct   = ti.var(dt=ti.f32, shape=())
KdStruct   = ti.var(dt=ti.f32, shape=())
KsShear    = ti.var(dt=ti.f32, shape=())
KdShear    = ti.var(dt=ti.f32, shape=())
KsBend     = ti.var(dt=ti.f32, shape=())
KdBend     = ti.var(dt=ti.f32, shape=())
           
fov        = ti.var(dt=ti.f32, shape=())
near       = ti.var(dt=ti.f32, shape=())
far        = ti.var(dt=ti.f32, shape=())
collisionR = ti.var(dt=ti.f32, shape=())



@ti.func
def getNextNeighborKs(n):
    ks = 0.0
    if(n<4):
        ks = KsStruct
    else:
        if(n<8):
            ks = KsShear
        if(n<12) :
            ks = KsBend
    return ks
    
@ti.func
def getNextNeighborKd(n):
    kd = 0.0
    if(n<4):
        kd = KdStruct
    else:
        if(n<8):
            kd = KdShear
        else :
            kd = KdBend
    return kd
    
@ti.func
def getNextNeighborX(n):
    cood = 0
    if (n == 0) or (n == 4) or (n == 7):
        cood = 1
    if (n == 2) or (n == 5) or (n == 6):
        cood = -1
    if (n == 8):
        cood =  2
    if (n ==10):
        cood = -2
    return cood

@ti.func
def getNextNeighborY(n):
    cood = 0
    if (n == 1) or (n == 4) or (n == 5):
        cood = -1
    if (n == 3) or (n == 6) or (n == 7):
        cood = 1
    if (n == 9):
        cood = -2
    if (n ==11):
        cood = 2
    return cood
    
@ti.func
def get_length3(v):
    return ti.sqrt(v.x*v.x+v.y*v.y+v.z*v.z)

@ti.func
def get_length2(v):
    return ti.sqrt(v.x*v.x+ v.y*v.y)



@ti.func
def SolveConjugateGradient(A, x, b):

    r = b-A@x
    d = r
    q = ti.Vector([0.0, 0.0, 0.0])
    alpha_new = 0.0
    alpha = 0.0
    beta  = 0.0
    delta_old = 0.0
    delta_new = r.dot(r)
    delta0    = delta_new
    
    for i in range(0,16):
        q = A@d
        alpha = delta_new/d.dot(q)
        x = x + alpha*d
        r = r - alpha*q
        delta_old = delta_new
        delta_new =  r.dot(r)
        beta = delta_new/delta_old
        d = r + beta*d
        i += 1
        if delta_new< 0.0000001:
            break
    return x   
    
@ti.kernel
def init_cloth():
    for i in range(clothResX+1):
        for j in range(clothResY+1):
            pos[0, i, j]       = ti.Vector([clothWid * (i / clothResX) - clothWid / 2.0, 5.0, clothHgt * (j / clothResY)- clothHgt/2.0])
            vel[0, i, j]       = ti.Vector([0.0, 0.0, 0.0])
            F[0, i, j]         = ti.Vector([0.0, 0.0, 0.0])

        
@ti.kernel
def reset_cloth():
    for i in range(clothResX+1):
        for j in range(clothResY+1):
            pos.grad[0, i, j]  = ti.Vector([0.0, 0.0, 0.0])
            vel.grad[0, i, j]       = ti.Vector([0.0, 0.0, 0.0])
            F.grad[0, i, j]         = ti.Vector([0.0, 0.0, 0.0])


@ti.func
def compute_force(t, coord, jacobian):
    p1           = pos[t-1, coord]
    v1           = vel[t-1, coord]
    F[t, coord]     = gravity*mass + vel[t-1, coord]*damping
    E            = ti.Matrix.identity(ti.f32, 3)
    J[t, coord]     = ti.Matrix([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] )

    for k in range(0,12): 
        ks             = getNextNeighborKs(k)
        kd             = getNextNeighborKd(k)
        
        coord_offcet   = ti.Vector([getNextNeighborX(k), getNextNeighborY(k)])
        coord_neigh    = coord + coord_offcet
        

        if (coord_neigh.x >= 0) and (coord_neigh.x <= clothResX) and (coord_neigh.y >= 0) and (coord_neigh.y <= clothResY):
        
            rest_length = get_length2(coord_offcet * ti.Vector([clothWid / clothResX, clothHgt / clothResY]))
            
            p2          = pos[t-1, coord_neigh]
            
            if t<2: t=2
            v2          = (p2 - pos[t-2, coord_neigh]) / deltaT

            deltaP      = p1 - p2
            deltaV      = v1 - v2 
            dist        = get_length3(deltaP)
            
            if jacobian > 0:
                dist2       = dist*dist
                lo_l        = rest_length/dist
                J[t, coord]   += ks*(lo_l * (E - deltaP.outer_product(deltaP) / dist2) - E)
    
            leftTerm    = -ks * (dist-rest_length)
            rightTerm   = kd * (deltaV.dot(deltaP)/dist)
            F[t, coord]   += deltaP.normalized() * (leftTerm + rightTerm) 		 			

@ti.func
def collision(t, coord):
    #collosoin
    if(pos[t, coord].y<0):
        pos[t, coord].y=0
        
    offcet = pos[t, coord] - collisionC
    dist   = get_length3(offcet)
    if(dist < collisionR):  
        delta0         = (collisionR - dist)
        pos[t, coord]    += offcet.normalized() *delta0
#         pos_pre[coord] = pos[coord]
        vel[t, coord]     = ti.Vector([0.0, 0.0, 0.0])
        
    

                   

        
@ti.kernel
def integrator_implicit(t: ti.i32):
    for i in range(clothResX+1):
        for j in range(clothResY+1):
            coord  = ti.Vector([i, j])
            compute_force(t, coord, 1)
            index  = j * (clothResX+1) + i

            if (index != 0) and (index != clothResX):
    #             collision(t, coord)
    #             tmp             = pos[coord]
                M     = ti.Matrix.identity(ti.f32, 3)*mass
                A     = M - deltaT * deltaT* J[t, coord]
                b     = M@vel[t-1, coord] + deltaT*F[t, coord]

                vel[t, coord] = SolveConjugateGradient(A, ti.Vector([0.0, 0.0, 0.0]), b)
                pos[t, coord]  = pos[t-1, coord] + vel[t, coord] * deltaT
    #             pos_pre[coord]  = tmp   
            
@ti.kernel
def compute_loss(t: ti.i32):
#     ti.atomic_add(loss[None], dt * (target_v[t][0] - v[t, head_id][0]))**2
#     target_v = ti.Vector()
#     loss[None] = dt * (target_v[t][0] - v[t, head_id][0])**2
    loss[None] = pos[t, clothResX//2,clothResY//2][0]
    
def forward():
    steps = 1024
    for t in range(1, steps):
        integrator_implicit(t)
    compute_loss(steps-1)
    
def optimize():
    for iter in range(100):
        reset_cloth()
        with ti.Tape(loss):
            forward()
        print("iter:",iter," loss:",loss[None])
        learning_rate = 5
        v[0, 0, clothResY] -= learning_rate * v.grad[0, 0, clothResY]
        v[0, clothResX, clothResY] -= learning_rate * v.grad[0, clothResX, clothResY]
            


 
pointSize  = 2       
eye        = ti.Vector([3.0, 3.0, 3.0])
target     = ti.Vector([0.0, 3.0, 0.0])
up         = ti.Vector([0.0, 1.0, 0.0])
gravity    = ti.Vector([0.0, -0.00981, 0.0])
collisionC = ti.Vector([0.0, 3.0, 0.0])
collisionR = 1.0
mass       = 1.0
deltaT     = 0.05
damping    = -0.0125
fov        = 1.5
near       = 1.0
far        = 1000.0
           
KsStruct   = 50.0
KdStruct   = -0.25
KsShear    = 50.0
KdShear    = -0.25
KsBend     = 50.0
KdBend     = -0.25

mode = 1
# reset_cloth()
init_cloth()
frame = 0

optimize()
@JJHu1993 JJHu1993 added the potential bug Something that looks like a bug but not yet confirmed label Jul 23, 2020
@JJHu1993
Copy link
Author

I change the index of for loop and add ti.static to the index. Now the code works well without autodiff. If I use autodiff, the error still occurs: stmt 35795 cannot have operand 25249

code:

import sys
import os
import taichi as ti
import time
import math
import numpy as np

ti.init(advanced_optimization=False)
# ti.init()

imgSize = 720
screenRes = ti.Vector([imgSize, imgSize])
img = ti.Vector(3, dt=ti.f32, shape=[imgSize,imgSize])
depth = ti.var(dt=ti.f32, shape=[imgSize,imgSize])
gui = ti.GUI('Cloth', res=(imgSize,imgSize))


clothWid  = 4.0
clothHgt  = 4.0
clothResX = 31
clothResY = 31
max_steps = 1024

# pos_pre    = ti.Vector(3, dt=ti.f32, shape=(clothResX+1, clothResY+1), needs_grad=True)
tmp        = ti.Vector(3, dt=ti.f32, shape=(clothResX+1, clothResY+1))
pos        = ti.Vector(3, dt=ti.f32, shape=(max_steps, clothResX+1, clothResY+1), needs_grad=True)
vel        = ti.Vector(3, dt=ti.f32, shape=(max_steps, clothResX+1, clothResY+1), needs_grad=True)
F          = ti.Vector(3, dt=ti.f32, shape=(max_steps, clothResX+1, clothResY+1), needs_grad=True)
J          = ti.Matrix(3, 3, dt=ti.f32, shape=(max_steps, clothResX+1, clothResY+1), needs_grad=True)
loss       = ti.var(dt=ti.f32, shape=(), needs_grad=True)

eye        = ti.Vector(3, dt=ti.f32, shape=())
target     = ti.Vector(3, dt=ti.f32, shape=())
up         = ti.Vector(3, dt=ti.f32, shape=())
gravity    = ti.Vector(3, dt=ti.f32, shape=())
collisionC = ti.Vector(3, dt=ti.f32, shape=())

mass       = ti.var(dt=ti.i32, shape=())
damping    = ti.var(dt=ti.i32, shape=())
pointSize  = ti.var(dt=ti.i32, shape=())
           
deltaT     = ti.var(dt=ti.f32, shape=())
KsStruct   = ti.var(dt=ti.f32, shape=())
KdStruct   = ti.var(dt=ti.f32, shape=())
KsShear    = ti.var(dt=ti.f32, shape=())
KdShear    = ti.var(dt=ti.f32, shape=())
KsBend     = ti.var(dt=ti.f32, shape=())
KdBend     = ti.var(dt=ti.f32, shape=())
           
fov        = ti.var(dt=ti.f32, shape=())
near       = ti.var(dt=ti.f32, shape=())
far        = ti.var(dt=ti.f32, shape=())
collisionR = ti.var(dt=ti.f32, shape=())



@ti.func
def getNextNeighborKs(n):
    ks = 0.0
    if(n<4):
        ks = KsStruct
    else:
        if(n<8):
            ks = KsShear
        if(n<12) :
            ks = KsBend
    return ks
    
@ti.func
def getNextNeighborKd(n):
    kd = 0.0
    if(n<4):
        kd = KdStruct
    else:
        if(n<8):
            kd = KdShear
        else :
            kd = KdBend
    return kd
    
@ti.func
def getNextNeighborX(n):
    cood = 0
    if (n == 0) or (n == 4) or (n == 7):
        cood = 1
    if (n == 2) or (n == 5) or (n == 6):
        cood = -1
    if (n == 8):
        cood =  2
    if (n ==10):
        cood = -2
    return cood

@ti.func
def getNextNeighborY(n):
    cood = 0
    if (n == 1) or (n == 4) or (n == 5):
        cood = -1
    if (n == 3) or (n == 6) or (n == 7):
        cood = 1
    if (n == 9):
        cood = -2
    if (n ==11):
        cood = 2
    return cood
    
@ti.func
def get_length3(v):
    return ti.sqrt(v.x*v.x+v.y*v.y+v.z*v.z)

@ti.func
def get_length2(v):
    return ti.sqrt(v.x*v.x+ v.y*v.y)



@ti.func
def SolveConjugateGradient(A, x, b):

    r = b-A@x
    d = r
    q = ti.Vector([0.0, 0.0, 0.0])
    alpha_new = 0.0
    alpha = 0.0
    beta  = 0.0
    delta_old = 0.0
    delta_new = r.dot(r)
    delta0    = delta_new
    
    for i in ti.static(range(0,16)):
        q = A@d
        alpha = delta_new/d.dot(q)
        x = x + alpha*d
        r = r - alpha*q
        delta_old = delta_new
        delta_new =  r.dot(r)
        beta = delta_new/delta_old
        d = r + beta*d
#         i += 1
        if delta_new< 0.0000001:
            break
    return x   
    
@ti.kernel
def init_cloth():
    for i,j in tmp:
        pos[0, i, j]       = ti.Vector([clothWid * (i / clothResX) - clothWid / 2.0, 5.0, clothHgt * (j / clothResY)- clothHgt/2.0])
        vel[0, i, j]       = ti.Vector([0.0, 0.0, 0.0])
        F[0, i, j]         = ti.Vector([0.0, 0.0, 0.0])

        
@ti.kernel
def reset_cloth():
    for i,j in tmp:
        pos.grad[0, i, j]  = ti.Vector([0.0, 0.0, 0.0])
        vel.grad[0, i, j]       = ti.Vector([0.0, 0.0, 0.0])
        F.grad[0, i, j]         = ti.Vector([0.0, 0.0, 0.0])


@ti.func
def compute_force(t, coord, jacobian):
    p1           = pos[t-1, coord]
    v1           = vel[t-1, coord]
    F[t, coord]     = gravity*mass + vel[t-1, coord]*damping
    E            = ti.Matrix.identity(ti.f32, 3)
    J[t, coord]     = ti.Matrix([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] )

    for k in ti.static(range(0,12)): 
        ks             = getNextNeighborKs(k)
        kd             = getNextNeighborKd(k)
        
        coord_offcet   = ti.Vector([getNextNeighborX(k), getNextNeighborY(k)])
        coord_neigh    = coord + coord_offcet
        

        if (coord_neigh.x >= 0) and (coord_neigh.x <= clothResX) and (coord_neigh.y >= 0) and (coord_neigh.y <= clothResY):
        
            rest_length = get_length2(coord_offcet * ti.Vector([clothWid / clothResX, clothHgt / clothResY]))
            
            p2          = pos[t-1, coord_neigh]
            
            if t<2: t=2
            v2          = (p2 - pos[t-2, coord_neigh]) / deltaT

            deltaP      = p1 - p2
            deltaV      = v1 - v2 
            dist        = get_length3(deltaP)
            
            if jacobian > 0:
                dist2       = dist*dist
                lo_l        = rest_length/dist
                J[t, coord]   += ks*(lo_l * (E - deltaP.outer_product(deltaP) / dist2) - E)
    
            leftTerm    = -ks * (dist-rest_length)
            rightTerm   = kd * (deltaV.dot(deltaP)/dist)
            F[t, coord]   += deltaP.normalized() * (leftTerm + rightTerm) 		 			

@ti.func
def collision(t, coord):
    #collosoin
    if(pos[t, coord].y<0):
        pos[t, coord].y=0
        
    offcet = pos[t, coord] - collisionC
    dist   = get_length3(offcet)
    if(dist < collisionR):  
        delta0         = (collisionR - dist)
        pos[t, coord]    += offcet.normalized() *delta0
#         pos_pre[coord] = pos[coord]
        vel[t, coord]     = ti.Vector([0.0, 0.0, 0.0])
        
    

                   
@ti.kernel
def integrator_explicit(t: ti.i32):
    
    for i,j in tmp:
        coord  = ti.Vector([i, j])
        compute_force(t, coord, 0)
        index  = j * (clothResX+1) + i

        if (index != 0) and (index != clothResX):
    #                 collision(coord)

    #                 tmp             = pos[coord]
            acc = F[t, coord] / mass
            pos[t, coord]  += vel[t, coord] * deltaT
            vel[t, coord]  += acc * deltaT
    #                 pos_pre[coord]  = tmp     

@ti.kernel
def integrator_implicit(t: ti.i32):
    for i,j in tmp:
        coord  = ti.Vector([i, j])
        compute_force(t, coord, 1)
        index  = j * (clothResX+1) + i

        if (index != 0) and (index != clothResX):
    #             collision(t, coord)
    #             tmp             = pos[coord]
            M     = ti.Matrix.identity(ti.f32, 3)*mass
            A     = M - deltaT * deltaT* J[t, coord]
            b     = M@vel[t-1, coord] + deltaT*F[t, coord]

            vel[t, coord] = SolveConjugateGradient(A, ti.Vector([0.0, 0.0, 0.0]), b)
            pos[t, coord]  = pos[t-1, coord] + vel[t, coord] * deltaT
    #             pos_pre[coord]  = tmp   
            
@ti.kernel
def compute_loss(t: ti.i32):
#     ti.atomic_add(loss[None], dt * (target_v[t][0] - v[t, head_id][0]))**2
#     target_v = ti.Vector()
#     loss[None] = dt * (target_v[t][0] - v[t, head_id][0])**2
    loss[None] = pos[t, clothResX//2,clothResY//2][0]
    
def forward():
    steps = 1024
    for t in range(1, steps):
        integrator_implicit(t)
    compute_loss(steps-1)
    
def optimize():
    for iter in range(100):
        reset_cloth()
        with ti.Tape(loss):
            forward()
        print("iter:",iter," loss:",loss[None])
        learning_rate = 5
        v[0, 0, clothResY] -= learning_rate * v.grad[0, 0, clothResY]
        v[0, clothResX, clothResY] -= learning_rate * v.grad[0, clothResX, clothResY]
            


 
pointSize  = 2       
eye        = ti.Vector([3.0, 3.0, 3.0])
target     = ti.Vector([0.0, 3.0, 0.0])
up         = ti.Vector([0.0, 1.0, 0.0])
gravity    = ti.Vector([0.0, -0.00981, 0.0])
collisionC = ti.Vector([0.0, 3.0, 0.0])
collisionR = 1.0
mass       = 1.0
deltaT     = 0.05
damping    = -0.0125
fov        = 1.5
near       = 1.0
far        = 1000.0
           
KsStruct   = 50.0
KdStruct   = -0.25
KsShear    = 50.0
KdShear    = -0.25
KsBend     = 50.0
KdBend     = -0.25

mode = 1
# reset_cloth()
init_cloth()
frame = 0

optimize()

@xumingkuan
Copy link
Contributor

xumingkuan commented Jul 24, 2020 •
edited
Loading

Thanks for reporting this. If the code works well without autodiff, it's very likely that some rules here are violated: https://taichi.readthedocs.io/en/stable/differentiable_programming.html

Or if it's the compiler's bug, could you please provide a minimal reproducible example so that it's easier to debug? Many thanks!

@erizmr erizmr added the autodiff This issue is related to automatic differentiation system label Feb 15, 2022
@tr654992156
Copy link

I meet the same problem. Under the autodiff mode, the code only works well with ti.static. But the ti.static has very high computation cost when the for-loop is unrolled many times. Do you have any solution？

@tr654992156
Copy link

Thanks for reporting this. If the code works well without autodiff, it's very likely that some rules here are violated: https://taichi.readthedocs.io/en/stable/differentiable_programming.html

Or if it's the compiler's bug, could you please provide a minimal reproducible example so that it's easier to debug? Many thanks!

I meet the same problem. Under the autodiff mode, the code only works well with ti.static. But the ti.static has very high computation cost when the for-loop is unrolled many times. Do you have any solution？

@xumingkuan xumingkuan removed their assignment Nov 25, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
autodiff This issue is related to automatic differentiation system potential bug Something that looks like a bug but not yet confirmed
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
4 participants
@JJHu1993 @xumingkuan @tr654992156 @erizmr