# some default values: 
pars = c(0.1, 0.12, 0.01); N = 1000;

if(.Platform$OS.type == 'windows'){
    #create R menu-items 
    winMenuAdd("Diffusion Model");
    winMenuAddItem("Diffusion Model", "Exact Mean and Variance",
                   "exactmeanvar(pars, menucall=T)");

    winMenuAddItem("Diffusion Model", "Pdf Mean and Variance", "pdfmeanvar(pars, menucall=T)");

    winMenuAddItem("Diffusion Model", "Simulated Mean and Variance", "simmeanvar(pars, N, menucall=T)"); 

    winMenuAddItem("Diffusion Model", "-", "")

    winMenuAddItem("Diffusion Model", "Help Diffusion Model", "helpmeanvar()");
} else {
    cat("The variable pars has the default value\npars = c(driftRate = ", pars[1],", boundarySeparation = ",pars[2],", diffusionVariance = ",pars[3],")\n\n")
    cat("Use","exactmeanvar(pars)", "to compute", "Exact Mean and Variance\n")
    cat("Use","pdfmeanvar(pars)", "to compute", "Pdf Mean and Variance\n")
    cat("Use","simmeanvar(pars)", "to compute", "Simulated Mean and Variance\n")
    cat("---------------------------------------------------------------\n")
    cat("Use","helpmeanvar()", "for", "Help Diffusion Model\n")
}

# the functions: 
# exactmeanvar 
# pdfmeanvar 
# simmeanvar 
# helpmeanvar

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
exactmeanvar = function(pars, menucall=F) {
    if (menucall)
    {
        pars = as.numeric(strsplit(winDialogString("Enter drift rate v, boundary separation a, and diffusion variance s2, seperated by semicolons:", '0.1; 0.12; 0.01'), ';')[[1]])
    }
    v  = pars[1];
    a  = pars[2];
    s2 = pars[3];

    if (a < 0)
    {
        cat("\n\n Boundary separation has to be greater than zero!\n\n");
        return();
    }
    if (s2 < 0)
    {
        cat("\n\n Diffusion variance has to be greater than zero!\n\n");
        return();
    }

    if (v == 0) #this requires different equation!
    {
        meanrtime = a^2 / (4*s2);
        varrtime  = a^4 / (24*(s2^2));
        sdrtime   = sqrt(varrtime);

        return(meanrtime, sdrtime, varrtime);
    }

    # calculate mean:
    d = (v*a)/s2;
    meanrtime = ( -(1/2)*a*(exp(-d)-1) ) / (v*(exp(-d)+1));
    # calculate variance:
    varrtime  = -(1/2)*a*((s2*exp(-2*d)) + (2*exp(-d)*a*v) - s2) * ( 1/(v^3*(exp(-2*d)+2*exp(-d)+1)) );

    sdrtime = sqrt(varrtime);

    return(list(meanrtime, sdrtime, varrtime));
} 
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

