text:
using GNSSTools
using Statistics
using ProgressMeter
using PyPlot
pygui(true)

# h_parms = [1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15]./16.3676e6^2  # TCXO
t_length = 1e-3
CN0s = Array(30:0.25:50)
iterations = 1000
T = 1e-3
f_d = 0

# L1 C/A
f_s = 2.048e6
sigtype = define_l1ca_code_type(; sig_freq=1*L1_freq)
replica = definesignal(sigtype, f_s, T)
signal = definesignal(sigtype, f_s, t_length;
                      prn=1, code_start_idx=1000, f_d=f_d, fd_rate=0,
                      CN0=36, include_phase_noise=false,
                      include_thermal_noise=true, include_adc=true,
                      include_carrier=true, receiver_h_parms=h_parms)
l1ca_SNRs = zeros(length(CN0s))
p = Progress(length(CN0s)*iterations, 1, "Processing L1 C/A...")
for i in 1:length(CN0s)
    SNR = 0
    for j in 1:iterations
        definesignal!(signal; CN0=CN0s[i], new_thermal_noise=true,
                      new_phase_noise=false)
        generatesignal!(signal)
        fd_est, n0_est, SNR_est = courseacquisition(signal, replica, 1)
        SNR += (10^(SNR_est/10))^2
        next!(p)
    end
    l1ca_SNRs[i] = 10log10(SNR/iterations)
end

# save("ch4_acquisition_performance.jld", "results", results)

# L%I
f_s = 10*2.048e6
sigtype = definesignaltype(definecodetype(l5i_codes, L5_chipping_rate),
                           L5_freq, "I")
replica = definesignal(sigtype, f_s, T)
signal = definesignal(sigtype, f_s, t_length;
                      prn=1, code_start_idx=1000, f_d=f_d, fd_rate=0,
                      CN0=36, include_phase_noise=false,
                      include_thermal_noise=true, include_adc=true,
                      include_carrier=true, receiver_h_parms=h_parms)
l5_SNRs = zeros(length(CN0s))
p = Progress(length(CN0s)*iterations, 1, "Processing L5...")
for i in 1:length(CN0s)
    SNR = 0
    for j in 1:iterations
        definesignal!(signal; CN0=CN0s[i], new_thermal_noise=true,
                      new_phase_noise=false)
        generatesignal!(signal)
        fd_est, n0_est, SNR_est = courseacquisition(signal, replica, 1)
        SNR += (10^(SNR_est/10))^2
        next!(p)
    end
    l5_SNRs[i] = 10log10(SNR/iterations)
end

fig = figure()
plot(CN0s, l1ca_SNRs, label="L1 C/A")
plot(CN0s, l5_SNRs, label="L5")
xlabel("C/N₀ (dB⋅Hz)")
ylabel("Acquired SNR (dB)")
legend()
savefig("figures/ch2_l1ca_vs_l5_snr_at_2x_sampling_rate.pdf", dpi=300)


text:
using Test, OrdinaryDiffEq

