#pragma OPENCL EXTENSION cl_khr_global_int32_base_atomics : enable
#pragma OPENCL EXTENSION cl_khr_local_int32_base_atomics : enable
#pragma OPENCL EXTENSION cl_amd_printf : enable

uint TausStep(uint* z, int S1, int S2, int S3, uint M)
{
    uint b = (((*z << S1) ^ *z) >> S2);
    return *z = (((*z & M) << S3) ^ b);
}

uint LCGStep(uint* z, uint A, uint C)
{
    return *z = (A * *z + C);
}

float rand2(uint2* zp)
{
    uint z0 = (*zp).s0;
    uint z1 = (*zp).s1;
    float ret = 2.3283064365387e-10 * (TausStep(&z0, 13, 19, 12, 4294967294UL)
                                       ^ LCGStep(&z1, 1664525, 1013904223UL));
    (*zp).s0 = z0;
    (*zp).s1 = z1;
    return ret;
}

#undef record
#define record(flags, data) \
if(_record_flags & (flags)) \
{ \
    int idx = atomic_inc(&_local_record_idx); \
    if(idx >= _record_buffer_size) \
    { \
        _error_buffer[i + 1] = 1; \
        _local_record_idx = _record_buffer_size - 1; \
    } \
    float4 record; \
    record.s0 = t; \
    record.s1 = (data); \
    record.s2 = flags; \
    record.s3 = i; \
    _record_buffer[idx + _local_record_idx_start] = record; \
} \


