Commit d8407280 authored by Johannes Blaschke's avatar Johannes Blaschke
Browse files

start working on rejection sampling

parent fecf6229

1d_diffusion.ipynb

+14 −14
Original line number Diff line number Diff line 
%% Cell type:code id: tags:



``` julia

include("src/fp.jl")

```



%% Output



    Main.FP



%% Cell type:code id: tags:



``` julia

using .FP

```



%% Cell type:code id: tags:



``` julia

N = 10

D = 1.;

h = 1e-3;

k = 1.;

```



%% Cell type:code id: tags:



``` julia

n = zeros(N);

f = FluxArray{Float64}(N);

```



%% Cell type:code id: tags:



``` julia

n[Integer(length(n)/2)] = 1

```



%% Output



    1



%% Cell type:code id: tags:



``` julia

fill_fr!(f);

fill_lr!(f, n);

```



%% Cell type:code id: tags:



``` julia

f.left

```



%% Output



    10-element Array{Float64,1}:

     0.0

     0.0

     0.0

     0.0

     0.4618638805059607

     0.13816584763737239

     0.0

     0.0

     0.0

     0.0

     0.0



%% Cell type:code id: tags:



``` julia

f.right

```



%% Output



    10-element Array{Float64,1}:

     0.0

     0.0

     0.0

     0.0

     0.5381361194940393

     0.8618341523626276

     0.0

     0.0

     0.0

     0.0

     0.0



%% Cell type:code id: tags:



``` julia

f.left .+ f.right

```



%% Output



    10-element Array{Float64,1}:

     0.0

     0.0

     0.0

     0.0

     1.0

     0.0

     0.0

     0.0

     0.0

     0.0



%% Cell type:code id: tags:



``` julia

fill_df!(f, D, h);

```



%% Cell type:code id: tags:



``` julia

f.diffusive

```



%% Output



    11-element Array{Float64,1}:

        0.0

        0.0

        0.0

        0.0

      461.8638805059607

     -538.1361194940392

      138.1658476373724

     -861.8341523626276

        0.0

        0.0

        0.0

        0.0

        0.0



%% Cell type:code id: tags:



``` julia

fill_sf!(f, D, h, k);

```



%% Cell type:code id: tags:



``` julia

f.stochastic

```



%% Output



    11-element Array{Float64,1}:

     -0.5929904221281879

     -0.0

     -0.42054232430683963

     -0.0

      0.0

     -4.280339552006784

      2.185681116243459

      0.0

     -0.0

      0.9582299402193841

      9.427963410111976

      0.0

      0.0

      0.0

     -0.0

     -0.31415706471809324

      0.6927609681908308



%% Cell type:code id: tags:



``` julia

apply_bc!(f);

```



%% Cell type:code id: tags:



``` julia

f.stochastic

```



%% Output



    11-element Array{Float64,1}:

      0.0

     -0.0

     -0.0

      0.0

     -4.280339552006784

      2.185681116243459

      0.0

     -0.0

      0.9582299402193841

      9.427963410111976

      0.0

      0.0

      0.0

     -0.0

      0.0



%% Cell type:code id: tags:



``` julia

```

src/fp.jl

+28 −0
Original line number Diff line number Diff line
@@ -6,6 +6,10 @@ import Base.size 
import Random.rand!

import Distributions.Normal



include("sample.jl")

import .Sample





"""

FluxArrays store flux data as Struct-of-Arrays. These arrays are allocated and

initialized (to zero) by the constructor (each array has length `n+1` where `n`
@@ -103,8 +107,12 @@ normal random numbers generated by `distribution`. `D` is the diffusion 
coefficient, `h` is the spatial discretization, and `k` is unknown (TODO).

"""

function fill_sf!(flx::FluxArray{T}, D::T, h::T, k::T) where T

    # Fill random weights

    rand!(flx.distribution, flx.stochastic);



    for i = 2:size_cc(flx)

        # Inner part: average over cell, ignore factor of 2

        # Outer part: set correct variance, 2 cancels

        flx.stochastic[i] *= sqrt(D*(flx.left[i] .+ flx.right[i-1])/(k*sqrt(h)));

        #    TODO: Ask Anthony about k, and if this is a nested sqrt ^^^^^^^^^^^

    end
@@ -123,6 +131,26 @@ function apply_bc!(flx::FluxArray{T}, 
end





function fill_ft!(flx::FluxArray{T}, h::T, k::T) where T

    #ratio of time to space step to propogate soln

    nu::T = k/h;



    # Compute total flux

    Ftot::Array{T} = flx.diffusive .+ flx.stochastic;

    Q::Array{T}    = -Ftot*nu;     #take time step to get real fluxes 

    

    for i=2:size_cc(flx)

        # Total particle number

        P::T  = flx.right[i-1] + flx.left[i]

        # Number in right box (i.e. number of particles going left)

        x0::T = flx.left[i];

        # Valid range of fluxes

        lower::T =  -x0;

        upper::T = P-x0;

    end

end





#function getFlux(n, h, N, k, D, FP, time)

#    #get the diffusive flux for FV solver

#

src/sample.jl

0 → 100644
+42 −0
Original line number Diff line number Diff line 
module Sample



export Proposal, rejection_sample



import Random.rand





struct Proposal{T}

    p::T;

    valid::Bool;

end





"""

Sample the given density using rejection sampling. The algorithm will attempt

`trials` many times. Returns a Proposal object.

"""

function rejection_sample(x::Array{T}, p::Array{T}, P::T, trials::Integer=1000) where T



    # Parameter for rejection

    M::T = max.(p);



    # Perform trials

    for trial = 1:trials

        # Using a uniform distribution proposal

        proposal::T = P * Random.rand();

        

        # Accept proposal?

        index::Integer = findmin(abs.(x.-proposal));

        F::T = p(index);

        U::T = Random.rand();

        G::T = T(1); # Ensure that 1 is of type T

        if U < F/(M*G)

            return Proposal(proposal, true);

        end

    end

    

    # Default return => if no sample is generated

    return Proposal(T(0), false);

end



end
 
 No newline at end of file