using SparseArrays
using LinearAlgebra, Random
N = 30
AA = sprand(MersenneTwister(12),  N, N, 0.5)
mm = sprand(MersenneTwister(123), N, N, 0.5)
A = DiffEqArrayOperator(AA)
M = DiffEqArrayOperator(mm'mm)
u0 = ones(N)
prob = ODEProblem(ODEFunction(A; mass_matrix=M), u0, (0.0, 1.0))

for alg in [Rosenbrock23(), Rosenbrock23(linsolve=LinSolveFactorize(lu))]
  sol = solve(prob, alg)
  @test sol.destats.njacs == 0
  @test sol.destats.nw == 1
end

ext: jl


text:
function euclideanDistance(point1::Array, point2::Array)

    nVar = length(point1)
    acc = 0
    for i=1:nVar
        acc += (point1[i]-point2[i])^2
    end
    return sqrt(acc)
end

function pairwiseDistance(array1::Array, array2::Array, metric::Function)

    len = size(array1)[1]

    matrix = zeros(len, len)

    for i=1:len
        for j=i+1:len
            matrix[i,j] = euclideanDistance(array1[i,:], array2[j,:])
            matrix[j,i] = matrix[i,j]
        end
    end

    return matrix
end

function getKnearest(row::Array, k::Integer)

    Knearest = collect(1:k)
    len = length(row)

    for i=1:k
        for j=i+1:k
            if row[Knearest[i]] < row[Knearest[j]]
                Knearest[j], Knearest[i] = Knearest[i],  Knearest[j]
            end
        end
    end

    for i=k+1:len
        if row[i] < row[Knearest[1]]
            Knearest[1] = i
            for j=1:k-1
                if row[Knearest[j]] < row[Knearest[j+1]] || break
                    Knearest[j], Knearest[j+1] = Knearest[j+1],  Knearest[j]
                end
            end
        end
    end

    return Knearest
end

function randomIndexSelection(len::Integer, n::Integer, rng::Random.AbstractRNG)

    if n>len
        n = len;
    end

    indexes = Array{UInt32}(undef, n)

    i=1
    while i<=n
        inside = false
        indexes[i] = (rand(rng, UInt32)%len) + 1;
        for j=1:i
            if i!=j
                inside = indexes[i] == indexes[j]
                inside && break
            end
        end
        if !inside
            i+=1
        end
    end
    return indexes
end

function randomIndexSelection2(len::Integer, n::Integer, rng::Random.AbstractRNG)

    if n>len
        n = len;
    end

    indexes = collect(UInt32,1:len)
    selected = Array{UInt32}(undef, n)

    i=1
    for i=1:n
        inside = false
        index = (rand(rng, UInt32)%len) + 1;
        len = len - 1
        selected[i] = indexes[index]
        indexes[index] = indexes[end]
        pop!(indexes)

    end
    return selected
end

function findprevIndexes(predicate::Function, A::Array, i::Integer)
    previous = Array{Integer}(undef, 0)

    while (aux = findprev(predicate, A, i)) != nothing
        push!(previous, aux)
        i -= 1
    end

    return previous
end

function findnextIndexes(predicate::Function, A::Array, i::Integer)
    next = Array{Integer}(undef, 0)

    while (aux = findnext(predicate, A, i)) != nothing
        push!(next, aux)
        i += 1
    end

    return next
end

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
@testset "pentland" begin
    img = generate_surface(SynthSphere(50), 0.5, [0.1, 0.0, 0.9])
    Z = Pentland(slant = 0.1106572211739, tilt = 0.0)(img)
    @test Z[76,76] ≈ 96.35408069226294
    @test isapprox(Z[1,1], 108.14578383093884; atol = 2e-15)
    @test isapprox(Z[151,151], 108.14578383093864; atol = 2e-15)
    @test Z[90,90] ≈ 93.80666569339293
    @test Z[101,36] ≈ 80.97642794759314
    @test isapprox(Z[22, 136], 108.14578383093865; atol = 2e-15)
    @test Z[42,42] ≈ 79.72206437513667
end


text:
#!/usr/bin/julia

###############################################################################
#
# Kpax3 - Bayesian cluster analysis of categorical data
#
# Copyright (c) 2016-2018 Alberto Pessia <alberto.pessia@gmail.com>
#
# Kpax3 is free software: you can redistribute it and/or modify it under the
# terms of the MIT License. Refer to LICENSE.md file for more info.
#
###############################################################################

###############################################################################
#
# CSV FILE SUPPORT
#
###############################################################################
#
# csv support is still in an early stage of development. The dataset can be
# analyzed without problems but the user API is not mature yet.
#
###############################################################################

###############################################################################
#
# INTRODUCTION
#
###############################################################################
#
# Welcome to the Kpax3 tutorial!
#
# In this tutorial you will learn how to:
# - load a generic categorical dataset from a csv file
# - setup the parameters
# - run the algorithm
# - save the results
# - plot the results
#
# If you didn't do it already, copy this file and the provided dataset to a
# location of your choice.
#
# The tutorial has been designed in such a way that you should just copy/paste
# the commands into your julia environment. Alternatively, it is also possible
# to run the script directly from the command line by typing:
#
# julia Kpax3_csv_tutorial.jl
#
# from within the directory where the file is located.
#
###############################################################################

###############################################################################
#
# NOTE FOR WINDOWS USERS
#
###############################################################################
#
# File paths should be written using the slash '/' character instead of the
# usual backslash '\' character. If you want to use the backslash, you have to
# escape it (i.e. use '\\'). Examples of valid file paths:
#
# "C:/Users/me/Documents/file.txt"
# "C:\\Users\\me\\Documents\\file.txt"
#
###############################################################################

###############################################################################
#
# LOAD THE PACKAGE
#
###############################################################################
#
# First of all, we need to load the package "Kpax3". If the package is not
# already installed in your system, please follow the installation instructions
# written in the README.md file.
#
# To load the package, open a julia session and run the command
using Kpax3;

###############################################################################
#
# INPUT/OUTPUT
#
###############################################################################
#
# The data we will use in this tutorial consist of 75 aligned protein genomes
# of ZIKA Virus (Human). Total length of each sequence is 3425 amino acids, but
# only 186 polymorphisms. Once converted to binary format, dataset will consist
# of 75 units and 382 binary variables.
#
# Data was retrieved from NCBI’s Virus Variation Resource database
# https://www.ncbi.nlm.nih.gov/genome/viruses/variation/Zika/
# on 2017-02-10 and aligned with MUSCLE
# http://www.drive5.com/muscle/
# and converted to a comma delimited file (csv) file.
#
# Edit the following variable with the path to the dataset
input_file = "path_to_zika_dataset.csv";

# Edit the following variable with the path to the output file.
# This file will contain intermediate output generated by Kpax3. Note that file
# extension is not necessary.
output_file = "path_to_output_file";

###############################################################################
#
# SETUP PARAMETERS
#
###############################################################################
#
# Kpax3 comes with default parameters values that should work in the majority
# of applications. At this point, you can already start analysing the dataset.
# Uncomment the following line if you do not want to customize your Kpax3 run
# and use default values:
#settings = KSettings(input_file, output_file);

# -----------------------------------------------------------------------------
# If any setting is omitted, its default value will be used.
#
# - Missing values
# First of all, we need to specify which characters are to be considered
# missing data.
#
# If the dataset is contained in a fasta file (genetic data), then the variable
# is called `miss` and the format is a UInt8 vector representing the
# characters:
#   Default value:
#     miss = UInt8['?', '*', '#', '-', 'b', 'j', 'x', 'z']
#
# If the dataset is contained in a csv file (generic categorical data), then
# the variable is called `misscsv` and the format is a String vector
# representing the values:
#  Default value:
#    misscsv = [""]
#
# - Ewens-Pitman distribution for the partition
# alpha is a Float64 variable constrained in (0, 1).
# theta is a Float64 variable constrained in (-alpha, infinity)
#
# Default values:
#   alpha = 0.5
#   theta = -0.25
#
# - Column status prior probabilities
# This is a non-negative Float64 vector of length 3. If it does not sum to 1,
# it is automatically normalized.
# Probabilities contained in this variable represent an initial guess for the
# proportion of sites possessing the following statuses:
# - uninformative (noise)
# - informative but not characteristic for any cluster (weak signal)
# - informative and characteristic for a cluster (strong signal)
#
# Default value:
#   gamma = [0.6; 0.35; 0.05]
#
# - Beta distribution hyper-parameter
# This is a Float64 variable and it's a bit tricky to explain and to properly
# setup. You can do it on your own, of course, but it requires some good
# understanding of the model. Read the Pessia and Corander (2017) paper to
# learn its usage. I suggest leaving it to its default value.
#
# Default value:
#   r = log(0.001) / log(0.95)
#
# - General settings
# Boolean variable 'verbose' to print diagnostics messages during run.
# Integer variable 'verbosestep' to print diagnostics messages every fixed
# number of steps.
#
# Default values:
#   verbose = false
#   verbosestep = 500
#
# - MCMC settings
# Integer variable 'T' to choose the length of the Markov Chain.
# Integer variable 'burnin' to choose the length of the burnin period.
# Integer variable 'tstep' to tune the thinning of the chain, i.e. samples will
# be collected every 'tstep' iterations.
# Float64 vector 'op' of length 3 with proportions of MCMC operators, i.e.
# [split-merge; gibbs; biased random walk]. To save time, we suggest using a
# very long chain with only split-merge and biased random walk. For particular
# datasets, gibbs sampler might be needed to reach convergence.
#
# Default values:
#   T = 100000
#   burnin = 10000
#   tstep = 1
#   op = [0.5; 0.0; 0.5]
#
# - Genetic Algorithm settings
# Integer variable 'popsize' to choose the population size.
# Integer variable 'maxiter' to choose the total number of iterations.
# Integer variable 'maxgap' to choose how many iterations without new
#   solutions we want to wait before stopping the search.
# Float64 variable 'xrate' for crossover rate.
# Float64 variable 'mrate' for mutation rate.
#
# Default values:
#   popsize = 20
#   maxiter = 20000
#   maxgap = 5000
#   xrate = 0.9
#   mrate = 0.005
# -----------------------------------------------------------------------------

# To save time in this tutorial, we will simulate a short chain. For a real
# data analysis, choose a longer chain in order to reduce the effect of
# autocorrelation between observations.
settings = KSettings(input_file,
                     output_file,
                     miss=UInt8['?', '*', '#', '-', 'b', 'j', 'x', 'z'],
                     misscsv=["", "?", "*", "#", "-", "B", "J", "X", "Z"],
                     alpha=0.5,
                     theta=-0.25,
                     gamma=[0.6; 0.35; 0.05],
                     r=log(0.001) / log(0.95),
                     verbose=true,
                     verbosestep=1000,
                     T=50000,
                     burnin=10000,
                     tstep=1,
                     op=[0.5; 0.0; 0.5],
                     popsize=20,
                     maxiter=20000,
                     maxgap=5000,
                     xrate=0.9,
                     mrate=0.005);

###############################################################################
#
# LOAD DATA
#
###############################################################################
#
# Loading data is straightforward.
zika_data = CategoricalData(settings);

###############################################################################
#
# INITIAL PARTITION
#
###############################################################################
#
# We need to start our search, or simulation, somewhere. Kpax3 offers two ways
# to input an initial partition of the units:
# 1) Provide a csv file with a known initial partition, for example a
#    solution from a previous run;
# 2) Heuristically find a good initial partition.
#
# Here, we will opt for the second option by using k-means. Choose a suitable
# initial number of clusters for your dataset.
using Clustering;
k = 10;
initial_partition = kmeans(convert(Array{Float64}, zika_data.data),
                           k).assignments;

# The csv file can be either a single column of integers (cluster indices)
# representing the partition, or a two-column file where the first column
# contains the units ids and the second column the cluster indices.
#initial_partition = "path_to_initial_partition.csv";

###############################################################################
#
# MAP ESTIMATE (OPTIMIZATION BY GENETIC ALGORITHM)
#
###############################################################################
#
# We will start by obtaining a point estimate equivalent to the maximum a
# posteriori (MAP). In general, MAP estimates provide good solutions to the
# problem of clustering genetic sequences. It is the recommended approach
# for huge datasets or if we are only interested in clustering the sequences.
#
# Keep in mind that the solution might be a local optimum, requiring multiple
# runs of the optimization algorithm to increase the chance of finding the
# global optimum.
#
# Fortunately for us, the genetic algorithm implemented by Kpax3 is good enough
# to get very close to the global optimum. Even if the solution is a local
# optimum, it will still be a good approximation to the exact solution.
map_solution = kpax3ga(zika_data, initial_partition, settings);

# Check the value of the log-posterior distribution (plus a constant)
# Compare different solutions by this value and pick the one with the maximum
# value.
using Printf;
@printf("Log-posterior value (plus a constant): %.4f\n", map_solution.logpp);

###############################################################################
#
# POSTERIOR DISTRIBUTION (MCMC)
#
###############################################################################
#
# To obtain a measure of uncertainty regarding the parameters of interest, we
# need to sample realizations of such parameters from the posterior
# distribution. We will do it by simulating a Markov Chain.
#
# NOTE
# This step is going to take a while.
#
# Depending on your dataset, you need to customize chain length and operators
# proportions. Obviously, the longer the chain (number of simulations) the
# better, at the cost of longer simulation times and increased size of hard
# drive memory required. If the dataset is too big for MCMC technique, MAP
# estimate (see previous section) is the recommended approach.
kpax3mcmc(zika_data, initial_partition, settings);

# Compute the estimate that minimize the posterior expected loss
mcmc_solution = kpax3estimate(zika_data, settings);

# If the MAP estimate is a local optimum, it might happen that the MCMC
# estimate will have a higher posterior probability
if mcmc_solution.logpp > map_solution.logpp
  # let's try again starting from a better solution. If no solution is found,
  # map_solution will be equal to mcmc_solution
  map_solution = kpax3ga(zika_data, mcmc_solution.R, settings);
  nothing;
end

###############################################################################
#
# SAVE RESULTS
#
###############################################################################
#
# We will save both the MAP estimate and the state that minimize the posterior
# expected loss. In general, MAP estimate is a better point estimate.
#
# Edit the path to the output files. No file extension needed.
map_output_file = "path_to_MAP_estimate_output_file";
mcmc_output_file = "path_to_MCMC_estimate_output_file";

writeresults(zika_data, map_solution, map_output_file, what=4, verbose=true);
writeresults(zika_data, mcmc_solution, mcmc_output_file, what=4, verbose=true);

# Output description:
# If what = 1:
#   filename_summary.txt        = Text file with the log posterior probability.
#                                 This can be used to compare different
#                                 parallel runs.
#   filename_partition.csv      = Best partition found by the algorithm.
#
# If what = 2 writes all the previous files plus:
#   filename_attributes.csv     = Complete columns classification matrix.
#
# If what = 3 writes all the previous files plus:
#   filename_characteristic.csv = Characteristic sites and their corresponding
#                                 values.
#
# If what = 4 writes all the previous files plus:
#   filename_dataset.txt        = Original dataset with highlighted clusters
#                                 and their corresponding characteristic
#                                 values. It is meaningful only if the data
#                                 values are single latin characters (a-zA-Z).
#
# To immediately check if the estimate is reliable, open
# "filename_characteristic.csv" and count how many characteristic sites have
# been found. If the ratio
# number of characteristic sites / total number of SNPs
# is greater than 0.1 and each cluster has at least a characteristic amino
# acid (the more the better), then the solution can be considered reliable.
#
# Suppose we are not satisfied with the current MAP solution and we think
# that units should be moved around, or a cluster should be split or two
# clusters merged. This could happen when metadata is available, so that we
# can immediately see if something is not right.
#
# After manually editing the partition file, we can load it again and check
# our hypothesis. No need to run the algorithm again. There is a function
# written specifically for this task.
#
# If needed, uncomment the following commands.
#my_partition = "path_to_manually_edited_partition.csv";
#my_solution = optimumstate(zika_data, my_partition, settings);
#@printf("Log-posterior value (plus a constant): %.4f\n", my_solution.logpp);
#
# You can save this solution the same way we did with MAP and MCMC estimates.

###############################################################################
#
# PLOT RESULTS
#
###############################################################################
#
# In the following, we will use the MAP estimate as the reference partition.
#
# Function signatures are shown in case you want to change default values.

# Choose the default plotting backend, such as GR. Install them with
#   ]add Plots GR
using Plots;
gr();

# For high quality figures, you can save all these plots to SVG format with the
# following command:
#   Plots.savefig(plotfunction(plot_parameters, size=(width, height),
#                 "path_to_figure.svg");
# where 'width' and 'height' are Float64 to be set according to your needs and
# 'plotfunction' is the corresponding desired plot.
#
# 'Plots.savefig' chooses file type automatically by the extension.

# Plot the complete dataset, highlighting clusters and polymorphic sites
# function kpax3plotd(x::AminoAcidData,
#                     state::AminoAcidState;
#                     clusterorder=:auto,
#                     clusterlabel=:auto)
kpax3plotd(zika_data, map_solution)

# Load posterior probabilities:
(k, pk) = readposteriork(output_file);
(id, P) = readposteriorP(output_file);
(site, aa, freq, C) = readposteriorC(output_file);

# Let's check how the Markov Chain performed by looking at trace and jump
# plots. We will do it separately for R (partition of the units) and C
# (partition of the columns).
#
# Compute the entropy of each sample and the average distance of entropies at
# different lags.
(entropy_R, avgd_R) = traceR(output_file, maxlag=50);
(entropy_C, avgd_C) = traceC(output_file, maxlag=50);

# function kpax3plottrace(entropy::Vector{Float64};
#                         maxlag=200,
#                         M=20000,
#                         main="")
kpax3plottrace(entropy_R)
kpax3plottrace(entropy_C)

# function kpax3plotdensity(entropy::Vector{Float64};
#                           maxlag=200,
#                           main="")
kpax3plotdensity(entropy_R)
kpax3plotdensity(entropy_C)

# function kpax3plotjump(avgd::Vector{Float64};
#                        main="")
kpax3plotjump(avgd_R)
kpax3plotjump(avgd_C)

# Posterior distribution of the number of clusters
# function kpax3plotk(k::Vector{Int},
#                     pk::Vector{Float64})
kpax3plotk(k, pk)

# Posterior distribution of adjacency matrix
# function kpax3plotp(R::Vector{Int},
#                     P::Matrix{Float64};
#                     clusterorder=:auto,
#                     clusterlabel=:auto)
kpax3plotp(map_solution.R, P)

# Posterior distribution of features types (not a very informative plot in this
# case because of few polymorphic sites compared to the length of the genome)
# function kpax3plotc(site::Vector{Int},
#                     freq::Vector{Float64},
#                     C::Matrix{Float64})
kpax3plotc(site, freq, C)


text:
@testset "Print" begin
    buf = IOBuffer()

    print(buf, dna"")
    @test String(take!(buf)) == ""

    print(buf, dna"ACGTN")
    @test String(take!(buf)) == "ACGTN"

    print(buf, rna"ACGUN")
    @test String(take!(buf)) == "ACGUN"

    print(buf, dna"A"^100)
    @test String(take!(buf)) == "A"^100

    print(buf, dna"A"^100, width=70)
    @test String(take!(buf)) == string("A"^70, '\n', "A"^30)
end

@testset "Show" begin
    buf = IOBuffer()

    show(buf, dna"")
    @test String(take!(buf)) == "< EMPTY SEQUENCE >"

    show(buf, dna"ATCG")
    @test String(take!(buf)) == "ATCG"
end


text:
struct Summary <: Output
  filename::String
  last_infections::Vector{Float32}
  num_infections::Vector{UInt32}
  Summary(filename::AbstractString, num_trajectories::Integer) = new(filename, fill(NaN32, num_trajectories), Vector{UInt32}(undef, num_trajectories))
end

function pushtrajectory!(s::Summary, trajectory_id::Integer, writelock::Base.AbstractLock, state::MocosSim.SimState, ::MocosSim.SimParams, ::DetectionCallback)
  try
    lock(writelock)
    s.last_infections[trajectory_id] = maximum(MocosSim.time, state.forest.inedges)
    s.num_infections[trajectory_id] = count(MocosSim.istransmission, state.forest.inedges)
  finally
    unlock(writelock)
  end
  nothing
end

function aftertrajectories(s::Summary, ::MocosSim.SimParams)
  file = jldopen(s.filename, "w")
  try
    file["num_infections"] = s.num_infections
    file["last_infections"] = s.last_infections
  finally
    close(file)
  end
  nothing
end

max_stars_repo_name: Novirmor/MocosSimLauncher

max_issues_repo_name: Novirmor/MocosSimLauncher

max_forks_repo_name: Novirmor/MocosSimLauncher


text:
using Test

m = test_page_model()
include("../src/components/GDP.jl")

gdp = add_comp!(m, GDP)

gdp[:pop0_initpopulation] = readpagedata(m, "data/pop0_initpopulation.csv")
gdp[:pop_population] = readpagedata(m, "test/validationdata/pop_population.csv")
gdp[:y_year] = Mimi.dim_keys(m.md, :time)
gdp[:y_year_0] = 2008.

p=load_parameters(m)
set_leftover_params!(m,p)

# run model
run(m)

# Generated data
gdp = m[:GDP, :gdp]

# Recorded data
gdp_compare = readpagedata(m, "test/validationdata/gdp.csv")

@test gdp ≈ gdp_compare rtol=100

cons_percap_consumption_0_compare = readpagedata(m, "test/validationdata/cons_percap_consumption_0.csv")
@test m[:GDP, :cons_percap_consumption_0] ≈ cons_percap_consumption_0_compare rtol=1e-2


text:
using Documenter, SmartGameFormat

makedocs(
    modules = [SmartGameFormat],
    clean = false,
    format = :html,
    sitename = "SmartGameFormat.jl",
    authors = "Christof Stocker",
    linkcheck = !("skiplinks" in ARGS),
    pages = Any[
        "Home" => "index.md",
        "Internals" => Any[
            "lexer.md",
            "parser.md",
        ]
    ],
    html_prettyurls = !("local" in ARGS),
)

deploydocs(
    repo = "github.com/Evizero/SmartGameFormat.jl.git",
    target = "build",
    julia = "0.6",
    deps = nothing,
    make = nothing,
)


text:
@muladd function perform_step!(integrator,cache::DRI1ConstantCache,f=integrator.f)
  @unpack a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE = cache
  @unpack t,dt,uprev,u,W,p = integrator

  # define three-point distributed random variables
  dW_scaled = W.dW / sqrt(dt)
  sq3dt = sqrt(3*dt)
  _dW = map(x -> calc_threepoint_random(sq3dt, NORMAL_ONESIX_QUANTILE, x), dW_scaled)
  chi1 = map(x -> (x^2-dt)/2, _dW) # diagonal of Ihat2
  if !(typeof(W.dW) <: Number)
    m = length(W.dW)
    # define two-point distributed random variables
    _dZ = map(x -> calc_twopoint_random(integrator.sqdt, x),  W.dZ)
    Ihat2 = zeros(eltype(W.dZ), m, m) # I^_(k,l)
    for k = 1:m
      for l = 1:m
        if k<l
          Ihat2[k, l] = (_dW[k]*_dW[l]-integrator.sqdt*_dZ[k])/2
        elseif l<k
          Ihat2[k, l] = (_dW[k]*_dW[l]+integrator.sqdt*_dZ[l])/2
        else k==l
          Ihat2[k, k] = chi1[k]
        end
      end
    end
  end

  # compute stage values
  k1 = integrator.f(uprev,p,t)
  g1 = integrator.g(uprev,p,t)

  # H_i^(0), stage 1
  # H01 = uprev
  # H_i^(0), stage 2
  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
    H02 = uprev + a021*k1*dt + b021*g1*_dW
  else
    H02 = uprev + a021*k1*dt + b021*g1.*_dW
  end

  # H_i^(0), stage 3 (requires H_i^(k) stage 2 in general)
  k2 = integrator.f(H02,p,t+c02*dt)
  H03 = uprev + a032*k2*dt + a031*k1*dt
  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
    H03 += b031*g1*_dW
  else
    H03 += b031*g1.*_dW
  end

  k3 = integrator.f(H03,p,t+c03*dt)

  # H_i^(k), stage 1
  # H11 = uprev
  # H_i^(k), stage 2 and 3
  if typeof(W.dW) <: Number
    H12 = uprev + a121*k1*dt + b121*g1*integrator.sqdt
    H13 = uprev + a131*k1*dt + b131*g1*integrator.sqdt
  else
    H12 = [uprev .+ a121*k1*dt for k=1:m]
    H13 = [uprev .+ a131*k1*dt for k=1:m]
    for k=1:m
      if is_diagonal_noise(integrator.sol.prob)
        tmp = zero(integrator.u)
        tmp[k] = g1[k]
        H12[k] += b121*tmp*integrator.sqdt
        H13[k] += b131*tmp*integrator.sqdt
      else
        H12[k] += b121*g1[:,k]*integrator.sqdt
        H13[k] += b131*g1[:,k]*integrator.sqdt
      end
    end
  end

  # H^_i^(k), stage 1
  # H21 = uprev
  # H^_i^(k), stage 2 and 3

  if typeof(W.dW) <: Number
    g2 = integrator.g(H12,p,t+c12*dt)
    g3 = integrator.g(H13,p,t+c13*dt)
    # for m=1:  H22 = uprev
  else
    g2 = [integrator.g(H12[k],p,t+c12*dt) for k=1:m]
    g3 = [integrator.g(H13[k],p,t+c13*dt) for k=1:m]
    H22 = [copy(uprev) for k=1:m]
    H23 = [copy(uprev) for k=1:m]
    # add inbounds for speed if working properly
    for k=1:m
      if is_diagonal_noise(integrator.sol.prob)
        tmp = zero(integrator.u)
        tmp[k] = b221*g1[k] + b222*g2[k][k] + b223*g3[k][k]
        H22[k] += tmp*integrator.sqdt
        tmp[k] = b231*g1[k] + b232*g2[k][k] + b233*g3[k][k]
        H23[k] += tmp*integrator.sqdt
      else
        H22[k] += (b221*g1[:,l]+b222*g2[l][:,l]+b223*g3[l][:,l])*integrator.sqdt
        H23[k] += (b231*g1[:,l]+b232*g2[l][:,l]+b233*g3[l][:,l])*integrator.sqdt
      end
    end
  end

  # add stages together Eq. (3)
  u = uprev + α1*k1*dt + α2*k2*dt + α3*k3*dt

  # add noise
  if typeof(W.dW) <: Number
    # lines 2 and 3
    u += g1*(_dW*beta11) + g2*(_dW*beta12+chi1*beta22/integrator.sqdt) + g3*(_dW*beta13+chi1*beta23/integrator.sqdt)
    # lines 4 and 5 are zero by construction
    # u += g1*(_dW*(beta31+beta32+beta33)+chi1*integrator.sqdt*(beta42+beta43))
  else
    if is_diagonal_noise(integrator.sol.prob)
      u += g1.*_dW*beta11
      for k=1:m
        u[k] += g2[k][k]*_dW[k]*beta12+g3[k][k]*_dW[k]*beta13
        u[k] += g2[k][k]*chi1[k]*beta22/integrator.sqdt + g3[k][k]*chi1[k]*beta23/integrator.sqdt
        u[k] +=  (m-1)*beta31*g1[k]*_dW[k]
        for l=1:m
          if l!=k
            tmpg = integrator.g(H22[l],p,t)
            u[k] += tmpg[k]*(_dW[k]*beta32 +  Ihat2[k,l]*beta42/integrator.sqdt)
            tmpg = integrator.g(H23[l],p,t)
            u[k] += tmpg[k]*(_dW[k]*beta33 + Ihat2[k,l]*beta43/integrator.sqdt)
          end
        end
      end
    else
      # non-diag noise
      for k=1:m
        g1k = @view g1[:,k]
        g2k = @view g2[k][:,k]
        g3k = @view g3[k][:,k]
        @.. u = u + g1k*_dW[k]*(beta11+beta31)
        @.. u = u + g2k*_dW[k]*beta12 + g3k*_dW[k]*beta13
        @.. u = u + g2k*chi1[k]*beta22/integrator.sqdt + g3k*chi1[k]*beta23/integrator.sqdt
        tmpg = integrator.g(H22[k],p,t)
        @.. u = u + tmpgk*_dW[k]*beta32 + tmpgk*integrator.sqdt*beta42
        tmpg = integrator.g(H23[k],p,t)
        @.. u = u + tmpgk*_dW[k]*beta33 + tmpgk*integrator.sqdt*beta43
      end
    end
  end

  if integrator.opts.adaptive && (typeof(W.dW) <: Number || is_diagonal_noise(integrator.sol.prob) || m==1)

    # schemes with lower convergence order
    if c03!=0.0
      # scheme has det. conv. order 3 and we look for det. conv. order 2 scheme
      rat = c02/c03
      δ₁ = integrator.opts.delta*(rat-1)
      δ₂ = -integrator.opts.delta*rat
      uhat = uprev + ((α1+δ₁)*k1 + (α2+integrator.opts.delta)*k2 + (α3+δ₂)*k3)*dt
    else
      # check against EM
      uhat = uprev + k1*dt
    end
    if is_diagonal_noise(integrator.sol.prob)
      uhat += g1 .* _dW
    else
      uhat += g1 * _dW
    end

    resids = calculate_residuals(u-uhat,uprev,u,integrator.opts.abstol,integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(resids, t)
  end

  integrator.u = u

end

@muladd function perform_step!(integrator,cache::DRI1Cache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack _dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg,uhat,tmp,resids = cache
  @unpack a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE = cache.tab

  m = length(W.dW)
  sq3dt = sqrt(3*dt)
  if typeof(W.dW) <: Union{SArray,Number}
    # tbd
    _dW = map(x -> calc_threepoint_random(sq3dt, NORMAL_ONESIX_QUANTILE, x), W.dW / sqrt(dt))
  else
    # define three-point distributed random variables
    sqrtdt = sqrt(dt)
    @.. chi1 = W.dW / sqrtdt
    calc_threepoint_random!(_dW, sq3dt, NORMAL_ONESIX_QUANTILE, chi1)
    map!(x -> (x^2-dt)/2, chi1, _dW)

    # define two-point distributed random variables
    if m > 1
      calc_twopoint_random!(_dZ,integrator.sqdt,W.dZ)
      for k = 1:m
        for l = 1:m
          if k<l
            Ihat2[k, l] = (_dW[k]*_dW[l]-integrator.sqdt*_dZ[k])/2
          elseif l<k
            Ihat2[k, l] = (_dW[k]*_dW[l]+integrator.sqdt*_dZ[l])/2
          else k==l
            Ihat2[k, k] = chi1[k]
          end
        end
      end
    end
  end

  # compute stage values
  integrator.f(k1,uprev,p,t)
  integrator.g(g1,uprev,p,t)

  # H_i^(0), stage 1
  # H01 = uprev
  # H_i^(0), stage 2
  if is_diagonal_noise(integrator.sol.prob)
    @.. H02 = uprev + a021*k1*dt + b021*g1*_dW
  else
    mul!(tmp1,g1,_dW)
    @.. H02 = uprev + dt*a021*k1 + b021*tmp1
  end

  integrator.f(k2,H02,p,t+c02*dt)
  if is_diagonal_noise(integrator.sol.prob)
    @.. H03 = uprev + a032*k2*dt + a031*k1*dt + b031*g1*_dW
  else
    @.. H03 = uprev + a032*k2*dt + a031*k1*dt + b031*tmp1
  end

  integrator.f(k3,H03,p,t+c03*dt)

  # H_i^(k), stages
  # H11 = uprev
  for k=1:m
    if is_diagonal_noise(integrator.sol.prob)
      fill!(tmpg,zero(eltype(integrator.u)))
      tmpg[k] = g1[k]
      @.. H12[k] = uprev + a121*k1*dt + b121*tmpg*integrator.sqdt
      @.. H13[k] = uprev + a131*k1*dt + b131*tmpg*integrator.sqdt
    else
      g1k = @view g1[:,k]
      @.. H12[k] = uprev + a121*k1*dt + b121*g1k*integrator.sqdt
      @.. H13[k] = uprev + a131*k1*dt + b131*g1k*integrator.sqdt
    end
    integrator.g(g2[k],H12[k],p,t+c12*dt)
    integrator.g(g3[k],H13[k],p,t+c13*dt)
  end

  # H^_i^(k), stages (rewritten)
  for k=1:m
    if is_diagonal_noise(integrator.sol.prob)
      fill!(tmpg,zero(eltype(integrator.u)))
      tmpg[k] = b221*g1[k]+b222*g2[k][k]+b223*g3[k][k]
      @.. H22[k] = uprev+tmpg*integrator.sqdt
      tmpg[k] = b231*g1[k]+b232*g2[k][k]+b233*g3[k][k]
      @.. H23[k] = uprev+tmpg*integrator.sqdt
    else
      g1k = @view g1[:,k]
      g2k = @view g2[k][:,k]
      g3k = @view g3[k][:,k]
      @.. H22[k] = uprev + (b221*g1k+b222*g2k+b223*g3k)*integrator.sqdt
      @.. H23[k] = uprev + (b231*g1k+b232*g2k+b233*g3k)*integrator.sqdt
    end
  end

  # add stages together Eq. (3)
  @.. u = uprev + α1*k1*dt + α2*k2*dt + α3*k3*dt

  # add noise
  if is_diagonal_noise(integrator.sol.prob)
    @.. u = u + g1*_dW*beta11
    for k=1:m
      u[k] = u[k] + g2[k][k]*_dW[k]*beta12 + g3[k][k]*_dW[k]*beta13
      u[k] = u[k]  + g2[k][k]*chi1[k]*beta22/integrator.sqdt + g3[k][k]*chi1[k]*beta23/integrator.sqdt
      if m > 1
        u[k] = u[k] + (m-1)*beta31*g1[k]*_dW[k]
        for l=1:m
          if l!=k
            integrator.g(tmpg,H22[l],p,t)
            u[k] = u[k] + tmpg[k]*(_dW[k]*beta32 + Ihat2[k,l]*beta42/integrator.sqdt)
            integrator.g(tmpg,H23[l],p,t)
            u[k] = u[k] + tmpg[k]*(_dW[k]*beta33 + Ihat2[k,l]*beta43/integrator.sqdt)
          end
        end
      end
    end
  else
    for k=1:m
      g1k = @view g1[:,k]
      g2k = @view g2[k][:,k]
      g3k = @view g3[k][:,k]
      tmpgk = @view tmpg[:,k]
      @.. u = u + g1k*_dW[k]*(beta11+(m-1)*beta31)
      @.. u = u + g2k*_dW[k]*beta12 + g3k*_dW[k]*beta13
      @.. u = u + g2k*chi1[k]*beta22/integrator.sqdt + g3k*chi1[k]*beta23/integrator.sqdt
      for l=1:m
        if l!=k
          integrator.g(tmpg,H22[l],p,t)
          @.. u = u + tmpgk*(_dW[k]*beta32 + Ihat2[k,l]*beta42/integrator.sqdt)
          integrator.g(tmpg,H23[l],p,t)
          @.. u = u + tmpgk*(_dW[k]*beta33 + Ihat2[k,l]*beta43/integrator.sqdt)
        end
      end
    end
  end

  if integrator.opts.adaptive && (typeof(W.dW) <: Number || is_diagonal_noise(integrator.sol.prob) || m==1)
    # schemes with lower convergence order
    if c03!=0.0
      # scheme has det. conv. order 3 and we look for det. conv. order 2 scheme
      rat = c02/c03
      δ₁ = integrator.opts.delta*(rat-1)
      δ₂ = -integrator.opts.delta*rat
      @.. uhat = uprev + ((α1+δ₁)*k1 + (α2+integrator.opts.delta)*k2 + (α3+δ₂)*k3)*dt
    else
      # check against EM
      @.. uhat = uprev + k1*dt
    end
    if is_diagonal_noise(integrator.sol.prob)
      @.. uhat = uhat + g1 .* _dW
    else
      @.. uhat = uhat + g1 * _dW
    end
    @.. tmp = u-uhat

    calculate_residuals!(resids, tmp, uprev, uhat, integrator.opts.abstol,
                                integrator.opts.reltol,integrator.opts.internalnorm, t)

    integrator.EEst = integrator.opts.internalnorm(resids, t)
  end
end

@muladd function perform_step!(integrator,cache::DRI1NMCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack _dW, chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,tmp1,tmpg,uhat,tmp,resids = cache
  @unpack a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE = cache.tab

  m = length(W.dW)
  sq3dt = sqrt(3*dt)
  if typeof(W.dW) <: Union{SArray,Number}
    # tbd
    _dW = map(x -> calc_threepoint_random(sq3dt, NORMAL_ONESIX_QUANTILE, x), W.dW / sqrt(dt))
  else
    # define three-point distributed random variables
    sqrtdt = sqrt(dt)
    @.. chi1 = W.dW / sqrtdt
    calc_threepoint_random!(_dW, sq3dt, NORMAL_ONESIX_QUANTILE, chi1)
    map!(x -> (x^2-dt)/2, chi1, _dW)
  end

  # compute stage values
  integrator.f(k1,uprev,p,t)
  integrator.g(g1,uprev,p,t)

  # H_i^(0), stage 1
  # H01 = uprev
  # H_i^(0), stage 2
  if is_diagonal_noise(integrator.sol.prob)
    @.. H02 = uprev + a021*k1*dt + b021*g1*_dW
  else
    mul!(tmp1,g1,_dW)
    @.. H02 = uprev + dt*a021*k1 + b021*tmp1
  end

  integrator.f(k2,H02,p,t+c02*dt)
  if is_diagonal_noise(integrator.sol.prob)
    @.. H03 = uprev + a032*k2*dt + a031*k1*dt + b031*g1*_dW
  else
    @.. H03 = uprev + a032*k2*dt + a031*k1*dt + b031*tmp1
  end

  integrator.f(k3,H03,p,t+c03*dt)

  # H_i^(k), stages
  # H11 = uprev

  if is_diagonal_noise(integrator.sol.prob)
    @.. H12 = uprev + a121*k1*dt + b121*g1*integrator.sqdt
    @.. H13 = uprev + a131*k1*dt + b131*g1*integrator.sqdt
  else
    mul!(tmp1,g1,integrator.sqdt)
    @.. H12 = uprev + a121*k1*dt + b121*tmp1
    @.. H13 = uprev + a131*k1*dt + b131*tmp1
  end

  # H^_i^(k), stages
  integrator.g(g2,H12,p,t+c12*dt)
  integrator.g(g3,H13,p,t+c13*dt)

  # add stages together Eq. (3)
  @.. u = uprev + α1*k1*dt + α2*k2*dt + α3*k3*dt

  # add noise
  if is_diagonal_noise(integrator.sol.prob)
    @.. u = u + (g1*beta11 + g2*beta12 + g3*beta13)*_dW +  g2*chi1*beta22/integrator.sqdt + g3*chi1*beta23/integrator.sqdt

  else
    @.. tmpg = g1*beta11 + g2*beta12 + g3*beta13
    mul!(tmp1, tmpg, _dW)
    @.. u = u + tmp1
    @.. tmpg = g2*beta22/integrator.sqdt + g3*beta23/integrator.sqdt
    mul!(tmp1, tmpg, chi1)
    @.. u = u + tmp1
  end

  if integrator.opts.adaptive && (typeof(W.dW) <: Number || is_diagonal_noise(integrator.sol.prob) || m==1)

    # schemes with lower convergence order
    if c03!=0.0
      # scheme has det. conv. order 3 and we look for det. conv. order 2 scheme
      rat = c02/c03
      δ₁ = integrator.opts.delta*(rat-1)
      δ₂ = -integrator.opts.delta*rat
      @.. uhat = uprev + ((α1+δ₁)*k1 + (α2+integrator.opts.delta)*k2 + (α3+δ₂)*k3)*dt
    else
      # check against EM
      @.. uhat = uprev + k1*dt
    end
    if is_diagonal_noise(integrator.sol.prob)
      @.. uhat = uhat + g1 .* _dW
    else
      @.. uhat = uhat + g1 * _dW
    end
    @.. tmp = u-uhat

    calculate_residuals!(resids, tmp, uprev, uhat, integrator.opts.abstol,
                                integrator.opts.reltol,integrator.opts.internalnorm, t)

    integrator.EEst = integrator.opts.internalnorm(resids, t)
  end
end

# Roessler SRK for first order weak approx
@muladd function perform_step!(integrator,cache::RDI1WMConstantCache,f=integrator.f)
  @unpack a021,b021,α1,α2,c02,beta11,NORMAL_ONESIX_QUANTILE = cache
  @unpack t,dt,uprev,u,W,p = integrator

  # define three-point distributed random variables
  dW_scaled = W.dW / sqrt(dt)
  sq3dt = sqrt(3*dt)
  _dW = map(x -> calc_threepoint_random(sq3dt, NORMAL_ONESIX_QUANTILE, x), dW_scaled)
  chi1 = map(x -> (x^2-dt)/2, _dW) # diagonal of Ihat2
  if !(typeof(W.dW) <: Number)
    m = length(W.dW)
    # define two-point distributed random variables
    _dZ = map(x -> calc_twopoint_random(integrator.sqdt, x),  W.dZ)
    Ihat2 = zeros(eltype(W.dZ), m, m) # I^_(k,l)
    for k = 1:m
      for l = 1:m
        if k<l
          Ihat2[k, l] = (_dW[k]*_dW[l]-integrator.sqdt*_dZ[k])/2
        elseif l<k
          Ihat2[k, l] = (_dW[k]*_dW[l]+integrator.sqdt*_dZ[l])/2
        else k==l
          Ihat2[k, k] = chi1[k]
        end
      end
    end
  end

  # compute stage values
  k1 = integrator.f(uprev,p,t)
  g1 = integrator.g(uprev,p,t)

  # H_i^(0), stage 1
  # H01 = uprev
  # H_i^(0), stage 2
  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
    H02 = uprev + a021*k1*dt + b021*g1*_dW
  else
    H02 = uprev + a021*k1*dt + b021*g1.*_dW
  end
  k2 = integrator.f(H02,p,t+c02*dt)

  # add stages together Eq. (3)
  u = uprev + α1*k1*dt + α2*k2*dt

  # add noise
  if typeof(W.dW) <: Number
    # lines 2 and 3
    u += g1*(_dW*beta11)
  else
      if is_diagonal_noise(integrator.sol.prob)
        u += g1.*_dW*beta11
      else
        # non-diag noise
        for k=1:m
          g1k = @view g1[:,k]
          @.. u = u + g1k*_dW[k]*(beta11)
        end
      end
  end

  integrator.u = u

end

@muladd function perform_step!(integrator,cache::RDI1WMCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack _dW,_dZ,chi1,Ihat2,tab,g1,k1,k2,H02,tmp1 = cache
  @unpack a021,b021,α1,α2,c02,beta11,NORMAL_ONESIX_QUANTILE = cache.tab

  m = length(W.dW)
  sq3dt = sqrt(3*dt)
  if typeof(W.dW) <: Union{SArray,Number}
    # tbd
    _dW = map(x -> calc_threepoint_random(sq3dt, NORMAL_ONESIX_QUANTILE, x), W.dW / sqrt(dt))
  else
    # define three-point distributed random variables
    @.. chi1 = W.dW / sqrt(dt)
    calc_threepoint_random!(_dW, sq3dt, NORMAL_ONESIX_QUANTILE, chi1)
    map!(x -> (x^2-dt)/2, chi1, _dW)

    # define two-point distributed random variables
    calc_twopoint_random!(_dZ,integrator.sqdt,W.dZ)

    for k = 1:m
      for l = 1:m
        if k<l
          Ihat2[k, l] = (_dW[k]*_dW[l]-integrator.sqdt*_dZ[k])/2
        elseif l<k
          Ihat2[k, l] = (_dW[k]*_dW[l]+integrator.sqdt*_dZ[l])/2
        else k==l
          Ihat2[k, k] = chi1[k]
        end
      end
    end
  end

  # compute stage values
  integrator.f(k1,uprev,p,t)
  integrator.g(g1,uprev,p,t)

  # H_i^(0), stage 1
  # H01 = uprev
  # H_i^(0), stage 2
  if is_diagonal_noise(integrator.sol.prob)
    @.. H02 = uprev + a021*k1*dt + b021*g1*_dW
  else
    mul!(tmp1,g1,_dW)
    @.. H02 = uprev + dt*a021*k1 + b021*tmp1
  end

  integrator.f(k2,H02,p,t+c02*dt)

  # add stages together Eq. (3)
  @.. u = uprev + α1*k1*dt + α2*k2*dt

  # add noise
  if is_diagonal_noise(integrator.sol.prob)
      @.. u = u + g1*_dW*beta11
  else
    for k=1:m
      g1k = @view g1[:,k]
      @.. u = u + g1k*_dW[k]*beta11
    end
  end
end

# Stratonovich sense

@muladd function perform_step!(integrator,cache::RSConstantCache,f=integrator.f)
  @unpack a021,a031,a032,a131,a141,b031,b032,b121,b131,b132,b141,b142,b143,b221,b231,b331,b332,b341,b342,α1,α2,α3,α4,c02,c03,c13,c14,beta11,beta12,beta13,beta14,beta22,beta23,NORMAL_ONESIX_QUANTILE = cache
  @unpack t,dt,uprev,u,W,p = integrator

  # define three-point distributed random variables
  dW_scaled = W.dW / sqrt(dt)
  sq3dt = sqrt(3*dt)
  _dW = map(x -> calc_threepoint_random(sq3dt, NORMAL_ONESIX_QUANTILE, x), dW_scaled)
  if !(typeof(W.dW) <: Number)
    m = length(W.dW)
    # define two-point distributed random variables
    _dZ = map(x -> calc_twopoint_random(integrator.sqdt, x),  W.dZ)
    Ihat2 = zeros(eltype(W.dZ), m, m) # I^_(k,l)
    for k = 1:m
      for l = 1:k-1
        Ihat2[k, l] = _dW[k]*_dZ[l]
        Ihat2[l, k] = -_dW[k]*_dZ[l]
      end
    end
  end
  # compute stage values
  k1 = integrator.f(uprev,p,t)
  g1 = integrator.g(uprev,p,t)

  # H_1^(0)
  # H01 = uprev
  # H_2^(0)
  H02 = uprev + a021*k1*dt

  # H_1^(k)
  # H11 = uprev
  # H_2^(k), stage 2
  if typeof(W.dW) <: Number
    H12 = uprev + b121*g1*_dW
  elseif is_diagonal_noise(integrator.sol.prob)
    H12 = Vector{typeof(uprev)}[uprev .+ b121*g1[k]*_dW[k] for k=1:m]
  else
    H12 = Vector{typeof(uprev)}[uprev .+ b121*g1[:,k]*_dW[k] for k=1:m]
  end

  if typeof(W.dW) <: Number
    g2 = integrator.g(H12,p,t)
  else
    g2 = [integrator.g(H12[k],p,t) for k=1:m]
  end

  # H_3^(k)
  if typeof(W.dW) <: Number
    H13 = uprev + a131*k1*dt + b131*g1*_dW + b132*g2*_dW
  else
    H13 = [zero(typeof(uprev)) for k=1:m]
    for k=1:m
      H13[k] += uprev .+ a131*k1*dt
      if is_diagonal_noise(integrator.sol.prob)
        H13[k] += b131*g1[k]*_dW[k] .+ b132*g2[k]*_dW[k]
        for l=1:m
          if l!=k
            H13[k] += b331*g1[l]*_dW[l] .+ b332*g2[l][l]*_dW[l]
          end
        end
      else
        H13 += b131*g1[:,k]*_dW[k] .+ b132*g2[k][:,k]*_dW[k]
        for l=1:m
          if l!=k
            H13[k] += b331*g1[:,l]*_dW[l] .+ b332*g2[l][:,l]*_dW[l]
          end
        end
      end
    end
  end

  if typeof(W.dW) <: Number
    g3 = integrator.g(H13,p,t+c13*dt)
  else
    g3 = [integrator.g(H13[k],p,t+c13*dt) for k=1:m]
  end
  # H_3^(k)
  if typeof(W.dW) <: Number
    H14 = uprev + a141*k1*dt + b141*g1*_dW + b142*g2*_dW + b143*g3*_dW
  else
    H14 = [zero(typeof(uprev)) for k=1:m]
    for k=1:m
      H14[k] += uprev .+ a141*k1*dt
      if is_diagonal_noise(integrator.sol.prob)
        H14[k] += b141*g1[k]*_dW[k] .+ b142*g2[k]*_dW[k] .+ b143*g3[k]*_dW[k]
        for l=1:m
          if l!=k
            H14[k] += b341*g1[l]*_dW[l] .+ b342*g2[l][l]*_dW[l]
          end
        end
      else
        H14 += b141*g1[:,k]*_dW[k] .+ b142*g2[k][:,k]*_dW[k] .+ b143*g3[k][:,k]*_dW[k]
        for l=1:m
          if l!=k
            H14[k] += b341*g1[:,l]*_dW[l] .+ b342*g2[l][:,l]*_dW[l]
          end
        end
      end
    end
  end

  if typeof(W.dW) <: Number
    g4 = integrator.g(H14,p,t+c14*dt)
  else
    g4 = [integrator.g(H14[k],p,t+c14*dt) for k=1:m]
  end

  # H_3^(0) (requires H_2^(k))
  k2 = integrator.f(H02,p,t+c02*dt)
  H03 = uprev + a032*k2*dt + a031*k1*dt

  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
    H03 += b031*g1*_dW
    if typeof(W.dW) <: Number
      H03 += b032*g2*_dW
    else
      for k=1:m
        H03 += b032*g2[k][:,k]*_dW[k]
      end
    end
  else
    H03 += b031*g1.*_dW
    for k=1:m
      H03 += b032*g2[k][k]*_dW[k]
    end
  end

  k3 = integrator.f(H03,p,t+c03*dt)

  # H_4^(0)
  # H04 = uprev; k4 = integrator.f(uprev,p,t+0*dt) = k1

  # H^_1^(k)
  # H21 = uprev
  # H^_2^(k) # H^_3^(k)

  if !(typeof(W.dW) <: Number)
    H22 = [uprev for k=1:m]
    H23 = [uprev for k=1:m]
    # add inbounds for speed if working properly
    for k=1:m
      for l=1:m
        if k == l
          continue
        end
        if is_diagonal_noise(integrator.sol.prob)
          @.. H22[k] += b221*g1[l]*Ihat2[k,l]/integrator.sqdt
          @.. H23[k] += b231*g1[l]*Ihat2[k,l]/integrator.sqdt
        else
          H22[k] += b221*g1[:,l]*Ihat2[k,l]/integrator.sqdt
          H23[k] += b231*g1[:,l]*Ihat2[k,l]/integrator.sqdt
        end
      end
    end
  end
  # H^_4^(k)
  # H24 = uprev

  # add stages together Eq. (5.1)
  u = uprev + (α1*k1 + α2*k2 + α3*k3 + α4*k1)*dt

  # add noise
  if typeof(W.dW) <: Number
    u += (g1*beta11+g2*beta12+g3*beta13+g4*beta14)*_dW # beta2 terms are zero by construction
  else
      if is_diagonal_noise(integrator.sol.prob)
        u += g1.*_dW*beta11
        for k=1:m
          u[k] += (g2[k][k]*beta12+g3[k][k]*beta13+g4[k][k]*beta14)*_dW[k]
          tmpg = integrator.g(H22[k],p,t)
          u[k] = u[k] + tmpg[k]*beta22*integrator.sqdt
          tmpg = integrator.g(H23[k],p,t)
          u[k] = u[k] + tmpg[k]*beta23*integrator.sqdt
        end
      else
        # non-diag noise
        for k=1:m
          g1k = @view g1[:,k]
          g2k = @view g2[k][:,k]
          g3k = @view g3[k][:,k]
          g4k = @view g4[k][:,k]
          @.. u = u + (g1k*beta11 + g2k*beta12 + g3k*beta13 + g4k*beta14)*_dW[k]
          tmpg = integrator.g(H22[k],p,t)
          @.. u = u + tmpg*beta22*integrator.sqdt
          tmpg = integrator.g(H23[k],p,t)
          @.. u = u + tmpg*beta23*integrator.sqdt
        end
      end
  end

  integrator.u = u
end

@muladd function perform_step!(integrator,cache::RSCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack _dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,g4,k1,k2,k3,H02,H03,H12,H13,H14,H22,H23,tmp1,tmpg = cache
  @unpack a021,a031,a032,a131,a141,b031,b032,b121,b131,b132,b141,b142,b143,b221,b231,b331,b332,b341,b342,α1,α2,α3,α4,c02,c03,c13,c14,beta11,beta12,beta13,beta14,beta22,beta23,NORMAL_ONESIX_QUANTILE = cache.tab

  m = length(W.dW)
  sq3dt = sqrt(3*dt)
  if typeof(W.dW) <: Union{SArray,Number}
    # tbd
    _dW = map(x -> calc_threepoint_random(sq3dt, NORMAL_ONESIX_QUANTILE, x), W.dW / sqrt(dt))
  else
    # define three-point distributed random variables
    @.. chi1 = W.dW / sqrt(dt)
    calc_threepoint_random!(_dW, sq3dt, NORMAL_ONESIX_QUANTILE, chi1)
    map!(x -> (x^2-dt)/2, chi1, _dW)

    # define two-point distributed random variables
    calc_twopoint_random!(_dZ,integrator.sqdt,W.dZ)

    for k = 1:m
      for l = 1:k-1
        Ihat2[k, l] = _dW[k]*_dZ[l]
        Ihat2[l, k] = -_dW[k]*_dZ[l]
      end
    end
  end

  # compute stage values
  integrator.f(k1,uprev,p,t)
  integrator.g(g1,uprev,p,t)

  # H_1^(0)
  # H01 = uprev
  # H_2^(0)
  @.. H02 = uprev + a021*k1*dt

  # H_1^(k)
  # H11 = uprev
  # H_2^(k), stage 2

  if is_diagonal_noise(integrator.sol.prob)
    for k=1:m
      @.. H12[k] = uprev + b121*g1[k]*_dW[k]
    end
  else
    for k=1:m
      g1k = @view g1[:,k]
      @.. H12[k] = uprev .+ b121*g1k*_dW[k]
    end
  end

  for k=1:m
    integrator.g(g2[k],H12[k],p,t)
  end

  # H_3^(k)
  for k=1:m
    @.. H13[k] = uprev+a131*k1*dt
    if is_diagonal_noise(integrator.sol.prob)
      @.. H13[k] = H13[k] + b131*g1[k]*_dW[k] + b132*g2[k]*_dW[k]
      for l=1:m
        if l!=k
          @.. H13[k] = H13[k] + b331*g1[l]*_dW[l] + b332*g2[l]*_dW[l]
        end
      end
    else
      g1k = @view g1[:,k]
      g2k = @view g2[k][:,k]
      @.. H13[k] = H13[k] + b131*g1k*_dW[k] + b132*g2k*_dW[k]
      for l=1:m
        if l!=k
          g1l = @view g1[:,l]
          g2l = @view g2[l][:,l]
          @.. H13[k] = H13[k] + b331*g1l*_dW[l] + b332*g2l*_dW[l]
        end
      end
    end
  end

  for k=1:m
    integrator.g(g3[k],H13[k],p,t+c13*dt)
  end

  # H_3^(k)
  for k=1:m
    @.. H14[k] = uprev + a141*k1*dt
    if is_diagonal_noise(integrator.sol.prob)
      @.. H14[k] = H14[k] + b141*g1[k]*_dW[k] + b142*g2[k]*_dW[k] + b143*g3[k]*_dW[k]
      for l=1:m
        if l!=k
          @.. H14[k] = H14[k] + b341*g1[l]*_dW[l] .+ b342*g2[l]*_dW[l]
        end
      end
    else
      g1k = @view g1[:,k]
      g2k = @view g2[k][:,k]
      g3k = @view g3[k][:,k]

      @.. H14[k] = H14[k] + b141*g1k*_dW[k] + b142*g2k*_dW[k] + b143*g3k*_dW[k]
      for l=1:m
        if l!=k
          g1l = @view g1[:,l]
          g2l = @view g2[l][:,l]
          @.. H14[k] = H14[k] + b341*g1l*_dW[l] + b342*g2l*_dW[l]
        end
      end
    end
  end

  for k=1:m
    integrator.g(g4[k],H14[k],p,t+c14*dt)
  end

  # H_3^(0) (requires H_2^(k))
  integrator.f(k2,H02,p,t+c02*dt)
  @.. H03 = uprev + a032*k2*dt + a031*k1*dt

  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
    mul!(tmp1,g1,_dW)
    @.. H03 = H03 + b031*tmp1
    for k=1:m
      g2k = @view g2[k][:,k]
      @.. H03 = H03 +  b032*g2k*_dW[k]
    end
  else
    @.. H03 = H03 + b031*g1*_dW
    for k=1:m
      @.. H03 = H03 + b032*g2[k][k]*_dW[k]
    end
  end

  integrator.f(k3,H03,p,t+c03*dt)

  # H_4^(0)
  # H04 = uprev; k4 = integrator.f(uprev,p,t+0*dt) = k1

  # H^_1^(k)
  # H21 = uprev
  # H^_2^(k) # H^_3^(k)

  if (!(typeof(W.dW) <: Number) || (m!=1))
    for k=1:m
     H22[k] = uprev
     H23[k] = uprev
    end
    for k=1:m
      for l=1:m
        if k == l
          continue
        end
        if is_diagonal_noise(integrator.sol.prob)
          @.. H22[k] = H22[k] + b221*g1[l]*Ihat2[k,l]/integrator.sqdt
          @.. H23[k] = H23[k] + b231*g1[l]*Ihat2[k,l]/integrator.sqdt
        else
          g1l = @view g1[:,l]
          @.. H22[k] = H22[k] + b221*g1l*Ihat2[k,l]/integrator.sqdt
          @.. H23[k] = H23[k] + b231*g1l*Ihat2[k,l]/integrator.sqdt
        end
      end
    end
  end
  # H^_4^(k)
  # H24 = uprev

  # add stages together Eq. (5.1)
  @.. u = uprev + (α1*k1 + α2*k2 + α3*k3 + α4*k1)*dt

  # add noise

  if is_diagonal_noise(integrator.sol.prob)
    @.. u = u + g1*_dW*beta11
    for k=1:m
      u[k] = u[k] + (g2[k][k]*beta12+g3[k][k]*beta13+g4[k][k]*beta14)*_dW[k]
      if m != 1
        integrator.g(tmpg,H22[k],p,t)
        u[k] = u[k] + tmpg[k]*beta22*integrator.sqdt
        integrator.g(tmpg,H23[k],p,t)
        u[k] = u[k] + tmpg[k]*beta23*integrator.sqdt
      end
    end
  else
    # non-diag noise
    for k=1:m
      g1k = @view g1[:,k]
      g2k = @view g2[k][:,k]
      g3k = @view g3[k][:,k]
      g4k = @view g4[k][:,k]
      tmpgk = @view tmpg[:,k]
      @.. u = u + (g1k*beta11 + g2k*beta12 + g3k*beta13 + g4k*beta14)*_dW[k]
      integrator.g(tmpg,H22[k],p,t)
      @.. u = u + tmpgk*beta22*integrator.sqdt
      integrator.g(tmpg,H23[k],p,t)
      @.. u = u + tmpgk*beta23*integrator.sqdt
    end
  end
end

# PL1WM
@muladd function perform_step!(integrator,cache::PL1WMConstantCache,f=integrator.f)
  @unpack NORMAL_ONESIX_QUANTILE = cache
  @unpack t,dt,uprev,u,W,p = integrator

  m = length(W.dW)
  sq3dt = sqrt(3*dt)
  # define three-point distributed random variables
  dW_scaled = W.dW / integrator.sqdt
  _dW = map(x -> calc_threepoint_random(sq3dt, NORMAL_ONESIX_QUANTILE, x), dW_scaled)
  chi1 = map(x -> (x^2-dt)/4, _dW)
  if !(typeof(W.dW) <: Number)
    # define two-point distributed random variables
    _dZ = map(x -> calc_twopoint_random(integrator.sqdt, x),  W.dZ)
    Ihat2 = zeros(eltype(W.dZ), m, m) # I^_(k,l)
    for k = 1:m
      Ihat2[k, k] = -dt
      for l = 1:k-1
        Ihat2[k, l] = _dZ[Int(1+1//2*(k-3)*k+l)]
      	Ihat2[l, k] = -Ihat2[k, l]
      end
    end
  end
  # compute stage values
  k1 = integrator.f(uprev,p,t)
  g1 = integrator.g(uprev,p,t)

  # Y, Yp, Ym
  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
    Y = uprev + k1*dt + g1*_dW
    if typeof(W.dW) <: Number
      Yp = uprev + k1*dt + g1*integrator.sqdt
      Ym = uprev + k1*dt - g1*integrator.sqdt
    else
      Yp = Vector{typeof(uprev)}[uprev .+ k1*dt .+ g1[:,k]*integrator.sqdt for k=1:m]
      Ym = Vector{typeof(uprev)}[uprev .+ k1*dt .- g1[:,k]*integrator.sqdt for k=1:m]
    end
  else
    Y = uprev + k1*dt + g1.*_dW
    Yp = Vector{typeof(uprev)}(undef,m)
    Ym = Vector{typeof(uprev)}(undef,m)
    for k=1:m
      Yp[k] = uprev .+ k1*dt .+ g1[k]*integrator.sqdt
      Ym[k] = uprev .+ k1*dt .- g1[k]*integrator.sqdt
    end
  end

  k2 = integrator.f(Y,p,t)

  # add stages together (Ch 15.1 Eq.(1.3))
  u = uprev + 1//2*(k1+k2)*dt

  # add noise
  if typeof(W.dW) <: Number
    g2p = integrator.g(Yp,p,t)
    g2m = integrator.g(Ym,p,t)
    u += 1//4*(g2p+g2m+2*g1)*_dW + (g2p-g2m)*chi1/integrator.sqdt #(1.1)
  else
    if is_diagonal_noise(integrator.sol.prob)
      for k=1:m
        tmpg1 = integrator.g(Yp[k],p,t)
        tmpg2 = integrator.g(Ym[k],p,t)
        @.. u += 1//4*(tmpg1[k]+tmpg2[k]+2*g1[k])*_dW[k]
        @.. u += (tmpg1[k]-tmpg2[k])*chi1[k]/integrator.sqdt
        for l=1:m
          if l!=k
            Ulp = @.. uprev + g1[l]*integrator.sqdt
            Ulm = @.. uprev - g1[l]*integrator.sqdt

            tmpg1 = integrator.g(Ulp,p,t)
            tmpg2 = integrator.g(Ulm,p,t)
            @.. u += 1//4*(tmpg1[k]+tmpg2[k]-2*g1[k])*_dW[k]
            @.. u += 1//4*(tmpg1[k]-tmpg2[k])*(_dW[k]*_dW[l] + Ihat2[l,k])/integrator.sqdt
          end
        end
      end
      else
        # non-diag noise
        for k=1:m
          tmpg1 = integrator.g(Yp[k],p,t)
          tmpg2 = integrator.g(Ym[k],p,t)

          u += 1//4*(tmpg1[:,k]+tmpg2[:,k]+2*g1[:,k])*_dW[k]
          u += (tmpg1[:,k]-tmpg2[:,k])*chi1[k]/integrator.sqdt
          for l=1:m
            if l!=k
              Ulp = uprev + g1[:,l]*integrator.sqdt
              Ulm = uprev - g1[:,l]*integrator.sqdt

              tmpg1 = integrator.g(Ulp,p,t)
              tmpg2 = integrator.g(Ulm,p,t)
              u += 1//4*(tmpg1[:,k]+tmpg2[:,k]-2*g1[:,k])*_dW[k]
              u += 1//4*(tmpg1[:,k]-tmpg2[:,k])*(_dW[k]*_dW[l] + Ihat2[l,k])/integrator.sqdt
            end
          end
        end
      end
  end
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::PL1WMCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack _dW,_dZ,chi1,Ihat2,tab,g1,k1,k2,Y,Yp,Ym,tmp1,tmpg1,tmpg2,Ulp,Ulm = cache
  @unpack NORMAL_ONESIX_QUANTILE = cache.tab

  m = length(W.dW)
  sq3dt = sqrt(3*dt)
  # define three-point distributed random variables
  @.. chi1 = W.dW / integrator.sqdt
  calc_threepoint_random!(_dW, sq3dt, NORMAL_ONESIX_QUANTILE, chi1)
  map!(x -> (x^2-dt)/4, chi1, _dW)
  if !(typeof(W.dW) <: Number) || m > 1
    # define two-point distributed random variables
    calc_twopoint_random!(_dZ,integrator.sqdt,W.dZ)
    for k = 1:m
      Ihat2[k, k] = -dt
      for l = 1:k-1
        Ihat2[k, l] = _dZ[Int(1+1//2*(k-3)*k+l)]
      	Ihat2[l, k] = -Ihat2[k, l]
      end
    end
  end
  # compute stage values
  integrator.f(k1,uprev,p,t)
  integrator.g(g1,uprev,p,t)

  # Y, Yp, Ym
  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
    mul!(tmp1,g1,_dW)
    @.. Y = uprev + k1*dt + tmp1
    if typeof(W.dW) <: Number
      @.. Yp = uprev + k1*dt + g1*integrator.sqdt
      @.. Ym = uprev + k1*dt - g1*integrator.sqdt
    else
      for k=1:m
        g1k = @view g1[:,k]
        @.. Yp[k] = uprev + k1*dt + g1k*integrator.sqdt
        @.. Ym[k] = uprev + k1*dt - g1k*integrator.sqdt
      end
    end
  else
    @.. Y = uprev + k1*dt + g1*_dW
    for k=1:m
      @.. Yp[k] = uprev + k1*dt + g1[k]*integrator.sqdt
      @.. Ym[k] = uprev + k1*dt - g1[k]*integrator.sqdt
    end
  end

  integrator.f(k2,Y,p,t)

  # add stages together (Ch 15.1 Eq.(1.3))
  @.. u = uprev + 1//2*(k1+k2)*dt

  # add noise
  if typeof(W.dW) <: Number
    integrator.g(tmpg1,Yp,p,t)
    integrator.g(tmpg2,Ym,p,t)
    @.. u = u + 1//4*(tmpg1+tmpg2+2*g1)*_dW + 1//4*(tmpg1-tmpg2)*chi1[k]/integrator.sqdt #(1.1)
  else
    if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
      # non-diag noise
      for k=1:m
        integrator.g(tmpg1,Yp[k],p,t)
        integrator.g(tmpg2,Ym[k],p,t)
        tmpg1k = @view tmpg1[:,k]
        tmpg2k = @view tmpg2[:,k]
        g1k = @view g1[:,k]
        @.. u = u + 1//4*(tmpg1k+tmpg2k+2*g1k)*_dW[k]
        @.. u = u + (tmpg1k-tmpg2k)*chi1[k]/integrator.sqdt
        for l=1:m
          if l !=k
            g1l = @view g1[:,l]
            @.. Ulp = uprev + g1l*integrator.sqdt
            @.. Ulm = uprev - g1l*integrator.sqdt

            integrator.g(tmpg1,Ulp,p,t)
            integrator.g(tmpg2,Ulm,p,t)
            tmpg1k = @view tmpg1[:,k]
            tmpg2k = @view tmpg2[:,k]
            u += 1//4*(tmpg1k+tmpg2k-2*g1k)*_dW[k]
            u += 1//4*(tmpg1k-tmpg2k)*(_dW[k]*_dW[l] + Ihat2[l,k])/integrator.sqdt
          end
        end
      end
    else
      for k=1:m
        integrator.g(tmpg1,Yp[k],p,t)
        integrator.g(tmpg2,Ym[k],p,t)
        @.. u = u + 1//4*(tmpg1[k]+tmpg2[k]+2*g1[k])*_dW[k]
        @.. u = u + (tmpg1[k]-tmpg2[k])*chi1[k]/integrator.sqdt
        for l=1:m
          if l !=k
            @.. Ulp = uprev + g1[l]*integrator.sqdt
            @.. Ulm = uprev - g1[l]*integrator.sqdt

            integrator.g(tmpg1,Ulp,p,t)
            integrator.g(tmpg2,Ulm,p,t)
            @.. u = u + 1//4*(tmpg1[k]+tmpg2[k]-2*g1[k])*_dW[k]
            @.. u = u + 1//4*(tmpg1[k]-tmpg2[k])*(_dW[k]*_dW[l] + Ihat2[l,k])/integrator.sqdt
          end
        end
      end
    end
  end
end

# PL1WM
@muladd function perform_step!(integrator,cache::PL1WMAConstantCache,f=integrator.f)
  @unpack NORMAL_ONESIX_QUANTILE = cache
  @unpack t,dt,uprev,u,W,p = integrator

  sq3dt = sqrt(3*dt)
  # define three-point distributed random variables
  dW_scaled = W.dW / integrator.sqdt
  _dW = map(x -> calc_threepoint_random(sq3dt, NORMAL_ONESIX_QUANTILE, x), dW_scaled)

  # compute stage values
  k1 = integrator.f(uprev,p,t)
  g1 = integrator.g(uprev,p,t)

  # Y
  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
    tmp1 = g1*_dW
  else
    tmp1 = g1.*_dW
  end
  Y = uprev + k1*dt + tmp1

  k2 = integrator.f(Y,p,t)

  # add stages together (Ch 15.1 Eq.(1.4))
  u = uprev + 1//2*(k1+k2)*dt + tmp1

  integrator.u = u
end

@muladd function perform_step!(integrator,cache::PL1WMACache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack _dW,chi1,tab,g1,k1,k2,Y,tmp1 = cache
  @unpack NORMAL_ONESIX_QUANTILE = cache.tab

  sq3dt = sqrt(3*dt)
  # define three-point distributed random variables
  @.. chi1 = W.dW / integrator.sqdt
  calc_threepoint_random!(_dW, sq3dt, NORMAL_ONESIX_QUANTILE, chi1)

  # compute stage values
  integrator.f(k1,uprev,p,t)
  integrator.g(g1,uprev,p,t)

  # Y, Yp, Ym
  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
    mul!(tmp1,g1,_dW)
    @.. Y = uprev + k1*dt + tmp1

  else
    @.. tmp1 = g1*_dW
    @.. Y = uprev + k1*dt + tmp1
  end

  integrator.f(k2,Y,p,t)

  # add stages together (Ch 15.1 Eq.(1.4))
  @.. u = uprev + 1//2*(k1+k2)*dt + tmp1

end

# NON
@muladd function perform_step!(integrator,cache::NONConstantCache,f=integrator.f)
  @unpack c01,c02,c03,c04,cj1,cj2,cj3,cj4,cjl2,cjl3,clj2,clj3,a0021,a0032,a0043,aj021,aj041,a0j21,a0j31,a0j32,a0j41,ajj21,ajj31,ajj32,ajj41,ajj42,ajj43,ajl31,ajl32,ajl41,ajl42,ajljj31,aljjl21,aljjl31,NORMAL_ONESIX_QUANTILE = cache
  @unpack t,dt,uprev,u,W,p = integrator

  m = length(W.dW)
  # define three-point distributed random variables
  sq3dt = sqrt(3*dt)
  dW_scaled = W.dW / integrator.sqdt
  _dW = map(x -> calc_threepoint_random(sq3dt, NORMAL_ONESIX_QUANTILE, x), dW_scaled)

  if !(typeof(W.dW) <: Number)
    # define two-point distributed random variables
    _dZ = map(x -> calc_twopoint_random(integrator.sqdt, x),  W.dZ)
    Ihat2 = zeros(eltype(W.dZ), m, m) # I^_(k,l)
    for k = 1:m
      Ihat2[k, k] = _dW[k]
      for l = 1:k-1
        Ihat2[k, l] = _dW[l]
        Ihat2[l, k] = _dZ[k]
      end
    end
  end

  # compute stage values
  # stage 1
  ktmp = integrator.f(uprev,p,t)
  gtmp = integrator.g(uprev,p,t)

  Y100 = ktmp*dt
  Y1jajb = [zero(u) for ja=1:m, jb=1:m]
  Y2jajb = [zero(u) for ja=1:m, jb=1:m]
  Y3jajb = [zero(u) for ja=1:m, jb=1:m]
  Y4jajb = [zero(u) for ja=1:m, jb=1:m]
  if typeof(W.dW) <: Number
    Y1jajb = gtmp*_dW
  else
    for ja=1:m
      for jb=1:m
        if is_diagonal_noise(integrator.sol.prob)
          tmpu = zero(integrator.u)
          tmpu[jb] = gtmp[jb]
          @.. Y1jajb[ja,jb] = Ihat2[ja,jb]*tmpu
          if ja!=jb
            @.. Y4jajb[ja,jb] = Ihat2[ja,jb]*tmpu
          end
        else
          gtmpjb = @view gtmp[:,jb]
          @.. Y1jajb[ja,jb] = Ihat2[ja,jb]*gtmpjb
          if ja!=jb
            @.. Y4jajb[ja,jb] = Ihat2[ja,jb]*gtmpjb
          end
        end
      end
    end
  end

  # stage 2
  Y200 = uprev + a0021*Y100

  if typeof(W.dW) <: Number
    Y200 += a0j21*Y1jajb
    Y2jajb = integrator.g(uprev + aj021*Y100 + ajj21*Y1jajb[1],p,t)*_dW
  else
    if !is_diagonal_noise(integrator.sol.prob)
      for ja=1:m
        Y200 += a0j21*Y1jajb[ja,ja]
        for jb=1:m
          if jb!=ja
            if jb<ja
              tmpu = uprev + aljjl21*Y1jajb[jb,ja]
            else
              tmpu = copy(uprev)
            end
          else
            tmpu = uprev + aj021*Y100 + ajj21*Y1jajb[ja,ja]
          end
          gtmp = integrator.g(tmpu,p,t)
          tmpjb = @view gtmp[:,jb]
          @.. Y2jajb[ja,jb] = Ihat2[ja,jb]*tmpjb
        end
      end
    else
      for ja=1:m
        Y200 += a0j21*Y1jajb[ja,ja]
        for jb=1:m
          if jb!=ja
            if jb<ja
              tmpu = uprev + aljjl21*Y1jajb[jb,ja]
            else
              tmpu = copy(uprev)
            end
          else
            tmpu = uprev + aj021*Y100 + ajj21*Y1jajb[ja,ja]
          end
          gtmp = integrator.g(tmpu,p,t)
          tmpu = zero(integrator.u)
          tmpu[jb] = gtmp[jb]
          @.. Y2jajb[ja,jb] = Ihat2[ja,jb]*tmpu
        end
      end
    end
  end
  ktmp = integrator.f(Y200,p,t)
  Y200 = ktmp*dt

  # stage 3
  Y300 = uprev + a0032*Y200
  if typeof(W.dW) <: Number
    Y300 += a0j31*Y1jajb + a0j32*Y2jajb
    Y3jajb = integrator.g(uprev + ajj31*Y1jajb[1] + ajj32*Y2jajb[1],p,t)*_dW
  else
    if !is_diagonal_noise(integrator.sol.prob)
      for ja=1:m
        Y300 += a0j31*Y1jajb[ja,ja] + a0j32*Y2jajb[ja,ja]
        for jb=1:m
          if jb!=ja
            if jb<ja
              tmpu = uprev + aljjl31*Y1jajb[jb,ja]
            else
              tmpu = uprev + ajljj31*Y1jajb[ja,ja]
            end
          else
            tmpu = uprev + ajj31*Y1jajb[ja,ja] + ajj32*Y2jajb[ja,ja]
            for l=1:m
              if l!=ja
                tmpu += ajl31*Y1jajb[l,l] + ajl32*Y2jajb[l,l]
              end
            end
          end
          gtmp = integrator.g(tmpu,p,t)
          tmpjb = @view gtmp[:,jb]
          @.. Y3jajb[ja,jb] = Ihat2[ja,jb]*tmpjb
        end
      end
    else
      for ja=1:m
        Y300 += a0j31*Y1jajb[ja,ja] + a0j32*Y2jajb[ja,ja]
        for jb=1:m
          if jb!=ja
            if jb<ja
              tmpu = uprev + aljjl31*Y1jajb[jb,ja]
            else
              tmpu = uprev + ajljj31*Y1jajb[ja,ja]
            end
          else
            tmpu = uprev + ajj31*Y1jajb[ja,ja] + ajj32*Y2jajb[ja,ja]
            for l=1:m
              if l!=ja
                tmpu += ajl31*Y1jajb[l,l] + ajl32*Y2jajb[l,l]
              end
            end
          end
          gtmp = integrator.g(tmpu,p,t)
          tmpu = zero(integrator.u)
          tmpu[jb] = gtmp[jb]
          @.. Y3jajb[ja,jb] = Ihat2[ja,jb]*tmpu
        end
      end
    end
  end
  ktmp = integrator.f(Y300,p,t)
  Y300 = ktmp*dt

  # stage 4
  Y400 = uprev + a0043*Y300
  if typeof(W.dW) <: Number
    Y400 += a0j41*Y1jajb
    Y4jajb = integrator.g(uprev + ajj41*Y1jajb[1] + ajj42*Y2jajb[1] + ajj43*Y3jajb[1] + aj041*Y100,p,t)*_dW
  else
    if !is_diagonal_noise(integrator.sol.prob)
      for ja=1:m
        Y400 += a0j41*Y1jajb[ja,ja]
        tmpu = uprev + ajj41*Y1jajb[ja,ja] + ajj42*Y2jajb[ja,ja] + ajj43*Y3jajb[ja,ja] + aj041*Y100
        for l=1:m
          if l!=ja
            tmpu += ajl41*Y1jajb[l,l] + ajl42*Y2jajb[l,l]
          end
        end
        gtmp = integrator.g(tmpu,p,t)
        tmpjb = @view gtmp[:,ja]
        @.. Y4jajb[ja,ja] = Ihat2[ja,ja]*tmpjb
      end
    else
      for ja=1:m
        Y400 += a0j41*Y1jajb[ja,ja]

        tmpu = uprev + ajj41*Y1jajb[ja,ja] + ajj42*Y2jajb[ja,ja] + ajj43*Y3jajb[ja,ja] + aj041*Y100
        for l=1:m
          if l!=ja
            tmpu += ajl41*Y1jajb[l,l] + ajl42*Y2jajb[l,l]
          end
        end
        gtmp = integrator.g(tmpu,p,t)
        tmpu = zero(integrator.u)
        tmpu[ja] = gtmp[ja]
        @.. Y4jajb[ja,ja] = Ihat2[ja,ja]*tmpu
      end
    end
  end

  ktmp = integrator.f(Y400,p,t)
  Y400 = ktmp*dt

  # add stages together
  u = uprev + c01*Y100+c02*Y200+c03*Y300+c04*Y400
  if typeof(W.dW) <: Number
    # to be checked
    u += cj1*Y1jajb + cj2*Y2jajb + cj3*Y3jajb + cj4*Y4jajb
  else
    if is_diagonal_noise(integrator.sol.prob)
      for ja=1:m
        u += cj1*Y1jajb[ja,ja] + cj2*Y2jajb[ja,ja] + cj3*Y3jajb[ja,ja] + cj4*Y4jajb[ja,ja]
        for jb=1:m
          if jb!=ja
            if ja>jb
              u += clj2 * Y2jajb[ja,jb] + clj3 * Y3jajb[ja,jb]
            else
              u += cjl2 * Y2jajb[ja,jb] + cjl3 * Y3jajb[ja,jb]
            end
          end
        end
      end
    else
      for ja=1:m
        u += cj1*Y1jajb[ja,ja] + cj2*Y2jajb[ja,ja] + cj3*Y3jajb[ja,ja] + cj4*Y4jajb[ja,ja]
        for jb=1:m
          if jb!=ja
            if ja>jb
              u += clj2 * Y2jajb[ja,jb] + clj3 * Y3jajb[ja,jb]
            else
              u += cjl2 * Y2jajb[ja,jb] + cjl3 * Y3jajb[ja,jb]
            end
          end
        end
      end
    end
  end
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::NONCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack _dW,_dZ,chi1,Ihat2,tab,gtmp,ktmp,Y100,Y200,Y300,Y400,Y1jajb,Y2jajb,Y3jajb,Y4jajb,tmpu = cache
  @unpack c01,c02,c03,c04,cj1,cj2,cj3,cj4,cjl2,cjl3,clj2,clj3,a0021,a0032,a0043,aj021,aj041,a0j21,a0j31,a0j32,a0j41,ajj21,ajj31,ajj32,ajj41,ajj42,ajj43,ajl31,ajl32,ajl41,ajl42,ajljj31,aljjl21,aljjl31,NORMAL_ONESIX_QUANTILE = cache.tab

  m = length(W.dW)
  # define three-point distributed random variables
  @.. chi1 = W.dW / integrator.sqdt
  sq3dt = sqrt(3*dt)
  calc_threepoint_random!(_dW, sq3dt, NORMAL_ONESIX_QUANTILE, chi1)

  if !(typeof(W.dW) <: Number)
    # define two-point distributed random variables
    calc_twopoint_random!(_dZ,integrator.sqdt,W.dZ)
    for k = 1:m
      Ihat2[k, k] = _dW[k]
      for l = 1:k-1
       Ihat2[k, l] = _dW[l]
       Ihat2[l, k] = _dZ[k]
      end
    end
  end

  # compute stage values
  # stage 1
  integrator.f(ktmp,uprev,p,t)
  integrator.g(gtmp,uprev,p,t)

  @.. Y100 = ktmp*dt
  for ja=1:m
    for jb=1:m
      if is_diagonal_noise(integrator.sol.prob)
        fill!(tmpu,zero(eltype(integrator.u)))
        tmpu[jb] = gtmp[jb]
        @.. Y1jajb[ja,jb] = Ihat2[ja,jb]*tmpu
        if ja!=jb
          @.. Y4jajb[ja,jb] = Ihat2[ja,jb]*tmpu
        end
      else
        gtmpjb = @view gtmp[:,jb]
        @.. Y1jajb[ja,jb] = Ihat2[ja,jb]*gtmpjb
        if ja!=jb
          @.. Y4jajb[ja,jb] = Ihat2[ja,jb]*gtmpjb
        end
      end
    end
  end

  # stage 2
  @.. Y200 = uprev + a0021*Y100

  if !is_diagonal_noise(integrator.sol.prob)
    for ja=1:m
      @.. Y200 = Y200 + a0j21*Y1jajb[ja,ja]
      for jb=1:m
        if jb!=ja
          if jb<ja
            @.. tmpu = uprev + aljjl21*Y1jajb[jb,ja]
          else
            @.. tmpu = uprev
          end
        else
          @.. tmpu = uprev + aj021*Y100 + ajj21*Y1jajb[ja,ja]
        end
        integrator.g(gtmp,tmpu,p,t)
        tmpjb = @view gtmp[:,jb]
        @.. Y2jajb[ja,jb] = Ihat2[ja,jb]*tmpjb
        end
    end
  else
    for ja=1:m
      @.. Y200 = Y200 + a0j21*Y1jajb[ja,ja]
      for jb=1:m
        if jb!=ja
          if jb<ja
            @.. tmpu = uprev + aljjl21*Y1jajb[jb,ja]
          else
            @.. tmpu = uprev
          end
        else
          @.. tmpu = uprev + aj021*Y100 + ajj21*Y1jajb[ja,ja]
        end
        integrator.g(gtmp,tmpu,p,t)
        fill!(tmpu,zero(eltype(integrator.u)))
        tmpu[jb] = gtmp[jb]
        @.. Y2jajb[ja,jb] = Ihat2[ja,jb]*tmpu
      end
    end
  end
  integrator.f(ktmp,Y200,p,t)
  @.. Y200 = ktmp*dt

  # stage 3
  @.. Y300 = uprev + a0032*Y200

  if !is_diagonal_noise(integrator.sol.prob)
      for ja=1:m
        @.. Y300 = Y300 + a0j31*Y1jajb[ja,ja] + a0j32*Y2jajb[ja,ja]
        for jb=1:m
          if jb!=ja
            if jb<ja
              @.. tmpu = uprev + aljjl31*Y1jajb[jb,ja]
            else
              @.. tmpu = uprev + ajljj31*Y1jajb[ja,ja]
            end
          else
            @.. tmpu = uprev + ajj31*Y1jajb[ja,ja] + ajj32*Y2jajb[ja,ja]
            for l=1:m
              if l!=ja
                @.. tmpu = tmpu + ajl31*Y1jajb[l,l] + ajl32*Y2jajb[l,l]
              end
            end
          end
          integrator.g(gtmp,tmpu,p,t)
          tmpjb = @view gtmp[:,jb]
          @.. Y3jajb[ja,jb] = Ihat2[ja,jb]*tmpjb
        end
      end
  else
    for ja=1:m
      @.. Y300 = Y300 + a0j31*Y1jajb[ja,ja] + a0j32*Y2jajb[ja,ja]
      for jb=1:m
        if jb!=ja
          if jb<ja
            @.. tmpu = uprev + aljjl31*Y1jajb[jb,ja]
          else
            @.. tmpu = uprev + ajljj31*Y1jajb[ja,ja]
          end
        else
          @.. tmpu = uprev + ajj31*Y1jajb[ja,ja] + ajj32*Y2jajb[ja,ja]
          for l=1:m
            if l!=ja
              @.. tmpu = tmpu + ajl31*Y1jajb[l,l] + ajl32*Y2jajb[l,l]
            end
          end
        end
        integrator.g(gtmp,tmpu,p,t)
        fill!(tmpu,zero(eltype(integrator.u)))
        tmpu[jb] = gtmp[jb]
        @.. Y3jajb[ja,jb] = Ihat2[ja,jb]*tmpu
      end
    end
  end

  integrator.f(ktmp,Y300,p,t)
  @.. Y300 = ktmp*dt

  # stage 4
  @.. Y400 = uprev + a0043*Y300

  if !is_diagonal_noise(integrator.sol.prob)
    for ja=1:m
      @.. Y400 = Y400 + a0j41*Y1jajb[ja,ja]

      @.. tmpu = uprev + ajj41*Y1jajb[ja,ja] + ajj42*Y2jajb[ja,ja] + ajj43*Y3jajb[ja,ja] + aj041*Y100
      for l=1:m
        if l!=ja
          @.. tmpu = tmpu + ajl41*Y1jajb[l,l] + ajl42*Y2jajb[l,l]
        end
      end
      integrator.g(gtmp,tmpu,p,t)
      tmpjb = @view gtmp[:,ja]
      @.. Y4jajb[ja,ja] = Ihat2[ja,ja]*tmpjb
    end
  else
    for ja=1:m
      @.. Y400 = Y400 + a0j41*Y1jajb[ja,ja]

      @.. tmpu = uprev + ajj41*Y1jajb[ja,ja] + ajj42*Y2jajb[ja,ja] + ajj43*Y3jajb[ja,ja] + aj041*Y100
      for l=1:m
        if l!=ja
          @.. tmpu = tmpu + ajl41*Y1jajb[l,l] + ajl42*Y2jajb[l,l]
        end
      end
      integrator.g(gtmp,tmpu,p,t)
      fill!(tmpu,zero(eltype(integrator.u)))
      tmpu[ja] = gtmp[ja]
      @.. Y4jajb[ja,ja] = Ihat2[ja,ja]*tmpu
    end
  end

  integrator.f(ktmp,Y400,p,t)
  @.. Y400 = ktmp*dt

  # add stages together
  @.. u = uprev + c01*Y100+c02*Y200+c03*Y300+c04*Y400
  if typeof(W.dW) <: Number
    @.. u = u + cj1*Y1jj + cj2*Y2jj + cj3*Y3jj + cj4*Y4jj
  else
    if is_diagonal_noise(integrator.sol.prob)
      for ja=1:m
        @.. u = u + cj1*Y1jajb[ja,ja] + cj2*Y2jajb[ja,ja] + cj3*Y3jajb[ja,ja] + cj4*Y4jajb[ja,ja]
        for jb=1:m
          if jb!=ja
            if ja>jb
              @.. u = u + clj2 * Y2jajb[ja,jb] + clj3 * Y3jajb[ja,jb]
            else
              @.. u = u + cjl2 * Y2jajb[ja,jb] + cjl3 * Y3jajb[ja,jb]
            end
          end
        end
      end
    else
      for ja=1:m
        @.. u = u + cj1*Y1jajb[ja,ja] + cj2*Y2jajb[ja,ja] + cj3*Y3jajb[ja,ja] + cj4*Y4jajb[ja,ja]
        for jb=1:m
          if jb!=ja
            if ja>jb
              @.. u = u + clj2 * Y2jajb[ja,jb] + clj3 * Y3jajb[ja,jb]
            else
              @.. u = u + cjl2 * Y2jajb[ja,jb] + cjl3 * Y3jajb[ja,jb]
            end
          end
        end
      end
    end
  end
end

# COM
@muladd function perform_step!(integrator,cache::COMConstantCache,f=integrator.f)
  @unpack c01,c02,c03,c04,cj1,cj2,cj3,cj4,a0021,a0032,a0043,aj021,aj041,a0j21,a0j31,a0j32,a0j41,ajj21,ajj31,ajj32,ajj41,ajj42,ajj43,ajl31,ajl32,ajl41,ajl42,NORMAL_ONESIX_QUANTILE = cache
  @unpack t,dt,uprev,u,W,p = integrator

  m = length(W.dW)
  _dW = W.dW

  # define three-point distributed random variables
  #sq3dt = sqrt(3*dt)
  #dW_scaled = W.dW / integrator.sqdt
  #_dW = map(x -> calc_threepoint_random(sq3dt, NORMAL_ONESIX_QUANTILE, x), dW_scaled)

  # compute stage values
  # stage 1
  ktmp = integrator.f(uprev,p,t)
  gtmp = integrator.g(uprev,p,t)

  Y10 = ktmp*dt

  if typeof(W.dW) <: Number
    Y1j = gtmp*_dW
  else
    if is_diagonal_noise(integrator.sol.prob)
      Y1j = gtmp.*_dW
    else
      Y1j = gtmp.*_dW'
    end
  end

  # stage 2
  Y20 = uprev + a0021*Y10

  if typeof(W.dW) <: Number
    Y20 += a0j21*Y1j
    Y2j = integrator.g(uprev + aj021*Y10 + ajj21*Y1j,p,t)*_dW
  else
    if !is_diagonal_noise(integrator.sol.prob)
      Y2j = similar(integrator.sol.prob.noise_rate_prototype)
      for j = 1:m
        Y20 += a0j21*Y1j[:,j]
        tmpu = uprev + aj021*Y10 + ajj21*Y1j[:,j]
        Y2j[:,j] = integrator.g(tmpu,p,t)[:,j]*_dW[j]
      end
    else
      Y20 += a0j21*Y1j
      Y2j = similar(uprev)
      for j = 1:m
        tmpu = uprev + aj021*Y10
        tmpu[j] += ajj21*Y1j[j]
        Y2j[j] = integrator.g(tmpu,p,t)[j]*_dW[j]
      end
    end
  end
  ktmp = integrator.f(Y20,p,t)
  Y20 = ktmp*dt

  # stage 3
  Y30 = uprev + a0032*Y20
  if typeof(W.dW) <: Number
    Y30 += a0j31*Y1j + a0j32*Y2j
    Y3j = integrator.g(uprev + ajj31*Y1j + ajj32*Y2j,p,t)*_dW
  else
    if !is_diagonal_noise(integrator.sol.prob)
      Y3j = similar(integrator.sol.prob.noise_rate_prototype)
      for j = 1:m
        tmpu = uprev + ajj31*Y1j[:,j] + ajj32*Y2j[:,j]
        Y30 += a0j31*Y1j[:,j] + a0j32*Y2j[:,j]
        for l = 1:m
          if l!=j
            tmpu += (ajl31*Y1j[:,l] + ajl32*Y2j[:,l])
          end
        end
        Y3j[:,j] = integrator.g(tmpu,p,t)[:,j]*_dW[j]
      end
    else
      Y30 += a0j31*Y1j + a0j32*Y2j
      Y3j = similar(uprev)
      tmpu2 = uprev + ajl31*Y1j + ajl32*Y2j
      for j = 1:m
        tmpu = copy(tmpu2)
        tmpu[j] += ajj31*Y1j[j] + ajj32*Y2j[j] - (ajl31*Y1j[j] + ajl32*Y2j[j])
        Y3j[j] = integrator.g(tmpu,p,t)[j]*_dW[j]
      end
    end
  end
  ktmp = integrator.f(Y30,p,t)
  Y30 = ktmp*dt

  # stage 4
  Y40 = uprev + a0043*Y30
  if typeof(W.dW) <: Number
    Y40 += a0j41*Y1j
    Y4j = integrator.g(uprev + ajj41*Y1j + ajj42*Y2j + ajj43*Y3j + aj041*Y10,p,t)*_dW
  else
    if !is_diagonal_noise(integrator.sol.prob)
      Y4j = similar(integrator.sol.prob.noise_rate_prototype)
      for j = 1:m
        tmpu = uprev + ajj41*Y1j[:,j] + ajj42*Y2j[:,j] + ajj43*Y3j[:,j] + aj041*Y10
        Y40 += a0j41*Y1j[:,j]
        for l = 1:m
          if l!=j
            tmpu += (ajl41*Y1j[:,l] + ajl42*Y2j[:,l])
          end
        end
        Y4j[:,j] = integrator.g(tmpu,p,t)[:,j]*_dW[j]
    end
    else
      Y40 += a0j41*Y1j
      Y4j = similar(uprev)
      tmpu2 = uprev + ajl41*Y1j + ajl42*Y2j + aj041*Y10
      for j = 1:m
        tmpu = copy(tmpu2)
        tmpu[j] += ajj41*Y1j[j] + ajj42*Y2j[j] + ajj43*Y3j[j] - (ajl41*Y1j[j] + ajl42*Y2j[j])
        Y4j[j] = integrator.g(tmpu,p,t)[j]*_dW[j]
      end
    end
  end
  ktmp = integrator.f(Y40,p,t)
  Y40 = ktmp*dt

  # add stages together
  u = uprev + c01*Y10+c02*Y20+c03*Y30+c04*Y40

  if typeof(W.dW) <: Number
    u += cj1*Y1j + cj2*Y2j + cj3*Y3j + cj4*Y4j
  else
    if is_diagonal_noise(integrator.sol.prob)
      u += @. cj1*Y1j + cj2*Y2j + cj3*Y3j + cj4*Y4j
    else
      for j = 1:m
        u += @. cj1*Y1j[:,j] + cj2*Y2j[:,j] + cj3*Y3j[:,j] + cj4*Y4j[:,j]
      end
    end
  end
  integrator.u = u
end

# COM
@muladd function perform_step!(integrator,cache::COMCache,f=integrator.f)
  @unpack _dW,tab,gtmp,ktmp,Y10,Y20,Y30,Y40,Y1j,Y2j,Y3j,Y4j,tmpu,tmpu2 = cache
  @unpack c01,c02,c03,c04,cj1,cj2,cj3,cj4,a0021,a0032,a0043,aj021,aj041,a0j21,a0j31,a0j32,a0j41,ajj21,ajj31,ajj32,ajj41,ajj42,ajj43,ajl31,ajl32,ajl41,ajl42,NORMAL_ONESIX_QUANTILE = cache.tab
  @unpack t,dt,uprev,u,W,p = integrator

  m = length(W.dW)
  @.. _dW =  W.dW

  # define three-point distributed random variables
  #sq3dt = sqrt(3*dt)
  #dW_scaled = W.dW / integrator.sqdt
  #_dW = map(x -> calc_threepoint_random(sq3dt, NORMAL_ONESIX_QUANTILE, x), dW_scaled)

  # compute stage values
  # stage 1
  integrator.f(ktmp,uprev,p,t)
  integrator.g(gtmp,uprev,p,t)

  @.. Y10 = ktmp*dt

  if is_diagonal_noise(integrator.sol.prob)
    @.. Y1j = gtmp*_dW
  else
    for j = 1:m
      Y1jj = @view Y1j[:,j]
      gtmpj = @view gtmp[:,j]
      @.. Y1jj = gtmpj*_dW[j]
    end
  end

  # stage 2
  @.. Y20 = uprev + a0021*Y10

  if !is_diagonal_noise(integrator.sol.prob)
    fill!(Y2j, zero(eltype(integrator.sol.prob.noise_rate_prototype)))
    for j = 1:m
      Y1jj = @view Y1j[:,j]
      Y2jj = @view Y2j[:,j]
      @.. Y20 = Y20 + a0j21*Y1jj
      @.. tmpu = uprev + aj021*Y10 + ajj21*Y1jj
      integrator.g(gtmp,tmpu,p,t)
      gtmpj = @view gtmp[:,j]
      @.. Y2jj = gtmpj*_dW[j]
    end
  else
    fill!(Y2j, zero(eltype(integrator.u)))
    @.. Y20 = Y20 + a0j21*Y1j
    @.. tmpu = uprev + aj021*Y10
    for j = 1:m
      fill!(tmpu2,zero(eltype(integrator.u)))
      tmpu2[j] = ajj21*Y1j[j]
      @.. tmpu2 = tmpu + tmpu2
      integrator.g(gtmp,tmpu2,p,t)
      Y2j[j]= gtmp[j]*_dW[j]
    end
  end
  integrator.f(ktmp,Y20,p,t)
  @.. Y20 = ktmp*dt

  # stage 3
  @.. Y30 = uprev + a0032*Y20

  if !is_diagonal_noise(integrator.sol.prob)
    fill!(Y3j, zero(eltype(integrator.sol.prob.noise_rate_prototype)))
    for j = 1:m
      Y1jj = @view Y1j[:,j]
      Y2jj = @view Y2j[:,j]
      Y3jj = @view Y3j[:,j]
      @.. tmpu = uprev + ajj31*Y1jj + ajj32*Y2jj
      @.. Y30 = Y30 + a0j31*Y1jj + a0j32*Y2jj
      for l = 1:m
        if l!=j
          Y1jl = @view Y1j[:,l]
          Y2jl = @view Y2j[:,l]
          @.. tmpu = tmpu + (ajl31*Y1jl + ajl32*Y2jl)
        end
      end
      integrator.g(gtmp,tmpu,p,t)
      gtmpj = @view gtmp[:,j]
      @.. Y3jj = gtmpj*_dW[j]
    end
  else
    @.. Y30 = Y30 + a0j31*Y1j + a0j32*Y2j
    fill!(Y3j, zero(eltype(integrator.u)))
    @.. tmpu2 = uprev + ajl31*Y1j + ajl32*Y2j
    for j = 1:m
      @.. tmpu = tmpu2
      tmpu[j] += ajj31*Y1j[j] + ajj32*Y2j[j] - (ajl31*Y1j[j] + ajl32*Y2j[j])
      integrator.g(gtmp,tmpu,p,t)
      Y3j[j] = gtmp[j]*_dW[j]
    end
  end
  integrator.f(ktmp,Y30,p,t)
  @.. Y30 = ktmp*dt

  # stage 4
  @.. Y40 = uprev + a0043*Y30

  if !is_diagonal_noise(integrator.sol.prob)
    fill!(Y4j, zero(eltype(integrator.sol.prob.noise_rate_prototype)))
    for j = 1:m
      Y1jj = @view Y1j[:,j]
      Y2jj = @view Y2j[:,j]
      Y3jj = @view Y3j[:,j]
      Y4jj = @view Y4j[:,j]
      @.. tmpu = uprev + ajj41*Y1jj + ajj42*Y2jj + ajj43*Y3jj + aj041*Y10
      @.. Y40 = Y40 + a0j41*Y1jj
      for l = 1:m
        if l!=j
          Y1jl = @view Y1j[:,l]
          Y2jl = @view Y2j[:,l]
          @.. tmpu = tmpu + (ajl41*Y1jl + ajl42*Y2jl)
        end
      end
      integrator.g(gtmp,tmpu,p,t)
      gtmpj = @view gtmp[:,j]
      @.. Y4jj = gtmpj*_dW[j]
  end
  else
    @.. Y40 = Y40 + a0j41*Y1j
    fill!(Y4j, zero(eltype(integrator.u)))
    @.. tmpu2 = uprev + ajl41*Y1j + ajl42*Y2j + aj041*Y10
    for j = 1:m
      @.. tmpu = tmpu2
      tmpu[j] += ajj41*Y1j[j] + ajj42*Y2j[j] + ajj43*Y3j[j] - (ajl41*Y1j[j] + ajl42*Y2j[j])
      integrator.g(gtmp,tmpu,p,t)
      Y4j[j] = gtmp[j]*_dW[j]
    end
  end
  integrator.f(ktmp,Y40,p,t)
  @.. Y40 = ktmp*dt

  # add stages together
  @.. u = uprev + c01*Y10+c02*Y20+c03*Y30+c04*Y40

  if is_diagonal_noise(integrator.sol.prob)
    @. u = u + cj1*Y1j + cj2*Y2j + cj3*Y3j + cj4*Y4j
  else
    for j = 1:m
      Y1jj = @view Y1j[:,j]
      Y2jj = @view Y2j[:,j]
      Y3jj = @view Y3j[:,j]
      Y4jj = @view Y4j[:,j]
      @.. u = u + cj1*Y1jj + cj2*Y2jj + cj3*Y3jj + cj4*Y4jj
    end
  end
end

# NON2
@muladd function perform_step!(integrator,cache::NON2ConstantCache,f=integrator.f)
  @unpack c01,c02,c03,c04,cj1,cj2,cj3,cj4,a0021,a0032,a0043,aj021,aj041,a0j21,a0j31,a0j32,a0j41,ajj21,ajj31,ajj32,ajj41,ajj42,ajj43,ajl31,ajl32,ajl41,ajl42,γ,NORMAL_ONESIX_QUANTILE = cache
  @unpack t,dt,uprev,u,W,p = integrator

  m = length(W.dW)
  #_dW = W.dW

  # define three-point distributed random variables
  sq3dt = sqrt(3*dt)
  chi1 = W.dW / integrator.sqdt
  _dW = map(x -> calc_threepoint_random(sq3dt, NORMAL_ONESIX_QUANTILE, x), chi1)

  if !(typeof(W.dW) <: Number)
    # define two-point distributed random variables
    _dZ = map(x -> calc_twopoint_random(integrator.sqdt, x),  W.dZ)
    Ihat2 = zeros(eltype(W.dZ), m, m) # I^_(k,l)
    for j = 1:m
      for l = 1:j-1
        Ihat2[j, l] = -_dZ[j]*_dW[l]/integrator.sqdt
        Ihat2[l, j] = _dW[l]*_dZ[j]/integrator.sqdt
      end
    end
  end
  # compute stage values
  # stage 1
  ktmp = integrator.f(uprev,p,t)
  gtmp = integrator.g(uprev,p,t)

  # store gtmp for stage values Y^(k(j)j)
  gtmp1 = copy(gtmp)

  Y10 = ktmp*dt

  if typeof(W.dW) <: Number
    Y1j = gtmp*_dW
  else
    if is_diagonal_noise(integrator.sol.prob)
      Y1j = gtmp.*_dW
    else
      Y1j = gtmp.*_dW'
    end
  end

  # stage 2
  Y20 = uprev + a0021*Y10

  if typeof(W.dW) <: Number
    Y20 += a0j21*Y1j
    Y2j = integrator.g(uprev + aj021*Y10 + ajj21*Y1j,p,t)*_dW
  else
    if !is_diagonal_noise(integrator.sol.prob)
      Y2j = similar(integrator.sol.prob.noise_rate_prototype)
      for j = 1:m
        Y20 += a0j21*Y1j[:,j]
        tmpu = uprev + aj021*Y10 + ajj21*Y1j[:,j]
        Y2j[:,j] = integrator.g(tmpu,p,t)[:,j]*_dW[j]
      end
    else
      Y20 += a0j21*Y1j
      Y2j = similar(uprev)
      for j = 1:m
        tmpu = uprev + aj021*Y10
        tmpu[j] += ajj21*Y1j[j]
        Y2j[j] = integrator.g(tmpu,p,t)[j]*_dW[j]
      end
    end
  end
  ktmp = integrator.f(Y20,p,t)
  Y20 = ktmp*dt

  # stage 3
  Y30 = uprev + a0032*Y20
  if typeof(W.dW) <: Number
    Y30 += a0j31*Y1j + a0j32*Y2j
    Y3j = integrator.g(uprev + ajj31*Y1j + ajj32*Y2j,p,t)*_dW
  else
    if !is_diagonal_noise(integrator.sol.prob)
      Y3j = similar(integrator.sol.prob.noise_rate_prototype)
      for j = 1:m
        tmpu = uprev + ajj31*Y1j[:,j] + ajj32*Y2j[:,j]
        Y30 += a0j31*Y1j[:,j] + a0j32*Y2j[:,j]
        for l = 1:m
          if l!=j
            tmpu += (ajl31*Y1j[:,l] + ajl32*Y2j[:,l])
          end
        end
        Y3j[:,j] = integrator.g(tmpu,p,t)[:,j]*_dW[j]
      end
    else
      Y30 += a0j31*Y1j + a0j32*Y2j
      Y3j = similar(uprev)
      tmpu2 = uprev + ajl31*Y1j + ajl32*Y2j
      for j = 1:m
        tmpu = copy(tmpu2)
        tmpu[j] += ajj31*Y1j[j] + ajj32*Y2j[j] - (ajl31*Y1j[j] + ajl32*Y2j[j])
        Y3j[j] = integrator.g(tmpu,p,t)[j]*_dW[j]
      end
    end
  end
  ktmp = integrator.f(Y30,p,t)
  Y30 = ktmp*dt

  # stage 4
  Y40 = uprev + a0043*Y30
  if typeof(W.dW) <: Number
    Y40 += a0j41*Y1j
    Y4j = integrator.g(uprev + ajj41*Y1j + ajj42*Y2j + ajj43*Y3j + aj041*Y10,p,t)*_dW
  else
    if !is_diagonal_noise(integrator.sol.prob)
      Y4j = similar(integrator.sol.prob.noise_rate_prototype)
      for j = 1:m
        tmpu = uprev + ajj41*Y1j[:,j] + ajj42*Y2j[:,j] + ajj43*Y3j[:,j] + aj041*Y10
        Y40 += a0j41*Y1j[:,j]
        for l = 1:m
          if l!=j
            tmpu += (ajl41*Y1j[:,l] + ajl42*Y2j[:,l])
          end
        end
        Y4j[:,j] = integrator.g(tmpu,p,t)[:,j]*_dW[j]
    end
    else
      Y40 += a0j41*Y1j
      Y4j = similar(uprev)
      tmpu2 = uprev + ajl41*Y1j + ajl42*Y2j + aj041*Y10
      for j = 1:m
        tmpu = copy(tmpu2)
        tmpu[j] += ajj41*Y1j[j] + ajj42*Y2j[j] + ajj43*Y3j[j] - (ajl41*Y1j[j] + ajl42*Y2j[j])
        Y4j[j] = integrator.g(tmpu,p,t)[j]*_dW[j]
      end
    end
  end
  ktmp = integrator.f(Y40,p,t)
  Y40 = ktmp*dt

  # add stages together
  u = uprev + c01*Y10+c02*Y20+c03*Y30+c04*Y40

  if typeof(W.dW) <: Number
    u += cj1*Y1j + cj2*Y2j + cj3*Y3j + cj4*Y4j
  else
    if is_diagonal_noise(integrator.sol.prob)
      u += @. cj1*Y1j + cj2*Y2j + cj3*Y3j + cj4*Y4j
    else
      for j = 1:m
        u += @. cj1*Y1j[:,j] + cj2*Y2j[:,j] + cj3*Y3j[:,j] + cj4*Y4j[:,j]
        #add stage values for non-commutative processes, Y3^(k(j)j) and Y4^(k(j)j)
        for k = 1:m
          η2 = @view Ihat2[k,:]
          tmp = gtmp1*η2/(4*γ)
          Y3kj = integrator.sqdt*integrator.g(uprev + tmp,p,t)
          Y4kj = integrator.sqdt*integrator.g(uprev - tmp,p,t)
          u += @. γ*Y3kj[:,k] - γ*Y4kj[:,k]
        end
      end
    end
  end
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::NON2Cache,f=integrator.f)
  @unpack _dW,_dZ, chi1,Ihat2,tab,gtmp,gtmp1,ktmp,Y10,Y20,Y30,Y40,Y1j,Y2j,Y3j,Y4j,tmpu,tmpu2 = cache
  @unpack c01,c02,c03,c04,cj1,cj2,cj3,cj4,a0021,a0032,a0043,aj021,aj041,a0j21,a0j31,a0j32,a0j41,ajj21,ajj31,ajj32,ajj41,ajj42,ajj43,ajl31,ajl32,ajl41,ajl42,γ,NORMAL_ONESIX_QUANTILE = cache.tab
  @unpack t,dt,uprev,u,W,p = integrator

  m = length(W.dW)
  @.. _dW =  W.dW

  # define three-point distributed random variables
  @.. chi1 = W.dW / integrator.sqdt
  sq3dt = sqrt(3*dt)
  calc_threepoint_random!(_dW, sq3dt, NORMAL_ONESIX_QUANTILE, chi1)

  if !(typeof(W.dW) <: Number)
    # define two-point distributed random variables
    calc_twopoint_random!(_dZ,integrator.sqdt,W.dZ)
    for j = 1:m
      for l = 1:j-1
        Ihat2[j, l] = -_dZ[j]*_dW[l]/integrator.sqdt
        Ihat2[l, j] = _dW[l]*_dZ[j]/integrator.sqdt
      end
    end
  end

  # compute stage values
  # stage 1
  integrator.f(ktmp,uprev,p,t)
  integrator.g(gtmp,uprev,p,t)

  copyto!(gtmp1, gtmp)

  @.. Y10 = ktmp*dt

  if is_diagonal_noise(integrator.sol.prob)
    @.. Y1j = gtmp*_dW
  else
    for j = 1:m
      Y1jj = @view Y1j[:,j]
      gtmpj = @view gtmp[:,j]
      @.. Y1jj = gtmpj*_dW[j]
    end
  end

  # stage 2
  @.. Y20 = uprev + a0021*Y10

  if !is_diagonal_noise(integrator.sol.prob)
    fill!(Y2j, zero(eltype(integrator.sol.prob.noise_rate_prototype)))
    for j = 1:m
      Y1jj = @view Y1j[:,j]
      Y2jj = @view Y2j[:,j]
      @.. Y20 = Y20 + a0j21*Y1jj
      @.. tmpu = uprev + aj021*Y10 + ajj21*Y1jj
      integrator.g(gtmp,tmpu,p,t)
      gtmpj = @view gtmp[:,j]
      @.. Y2jj = gtmpj*_dW[j]
    end
  else
    fill!(Y2j, zero(eltype(integrator.u)))
    @.. Y20 = Y20 + a0j21*Y1j
    @.. tmpu = uprev + aj021*Y10
    for j = 1:m
      fill!(tmpu2,zero(eltype(integrator.u)))
      tmpu2[j] = ajj21*Y1j[j]
      @.. tmpu2 = tmpu + tmpu2
      integrator.g(gtmp,tmpu2,p,t)
      Y2j[j]= gtmp[j]*_dW[j]
    end
  end
  integrator.f(ktmp,Y20,p,t)
  @.. Y20 = ktmp*dt

  # stage 3
  @.. Y30 = uprev + a0032*Y20

  if !is_diagonal_noise(integrator.sol.prob)
    fill!(Y3j, zero(eltype(integrator.sol.prob.noise_rate_prototype)))
    for j = 1:m
      Y1jj = @view Y1j[:,j]
      Y2jj = @view Y2j[:,j]
      Y3jj = @view Y3j[:,j]
      @.. tmpu = uprev + ajj31*Y1jj + ajj32*Y2jj
      @.. Y30 = Y30 + a0j31*Y1jj + a0j32*Y2jj
      for l = 1:m
        if l!=j
          Y1jl = @view Y1j[:,l]
          Y2jl = @view Y2j[:,l]
          @.. tmpu = tmpu + (ajl31*Y1jl + ajl32*Y2jl)
        end
      end
      integrator.g(gtmp,tmpu,p,t)
      gtmpj = @view gtmp[:,j]
      @.. Y3jj = gtmpj*_dW[j]
    end
  else
    @.. Y30 = Y30 + a0j31*Y1j + a0j32*Y2j
    fill!(Y3j, zero(eltype(integrator.u)))
    @.. tmpu2 = uprev + ajl31*Y1j + ajl32*Y2j
    for j = 1:m
      @.. tmpu = tmpu2
      tmpu[j] += ajj31*Y1j[j] + ajj32*Y2j[j] - (ajl31*Y1j[j] + ajl32*Y2j[j])
      integrator.g(gtmp,tmpu,p,t)
      Y3j[j] = gtmp[j]*_dW[j]
    end
  end
  integrator.f(ktmp,Y30,p,t)
  @.. Y30 = ktmp*dt

  # stage 4
  @.. Y40 = uprev + a0043*Y30

  if !is_diagonal_noise(integrator.sol.prob)
    fill!(Y4j, zero(eltype(integrator.sol.prob.noise_rate_prototype)))
    for j = 1:m
      Y1jj = @view Y1j[:,j]
      Y2jj = @view Y2j[:,j]
      Y3jj = @view Y3j[:,j]
      Y4jj = @view Y4j[:,j]
      @.. tmpu = uprev + ajj41*Y1jj + ajj42*Y2jj + ajj43*Y3jj + aj041*Y10
      @.. Y40 = Y40 + a0j41*Y1jj
      for l = 1:m
        if l!=j
          Y1jl = @view Y1j[:,l]
          Y2jl = @view Y2j[:,l]
          @.. tmpu = tmpu + (ajl41*Y1jl + ajl42*Y2jl)
        end
      end
      integrator.g(gtmp,tmpu,p,t)
      gtmpj = @view gtmp[:,j]
      @.. Y4jj = gtmpj*_dW[j]
  end
  else
    @.. Y40 = Y40 + a0j41*Y1j
    fill!(Y4j, zero(eltype(integrator.u)))
    @.. tmpu2 = uprev + ajl41*Y1j + ajl42*Y2j + aj041*Y10
    for j = 1:m
      @.. tmpu = tmpu2
      tmpu[j] += ajj41*Y1j[j] + ajj42*Y2j[j] + ajj43*Y3j[j] - (ajl41*Y1j[j] + ajl42*Y2j[j])
      integrator.g(gtmp,tmpu,p,t)
      Y4j[j] = gtmp[j]*_dW[j]
    end
  end
  integrator.f(ktmp,Y40,p,t)
  @.. Y40 = ktmp*dt

  # add stages together
  @.. u = uprev + c01*Y10+c02*Y20+c03*Y30+c04*Y40

  if is_diagonal_noise(integrator.sol.prob)
    @. u = u + cj1*Y1j + cj2*Y2j + cj3*Y3j + cj4*Y4j
  else
    for j = 1:m
      Y1jj = @view Y1j[:,j]
      Y2jj = @view Y2j[:,j]
      Y3jj = @view Y3j[:,j]
      Y4jj = @view Y4j[:,j]
      @.. u = u + cj1*Y1jj + cj2*Y2jj + cj3*Y3jj + cj4*Y4jj

      #add stage values for non-commutative processes
      for k = 1:m
        η2 = @view Ihat2[k,:]
        gtmpk = @view gtmp[:,k]
        mul!(tmpu2,gtmp1,η2)
        @.. tmpu2 = tmpu2/(4*γ)
        @.. tmpu = uprev + tmpu2
        integrator.g(gtmp,tmpu,p,t)
        @.. u = u + γ*integrator.sqdt*gtmpk
        @.. tmpu = uprev - tmpu2
        integrator.g(gtmp,tmpu,p,t)
        @.. u = u - @. γ*integrator.sqdt*gtmpk
      end
    end
  end
end

# SIE / SME
@muladd function perform_step!(integrator,cache::SIESMEConstantCache,f=integrator.f)
  @unpack α1,α2,γ1,λ1,λ2,λ3,µ1,µ2,µ3,µ0,µbar0,λ0,λbar0,ν1,ν2,β2,β3,δ2,δ3 = cache
  @unpack t,dt,uprev,u,W,p = integrator

  W2 = @. (W.dW)^2/integrator.sqdt
  W3 = @. ν2*(W.dW)^3/dt

  # compute stage values
  k0 = integrator.f(uprev,p,t)
  g0 = integrator.g(uprev,p,t)

  # k1, g1, g2
  if is_diagonal_noise(integrator.sol.prob)
    k1 = integrator.f(uprev + λ0*k0*dt + ν1*g0.*W.dW + g0.*W3, p, t + µ0*dt)
    g1 = integrator.g(uprev + λbar0*k0*dt + β2*g0*integrator.sqdt + β3*g0.*W2, p, t + µbar0*dt)
    g2 = integrator.g(uprev + λbar0*k0*dt + δ2*g0*integrator.sqdt + δ3*g0.*W2, p, t + µbar0*dt)

  else
    # typeof(W.dW) <: Number
    k1 = integrator.f(uprev + λ0*k0*dt + ν1*g0*W.dW + g0*W3, p, t + µ0*dt)
    g1 = integrator.g(uprev + λbar0*k0*dt + β2*g0*integrator.sqdt + β3*g0*W2, p, t + µbar0*dt)
    g2 = integrator.g(uprev + λbar0*k0*dt + δ2*g0*integrator.sqdt + δ3*g0*W2, p, t + µbar0*dt)
  end

  # add stages together
  u = uprev + (α1*k0+α2*k1)*dt

  # add noise
  if typeof(W.dW) <: Number
    u += γ1*g0*W.dW + (λ1*W.dW + λ2*integrator.sqdt + λ3*W2)*g1 +  (µ1*W.dW + µ2*integrator.sqdt + µ3*W2)*g2
  else
    u += @. γ1*g0*W.dW + (λ1*W.dW + λ2*integrator.sqdt + λ3*W2)*g1 +  (µ1*W.dW + µ2*integrator.sqdt + µ3*W2)*g2
  end
  integrator.u = u
end

# SIE / SME
@muladd function perform_step!(integrator,cache::SIESMECache,f=integrator.f)
  @unpack W2,W3,tab,k0,k1,g0,g1,g2,tmpu = cache
  @unpack α1,α2,γ1,λ1,λ2,λ3,µ1,µ2,µ3,µ0,µbar0,λ0,λbar0,ν1,ν2,β2,β3,δ2,δ3 = cache.tab
  @unpack t,dt,uprev,u,W,p = integrator

  @. W2 = (W.dW)^2/integrator.sqdt
  @. W3 = ν2*(W.dW)^3/dt

  # compute stage values
  integrator.f(k0,uprev,p,t)
  integrator.g(g0,uprev,p,t)

  # k1, g1, g2
  if is_diagonal_noise(integrator.sol.prob)
    @.. tmpu = uprev + λ0*k0*dt + ν1*g0*W.dW + g0*W3
    integrator.f(k1, tmpu, p, t + µ0*dt)
    @.. tmpu = uprev + λbar0*k0*dt + β2*g0*integrator.sqdt + β3*g0*W2
    integrator.g(g1,tmpu, p, t + µbar0*dt)
    @.. tmpu = uprev + λbar0*k0*dt + δ2*g0*integrator.sqdt + δ3*g0*W2
    integrator.g(g2,tmpu, p, t + µbar0*dt)
  end

  # add stages together
  @.. u = uprev + (α1*k0+α2*k1)*dt

  # add noise
  if is_diagonal_noise(integrator.sol.prob)
     @. u = u + γ1*g0*W.dW + (λ1*W.dW + λ2*integrator.sqdt + λ3*W2)*g1 +  (µ1*W.dW + µ2*integrator.sqdt + µ3*W2)*g2
  end
end


text:
using Winston

x = linspace(pi, 3pi, 60)
c = cos(x)
s = sin(x)

p = FramedPlot(aspect_ratio=1)
setattr(p.frame1, draw_grid=true, tickdir=1)

setattr(p.x1, label="bottom", subticks=1)
setattr(p.y1, label="left", draw_spine=false)
setattr(p.x2, label="top", range=(10,1000), log=true)

setattr(p.y2, label="right", draw_ticks=false,
    ticklabels=["-1", "-1/2", "0", "1/2", "1"])

add(p, Curve(x, c, kind="dash"))
add(p, Curve(x, s))

file(p, "example6.png")


text:
BUILD_SHARED_LIBS = true
include("build.jl")

max_stars_repo_name: chemfiles/chemfiles-prebuilt

max_issues_repo_name: chemfiles/chemfiles-prebuilt

max_forks_repo_name: chemfiles/chemfiles-prebuilt