I am creating a path tracer and I got some strange results when calculating the strength of light.

I am using a rendering equation inspired by this diagram:
[![Rendering equation][1]][1]

1) When I use the cosine everything gets really dark and the edges of my spheres are almost black. That makes sense because I am almost multiplying by cos(PI/2) = 0. I have read that it suppose to cancel out with something so I have removed it which help. But what does it cancel out with?

2) I am using Lambertian diffuse BRDF which means that the "Fr" part of the equation should be "albedo / PI". I am also doing uniform hemisphere sampling which gave me a PDF of 1 / PI.
Now when I calculate BRDF / PDF I get "2 * albedo". That does not seem right. This means that the more and more bounces I make the bigger influence the light has which is nonsense. The light should be dimmer from the 3rd bounce than from the 1st bounce. Where did I make a mistake?

**EDIT 1:**
I am including pseudo-code. The cosine is the dot product. And after boun
There are some things I have not included because they do not contribute to the equation. 

    <!-- language: lang-cpp -->
    //Initialization
    Ray ray;
    float3 passthrought = (1,1,1);
    float3 color = (0,0,0);
    
    //Bouncing
    for(int I = 0; i < 4; i++){
        Object obj = ray.intersect(...);

        color += obj.emission * passthrought;

        float3 newDirection = sample_hemisphere();
        float pdf = 1 / (2 * PI);
        float3 brdf = obj.color / PI; 

        if(pdf > EPSILON)
            passthrought *= dot(normal, newDirection) * brdf / pdf;
        else
            passthrought *= 0;

        //Prepare new bounce
        ray.origin = ray.direction * ray.distance + ray.origin;
        ray.direction = newDirection;
    }


    float3 sample_hemisphere(){
        //Tranformation to world space
        float3 w = ray->worldNormal;
	    float3 axis = fabs(w.x) > 0.1f ? (float3)(0.0f, 1.0f, 0.0f) : (float3)(1.0f, 0.0f, 0.0f);
	    float3 u = normalize(cross(axis, w));
	    float3 v = cross(w, u);

        //Sampling from 2 random variables <0, 1>
	    float e1 = rand();
	    float e2 = rand();

	    float s = sqrt(1.0 - e1 * e1);
	    float phi = 2 * M_PI_F * e2;

	    return normalize(u * cos(2 * M_PI_F * e2) * s + v * sin(2 * M_PI_F * e2) * s + w * e1);
    }



  [1]: https://i.sstatic.net/uKxvp.png