pdfmeanvar = function(pars, menucall=F) {
#==========================================================================
# It is more efficient to program loops such as those occurring in this
# function in C, and then link the C code to R. Because this function is
# mainly for illustrative purposes, I have not done so here.
#==========================================================================
	if (menucall)
    {
        pars = as.numeric(strsplit(winDialogString("Enter drift rate v, boundary separation a, and diffusion variance s2, seperated by semicolons:", '0.1; 0.12; 0.01'), ';')[[1]])
    }
    v  = pars[1];
    a  = pars[2];
    s2 = pars[3];

    if (a < 0)
    {
        cat("\n\n Boundary separation has to be greater than zero!\n\n");
        return();
    }
    if (s2 < 0)
    {
        cat("\n\n Diffusion variance has to be greater than zero!\n\n");
        return();
    }
    # Set up dummy variables to prevent superfluous calculation in loops:
    dum1 = (pi*s2/a^2) * exp(v*0.5*a/s2);
    dum2 = 0.5*(v^2) / s2;
    dum3 = -0.5*(pi^2)*s2 / (a^2);

    infsum = array();
    prob   = array();
    Ex     = array();
    Ex2    = array();
    ok     = 0;
    t      = 0;
    i      = 1;
    while (ok == 0) #infinite integral over t
    {
        t      = t + 0.01; # t + 0.001 will evaluate once every ms.
        part1  = dum1 / exp(dum2*t);
        ok2    = 0;
        k      = 1;
        totsum = 0;
        while (ok2 == 0) #infinite sum over k
        {
            infsum[k] = k*exp(dum3*(k^2)*t)*sin(0.5*k*pi);
            totsum    = sum(infsum);
            tolerance = (10^-29)*totsum;
            # now check whether 2 consecutive values are less than the tolerance level:
            if ( (k > 100) && (infsum[k] < tolerance) && (infsum[k-1] < tolerance) )
            {
                ok2 = 1;
            }
            k = k + 1;
        }
        prob[i] = part1 * totsum;

        if (prob[i] < 0) #this should never happen, but may happen because of rounding
        {
            print("prob < 0!");
            flush.console();
            prob[i] = 0;
            #return();
        }

        Ex[i]   = t*prob[i];
        Ex2[i]  = (t^2)*prob[i];

        pdfsum  = sum(prob);
        tolerance = (10^-29)*pdfsum; #this method may not work for bimodal distributions

        if ( (prob[i] < tolerance) )
            ok = 1;

        i = i + 1;
    }

    meanrtime = sum(Ex)/pdfsum;
    varrtime  = (sum(Ex2)/pdfsum) - (meanrtime)^2;
    sdrtime   = sqrt(varrtime);

    return(list(meanrtime, sdrtime, varrtime));
} 
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
simmeanvar = function(pars, N, menucall=F) {
#========================================================================== 
# generate data from model by random walk approximation (Ratcliff & 
# Tuerlinckx, 2002, pp. 441-442) 
# function takes as parameters: v, a, s2, N (last one is number of draws)
# It is more efficient to program loops such as those occurring in this
# function in C, and then link the C code to R. Because this function is
# mainly for illustrative purposes, I have not done so here.
#==========================================================================
    if (menucall)
    {
        pars = as.numeric(strsplit(winDialogString("Enter drift rate v, boundary separation a, diffusion variance s2, seperated by semicolons:", '0.1; 0.12; 0.01'), ';')[[1]]);
        N    = as.numeric(winDialogString("Enter number of simulated values:", "1000"));
    }
    v  = pars[1];
    a  = pars[2];
    s  = sqrt(pars[3]);

    if (a < 0)
    {
        cat("\n\n Boundary separation has to be greater than zero!\n\n");
        return();
    }
    if (s < 0)
    {
        cat("\n\n Diffusion variance has to be greater than zero!\n\n");
        return();
    }

    h           = 0.05 * (1/1000); #gives stepsize h in ms
    delta       = s*sqrt(h);       #distance up or down
    pstepdown   = 0.5*(1 - (v*sqrt(h)/s)); # N.B. Ratcliff and Tuerlinckx (2002) mistakenly give sqrt(h/s)!
    j           = 1;
    rtime       = 0;
    as.vector(rtime); #initialize rtime vector
    for (i in 1:N)
    {
        position = 0.5*a;
        hit      = 0;
        RT       = 0;
        while (hit == 0)
        {
            if (runif(1) < pstepdown)
            {
                position = position - delta; # step down
            }
            else
            {
                position = position + delta; # step up
            }
            if ((position >= a) || (position <= 0))
            {
                rtime[j] = RT;
                j        = j + 1;
                hit      = 1;
            }
            RT = RT + h;
        }
    cat("Sim nr. ", i, " out of ", N, " completed \n");
    flush.console();
    }
    meanrtime = mean(rtime);
    sdrtime   = sd(rtime);
    varrtime  = var(rtime);
    return(list(meanrtime, sdrtime, varrtime));
}

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

helpmeanvar = function() {
    winDialog('ok', 'Select the help file called <help.txt>');
    filenaam = file.choose();
    file.show(filenaam);
} 

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++