__kernel void RS_step
(
    const float _t,
    __global int* _error_buffer,
    __global int* _record_flags_buffer,
    __global int* _record_idx,
    __global int* _record_idx_start,
    __global float4* _record_buffer,
    const int _record_buffer_size,
    __global float* _rwvalues_0_buf,
    __global float* _rwvalues_1_buf,
    __global float* _rwvalues_2_buf,
    __global float* _rwvalues_3_buf,
    __global float* _rwvalues_4_buf,
    __global float* _rwvalues_5_buf,
    __global float* _rwvalues_6_buf,
    __global float* _rwvalues_7_buf,
    __global float* _rwvalues_8_buf,
    __global float* _rwvalues_9_buf,
    __global float* _rwvalues_10_buf,

    __global float* _rovalues_0_buf,
    __global float* _rovalues_1_buf,
    __global float* _rovalues_2_buf,
    __global float* _rovalues_3_buf,
    __global float* _rovalues_4_buf,
    __global float* _rovalues_5_buf,
    __global float* _rovalues_6_buf,
    __global float* _rovalues_7_buf,
    __global float* _rovalues_8_buf,
    __global float* _rovalues_9_buf,
    __global float* _rovalues_10_buf,
    __global float* _rovalues_11_buf,
    __global float* _rovalues_12_buf,
    __global float* _rovalues_13_buf,
    __global float* _rovalues_14_buf,

    const float stdp_Ap,
    const float stdp_Am,
    const float stdp_o1_tau,
    const float stdp_o2_tau,
    const float stdp_r1_tau,
    const float stdp_r2_tau,
    __global uint2* _rand_state_buf,
    __global float* _dt_buf,
    __global int* _circ_buffer_start,
    __global int* _circ_buffer_end,
    __global float* _circ_buffer,
    __global int* _fired_syn_idx_buffer,
    __global int* _fired_syn_buffer,
    __global float* _static_excitatory_on_buf,
    __global float* _static_excitatory_weight_buf,
    __global float* _static_inhibitory_weight_buf,
    __global float* _stdp_weight_buf,
    __global float* _stdp_r1_buf,
    __global float* _stdp_r2_buf,
    __global float* _stdp_last_r_time_buf,
    __local int* _syn_thresh_0_arr,
    const int count
)
{
    int i = get_global_id(0);
    int _local_id = get_local_id(0);
    int _group_id = get_group_id(0);
    const float V_tol = 2.000000e-01;
    const float u_tol = 2.000000e-02;
    const float axon_delay = 1.000000e+00;
    const float rand_syn_r_tol = 1.000000e-02;
    const float static_excitatory_E = 0.000000e+00;
    const float static_excitatory_tau1 = 2.000000e+00;
    const float static_excitatory_tau2 = 5.000000e+01;
    const float static_excitatory_N = 1.000000e+01;
    const float static_excitatory_P = 2.500000e-01;
    const float static_inhibitory_E = -8.000000e+01;
    const float static_inhibitory_gsyn = 2.000000e+00;
    const float static_inhibitory_tau = 1.000000e+01;
    const float static_inhibitory_N = 1.000000e+01;
    const float static_inhibitory_P = 2.500000e-01;
    const float stdp_s_tol = 1.000000e-02;
    const float stdp_E = -8.000000e+01;
    const float stdp_gsyn = 2.000000e-01;
    const float stdp_tau = 1.000000e+01;
    const float stdp_N = 1.000000e+01;
    const float stdp_P = 2.500000e-01;
    const float static_excitatory_s1_tol = 0.1;
    const float static_excitatory_s2_tol = 0.1;
    const float static_inhibitory_s_tol = 0.1;
    __local int _syn_thresh_0_num;
    int _local_size = get_local_size(0);
    __local int _num_complete;

    if(_local_id == 0)
        _num_complete = 0;

    float _cur_time = 0;
    const float timestep = 1;
    int _record_flags = _record_flags_buffer[i];
    __local int _local_record_idx;
    __local int _local_record_idx_start;

    if(_local_id == 0)
    {
        _local_record_idx = _record_idx[_group_id];
        _local_record_idx_start = _record_idx_start[_group_id];
    }

    float _dt;
    float t = _t;
    float _dt_residual = 0;
    _dt = _dt_buf[i];
    float _rwvalues_0 = _rwvalues_0_buf[i];
    float V = _rwvalues_0;
    float _rwvalues_1 = _rwvalues_1_buf[i];
    float u = _rwvalues_1;
    float _rwvalues_2 = _rwvalues_2_buf[i];
    float rand_syn_r = _rwvalues_2;
    float _rwvalues_3 = _rwvalues_3_buf[i];
    float static_excitatory_s1 = _rwvalues_3;
    float _rwvalues_4 = _rwvalues_4_buf[i];
    float static_excitatory_s2 = _rwvalues_4;
    float _rwvalues_5 = _rwvalues_5_buf[i];
    float static_inhibitory_s = _rwvalues_5;
    float _rwvalues_6 = _rwvalues_6_buf[i];
    float stdp_s = _rwvalues_6;
    float _rwvalues_7 = _rwvalues_7_buf[i];
    float stdp_o1 = _rwvalues_7;
    float _rwvalues_8 = _rwvalues_8_buf[i];
    float stdp_o2 = _rwvalues_8;
    float _rwvalues_9 = _rwvalues_9_buf[i];
    float stdp_last_o_time = _rwvalues_9;
    float _rwvalues_10 = _rwvalues_10_buf[i];
    float stdp_active = _rwvalues_10;
    float _rovalues_0 = _rovalues_0_buf[i];
    float a = _rovalues_0;
    float _rovalues_1 = _rovalues_1_buf[i];
    float b = _rovalues_1;
    float _rovalues_2 = _rovalues_2_buf[i];
    float c = _rovalues_2;
    float _rovalues_3 = _rovalues_3_buf[i];
    float d = _rovalues_3;
    float _rovalues_4 = _rovalues_4_buf[i];
    float start_amp = _rovalues_4;
    float _rovalues_5 = _rovalues_5_buf[i];
    float end_amp = _rovalues_5;
    float _rovalues_6 = _rovalues_6_buf[i];
    float start_t = _rovalues_6;
    float _rovalues_7 = _rovalues_7_buf[i];
    float end_t = _rovalues_7;
    float _rovalues_8 = _rovalues_8_buf[i];
    float rand_syn_tau = _rovalues_8;
    float _rovalues_9 = _rovalues_9_buf[i];
    float rand_syn_weight = _rovalues_9;
    float _rovalues_10 = _rovalues_10_buf[i];
    float rand_syn_rate = _rovalues_10;
    float _rovalues_11 = _rovalues_11_buf[i];
    float rand_syn_frequency = _rovalues_11;
    float _rovalues_12 = _rovalues_12_buf[i];
    float rand_syn_sine_amp = _rovalues_12;
    float _rovalues_13 = _rovalues_13_buf[i];
    float static_excitatory_gsyn1 = _rovalues_13;
    float _rovalues_14 = _rovalues_14_buf[i];
    float static_excitatory_gsyn2 = _rovalues_14;
    uint2 _rand_state = _rand_state_buf[i];
    const int _syn_offset = 0 + i * 1250;
    int _syn_table_end = _fired_syn_idx_buffer[i + 0];

    if(_syn_table_end != _syn_offset)
    {
        for(int _syn_table_idx = _syn_offset; _syn_table_idx < _syn_table_end; _syn_table_idx++)
        {
            int syn_i = _fired_syn_buffer[_syn_table_idx];

            if(syn_i < 700)
            {
                int _g_syn_i = syn_i - 0 + i * 700;
                /* Load values */
                float static_excitatory_on = _static_excitatory_on_buf[_g_syn_i];
                float static_excitatory_weight = _static_excitatory_weight_buf[_g_syn_i];
                /* Syn code */
                float frac = (sqrt(-2 * log(rand2(&_rand_state) + 0.000001)) * cospi(2 * rand2(&_rand_state))) * sqrt(static_excitatory_N * static_excitatory_P * (1.0 - static_excitatory_P)) + static_excitatory_N * static_excitatory_P;
                frac = fmax(frac, 0.0f);
                static_excitatory_s1 += static_excitatory_on * frac * static_excitatory_gsyn1 * static_excitatory_weight;
                static_excitatory_s2 += static_excitatory_on * frac * static_excitatory_gsyn2 * static_excitatory_weight;
                //record(2, static_excitatory_weight, 100+syn_i); //flags, val, tag
                /* Save values */
            }
            else if(syn_i < 750)
            {
                int _g_syn_i = syn_i - 700 + i * 50;
                /* Load values */
                float static_inhibitory_weight = _static_inhibitory_weight_buf[_g_syn_i];
                /* Syn code */
                float frac = (sqrt(-2 * log(rand2(&_rand_state) + 0.000001)) * cospi(2 * rand2(&_rand_state))) * sqrt(static_inhibitory_N * static_inhibitory_P * (1.0 - static_inhibitory_P)) + static_inhibitory_N * static_inhibitory_P;
                frac = fmax(frac, 0.0f);
                static_inhibitory_s += frac * static_inhibitory_gsyn * static_inhibitory_weight;
                //record(2, static_inhibitory_weight, 100+syn_i); //flags, val, tag
                /* Save values */
            }
            else if(syn_i < 1250)
            {
                int _g_syn_i = syn_i - 750 + i * 500;
                /* Load values */
                float stdp_weight = _stdp_weight_buf[_g_syn_i];
                float stdp_r1 = _stdp_r1_buf[_g_syn_i];
                float stdp_r2 = _stdp_r2_buf[_g_syn_i];
                float stdp_last_r_time = _stdp_last_r_time_buf[_g_syn_i];
                /* Syn code */
                stdp_o1 = stdp_o1 * exp(-(t - stdp_last_o_time) / stdp_o1_tau);
                stdp_o2 = stdp_o2 * exp(-(t - stdp_last_o_time) / stdp_o2_tau);
                stdp_last_o_time = t;

                if(stdp_active > 0)
                    stdp_weight -= stdp_Am * (stdp_o1 - stdp_o2);

                stdp_weight = fmax(stdp_weight, 0.0f);
                float frac = (sqrt(-2 * log(rand2(&_rand_state) + 0.000001)) * cospi(2 * rand2(&_rand_state))) * sqrt(stdp_N * stdp_P * (1.0 - stdp_P)) + stdp_N * stdp_P;
                frac = fmax(frac, 0.0f);
                stdp_s += frac * stdp_gsyn * stdp_weight;
                stdp_r1 = stdp_r1 * exp(-(t - stdp_last_r_time) / stdp_r1_tau) + 1;
                stdp_r2 = stdp_r2 * exp(-(t - stdp_last_r_time) / stdp_r2_tau) + 1;
                stdp_last_r_time = t;
                /* Save values */
                _stdp_weight_buf[_g_syn_i] = stdp_weight;
                _stdp_r1_buf[_g_syn_i] = stdp_r1;
                _stdp_r2_buf[_g_syn_i] = stdp_r2;
                _stdp_last_r_time_buf[_g_syn_i] = stdp_last_r_time;
            }
        }

        _dt = 0.1f;
        _fired_syn_idx_buffer[i + 0] = _syn_offset;
    }

    {
        if(rand_syn_weight > 0)
        {
            if(rand2(&_rand_state) < timestep * (rand_syn_rate + rand_syn_sine_amp * sin(2.0 * M_PI_F * t / 1000.0 * rand_syn_frequency)) / 1000.0)
                rand_syn_r += rand_syn_weight;

            _dt = 0.1f;
        }
    }

    barrier(CLK_LOCAL_MEM_FENCE);

    while(_num_complete < _local_size)
    {
        bool _any_thresh = false;
        /* Threshold statuses */
        bool thresh_0_state = false;
        bool thresh_1_state = false;

        if(_local_id == 0)
        {
            _syn_thresh_0_num = 0;
        }

        barrier(CLK_LOCAL_MEM_FENCE);

        while(!_any_thresh && _cur_time < timestep)
        {
            t = _t + _cur_time;

            /* Post-thresh integrator code */
            /* Clamp the _dt not too overshoot the timestep */
            if(_cur_time < timestep && _cur_time + _dt >= timestep)
            {
                _dt_residual = _dt;
                _dt = timestep - _cur_time + 0.0001f;
                _dt_residual -= _dt;
            }

            float _error = 0;
            /* See where the thresholded states are before changing them (doesn't work for synapse states)*/
            bool thresh_0_pre_state = V > 0;
            bool thresh_1_pre_state = V > 0;
            bool syn_thresh_0_pre_state = V > 0;
            /* Declare local variables */
            float I;
            /* Pre-stage code */
            {
                record(1, V); //flags, val, tag
            }
            {
                if(t == start_t)
                    _dt = 0.1f;
                else if(t == end_t)
                    _dt = 0.1f;
            }
            /* Integrator code */
            /* Declare storage for first state estimate*/
            float _V_0 = V;
            float _u_0 = u;
            float _rand_syn_r_0 = rand_syn_r;
            float _static_excitatory_s1_0 = static_excitatory_s1;
            float _static_excitatory_s2_0 = static_excitatory_s2;
            float _static_inhibitory_s_0 = static_inhibitory_s;
            float _stdp_s_0 = stdp_s;
            /* First derivative stage */
            float _dV_dt_1;
            float _du_dt_1;
            float _drand_syn_r_dt_1;
            float _dstatic_excitatory_s1_dt_1;
            float _dstatic_excitatory_s2_dt_1;
            float _dstatic_inhibitory_s_dt_1;
            float _dstdp_s_dt_1;
            /* Second derivative stage */
            float _dV_dt_2;
            float _du_dt_2;
            float _drand_syn_r_dt_2;
            float _dstatic_excitatory_s1_dt_2;
            float _dstatic_excitatory_s2_dt_2;
            float _dstatic_inhibitory_s_dt_2;
            float _dstdp_s_dt_2;
            /* Compute the first derivatives */
            {
                I = 0;
            }
            {
                if(t >= start_t && t < end_t)
                {
                    I += (t - start_t) / (end_t - start_t) * (end_amp - start_amp) + start_amp;
                }
            }
            {
                I += rand_syn_r * (0 - V);
            }
            {
                I += (static_excitatory_s1 + static_excitatory_s2 * 1.0f / (1.0f + 1.0f / 3.57 * exp(-0.062 * V))) * (static_excitatory_E - V);
            }
            {
                I += static_inhibitory_s * (static_inhibitory_E - V);
            }
            {
                I += stdp_s * (stdp_E - V);
            }
            {
                _dV_dt_1 = (0.04f * V + 5) * V + 140 - u + I;
                _du_dt_1 = a * (b * V - u);
            }
            {
                _drand_syn_r_dt_1 = -rand_syn_r / rand_syn_tau;
            }
            {
                _dstatic_excitatory_s1_dt_1 = -static_excitatory_s1 / static_excitatory_tau1;
                _dstatic_excitatory_s2_dt_1 = -static_excitatory_s2 / static_excitatory_tau2;
            }
            {
                _dstatic_inhibitory_s_dt_1 = -static_inhibitory_s / static_inhibitory_tau;
            }
            {
                _dstdp_s_dt_1 = -stdp_s / stdp_tau;
            }
            /* Compute the first state estimate */
            _dV_dt_1 *= _dt;
            _V_0 += _dV_dt_1;
            _du_dt_1 *= _dt;
            _u_0 += _du_dt_1;
            _drand_syn_r_dt_1 *= _dt;
            _rand_syn_r_0 += _drand_syn_r_dt_1;
            _dstatic_excitatory_s1_dt_1 *= _dt;
            _static_excitatory_s1_0 += _dstatic_excitatory_s1_dt_1;
            _dstatic_excitatory_s2_dt_1 *= _dt;
            _static_excitatory_s2_0 += _dstatic_excitatory_s2_dt_1;
            _dstatic_inhibitory_s_dt_1 *= _dt;
            _static_inhibitory_s_0 += _dstatic_inhibitory_s_dt_1;
            _dstdp_s_dt_1 *= _dt;
            _stdp_s_0 += _dstdp_s_dt_1;
            /* Compute the derivatives again */
            {
                I = 0;
            }
            {
                if(t >= start_t && t < end_t)
                {
                    I += (t - start_t) / (end_t - start_t) * (end_amp - start_amp) + start_amp;
                }
            }
            {
                I += _rand_syn_r_0 * (0 - _V_0);
            }
            {
                I += (_static_excitatory_s1_0 + _static_excitatory_s2_0 * 1.0f / (1.0f + 1.0f / 3.57 * exp(-0.062 * _V_0))) * (static_excitatory_E - _V_0);
            }
            {
                I += _static_inhibitory_s_0 * (static_inhibitory_E - _V_0);
            }
            {
                I += _stdp_s_0 * (stdp_E - _V_0);
            }
            {
                _dV_dt_2 = (0.04f * _V_0 + 5) * _V_0 + 140 - _u_0 + I;
                _du_dt_2 = a * (b * _V_0 - _u_0);
            }
            {
                _drand_syn_r_dt_2 = -_rand_syn_r_0 / rand_syn_tau;
            }
            {
                _dstatic_excitatory_s1_dt_2 = -_static_excitatory_s1_0 / static_excitatory_tau1;
                _dstatic_excitatory_s2_dt_2 = -_static_excitatory_s2_0 / static_excitatory_tau2;
            }
            {
                _dstatic_inhibitory_s_dt_2 = -_static_inhibitory_s_0 / static_inhibitory_tau;
            }
            {
                _dstdp_s_dt_2 = -_stdp_s_0 / stdp_tau;
            }
            /* Compute the final change in state */
            _dV_dt_2 = (_dV_dt_1 + _dt * _dV_dt_2) / 2;
            _du_dt_2 = (_du_dt_1 + _dt * _du_dt_2) / 2;
            _drand_syn_r_dt_2 = (_drand_syn_r_dt_1 + _dt * _drand_syn_r_dt_2) / 2;
            _dstatic_excitatory_s1_dt_2 = (_dstatic_excitatory_s1_dt_1 + _dt * _dstatic_excitatory_s1_dt_2) / 2;
            _dstatic_excitatory_s2_dt_2 = (_dstatic_excitatory_s2_dt_1 + _dt * _dstatic_excitatory_s2_dt_2) / 2;
            _dstatic_inhibitory_s_dt_2 = (_dstatic_inhibitory_s_dt_1 + _dt * _dstatic_inhibitory_s_dt_2) / 2;
            _dstdp_s_dt_2 = (_dstdp_s_dt_1 + _dt * _dstdp_s_dt_2) / 2;
            /* Compute the error in this step */
            _error = max(_error, fabs(_dV_dt_1 - _dV_dt_2) / V_tol);
            _error = max(_error, fabs(_du_dt_1 - _du_dt_2) / u_tol);
            _error = max(_error, fabs(_drand_syn_r_dt_1 - _drand_syn_r_dt_2) / rand_syn_r_tol);
            _error = max(_error, fabs(_dstatic_excitatory_s1_dt_1 - _dstatic_excitatory_s1_dt_2) / static_excitatory_s1_tol);
            _error = max(_error, fabs(_dstatic_excitatory_s2_dt_1 - _dstatic_excitatory_s2_dt_2) / static_excitatory_s2_tol);
            _error = max(_error, fabs(_dstatic_inhibitory_s_dt_1 - _dstatic_inhibitory_s_dt_2) / static_inhibitory_s_tol);
            _error = max(_error, fabs(_dstdp_s_dt_1 - _dstdp_s_dt_2) / stdp_s_tol);
            /* Update state */
            V += _dV_dt_2;
            u += _du_dt_2;
            rand_syn_r += _drand_syn_r_dt_2;
            static_excitatory_s1 += _dstatic_excitatory_s1_dt_2;
            static_excitatory_s2 += _dstatic_excitatory_s2_dt_2;
            static_inhibitory_s += _dstatic_inhibitory_s_dt_2;
            stdp_s += _dstdp_s_dt_2;
            /* Advance and compute the new step size*/
            _cur_time += _dt;

            if(_error == 0)
                _dt = timestep;
            else
            {
                /* Approximate the cube root using Halley's Method (error is usually between 0 and 10)*/
                float cr = (1.0f + 2 * _error) / (2.0f + _error);
                float cr3 = cr * cr * cr;
                cr = cr * (cr3 + 2.0f * _error) / (2.0f * cr3 + _error);
                _dt *= 0.9f / cr;
                /* _dt *= 0.9f * rootn(_error, -3.0f); */
            }

            /* Detect thresholds */
            thresh_0_state = !thresh_0_pre_state && (V > 0);
            _any_thresh |= thresh_0_state;
            thresh_1_state = !thresh_1_pre_state && (V > 0);
            _any_thresh |= thresh_1_state;

            if(!syn_thresh_0_pre_state && (V > 0))
            {
                _any_thresh = true;
                _syn_thresh_0_arr[atomic_inc(&_syn_thresh_0_num)] = i;
            }

            /* Check exit condition */
            if(_cur_time >= timestep)
                atomic_inc(&_num_complete);
        }

        /* Handle thresholds */
        if(thresh_0_state)
        {
            float delay = 1.0f;
            V = c;
            u += d;
            delay = axon_delay;
            record(4, 0);
            _dt = 0.1f;
            int _idx_idx = 1 * i + 0;
            int _buff_start = _circ_buffer_start[_idx_idx];

            if(_buff_start != _circ_buffer_end[_idx_idx])
            {
                const int _circ_buffer_size = 150;
                int _end_idx;

                if(_buff_start < 0) //It is empty
                {
                    _circ_buffer_start[_idx_idx] = 0;
                    _circ_buffer_end[_idx_idx] = 1;
                    _end_idx = 1;
                }
                else
                {
                    _end_idx = _circ_buffer_end[_idx_idx] = (_circ_buffer_end[_idx_idx] + 1) % _circ_buffer_size;
                }

                int _buff_idx = (i * 1 + 0) * _circ_buffer_size + _end_idx - 1;
                _circ_buffer[_buff_idx] = t + delay;
            }
            else //It is full, error
            {
                _error_buffer[i + 1] = 2 + 0;
            }
        }

        if(thresh_1_state)
        {
            stdp_o1 = stdp_o1 * exp(-(t - stdp_last_o_time) / stdp_o1_tau) + 1;
            stdp_o2 = stdp_o2 * exp(-(t - stdp_last_o_time) / stdp_o2_tau) + 1;
            stdp_last_o_time = t;
        }

        barrier(CLK_LOCAL_MEM_FENCE);

        if(_syn_thresh_0_num > 0)
        {
            for(int _ii = 0; _ii < _syn_thresh_0_num; _ii++)
            {
                int nrn_id = _syn_thresh_0_arr[_ii];
                /* Declare locals */
                __local float stdp_s_local;
                __local float stdp_o1_local;
                __local float stdp_o2_local;
                __local float stdp_active_local;

                /* Init locals */
                if(i == nrn_id)
                {
                    stdp_s_local = stdp_s;
                    stdp_o1_local = stdp_o1;
                    stdp_o2_local = stdp_o2;
                    stdp_active_local = stdp_active;
                }

                barrier(CLK_LOCAL_MEM_FENCE);
                int _syn_offset = nrn_id * 500;

                for(int _g_syn_i = _syn_offset + _local_id; _g_syn_i < 500 + _syn_offset; _g_syn_i += _local_size)
                {
                    /* Load syn globals */
                    float stdp_weight = _stdp_weight_buf[_g_syn_i];
                    float stdp_r1 = _stdp_r1_buf[_g_syn_i];
                    float stdp_r2 = _stdp_r2_buf[_g_syn_i];
                    float stdp_last_r_time = _stdp_last_r_time_buf[_g_syn_i];
                    /* Thresh source */
                    stdp_r1 = stdp_r1 * exp(-(t - stdp_last_r_time) / stdp_r1_tau);
                    stdp_r2 = stdp_r2 * exp(-(t - stdp_last_r_time) / stdp_r2_tau);
                    stdp_last_r_time = t;

                    if(stdp_active_local > 0)
                        stdp_weight += stdp_Ap * (stdp_r1 - stdp_r2) * (stdp_o1_local - 1.0 - (stdp_o2_local - 1.0)); // Take the o'stdp_s_local from before the update

                    /* Save syn globals */
                    _stdp_weight_buf[_g_syn_i] = stdp_weight;
                    _stdp_r1_buf[_g_syn_i] = stdp_r1;
                    _stdp_r2_buf[_g_syn_i] = stdp_r2;
                    _stdp_last_r_time_buf[_g_syn_i] = stdp_last_r_time;
                }
            }
        }
    }

    if(_dt_residual > 0.1f)
        _dt = _dt_residual;

    if(_dt > timestep)
        _dt = timestep;

    _dt_buf[i] = _dt;
    _rwvalues_0_buf[i] = (float)(V);
    _rwvalues_1_buf[i] = (float)(u);
    _rwvalues_2_buf[i] = (float)(rand_syn_r);
    _rwvalues_3_buf[i] = (float)(static_excitatory_s1);
    _rwvalues_4_buf[i] = (float)(static_excitatory_s2);
    _rwvalues_5_buf[i] = (float)(static_inhibitory_s);
    _rwvalues_6_buf[i] = (float)(stdp_s);
    _rwvalues_7_buf[i] = (float)(stdp_o1);
    _rwvalues_8_buf[i] = (float)(stdp_o2);
    _rwvalues_9_buf[i] = (float)(stdp_last_o_time);
    _rwvalues_10_buf[i] = (float)(stdp_active);
    _rand_state_buf[i] = _rand_state;
    barrier(CLK_LOCAL_MEM_FENCE);

    if(_local_id == 0)
    {
        _record_idx[_group_id] = _local_record_idx;
    }
}