
/* Copyright 2012 Evan Miller */
/* All rights reserved. */

#ifdef cl_khr_fp64
typedef double   fcl_fp;
typedef double4  fcl_fp4;
typedef double16 fcl_fp16;
#else
typedef float    fcl_fp;
typedef float4   fcl_fp4;
typedef float16  fcl_fp16;
#endif

#ifndef M_SQRT1_2
#define M_SQRT1_2 M_SQRT1_2_F
#endif

#ifndef M_LN2
#define M_LN2 M_LN2_F
#endif

#ifdef cl_khr_fp64
#define FCL_1_OVER_SQRT_2_PI 0.3989422804014326
#else
#define FCL_1_OVER_SQRT_2_PI 0.3989422804014326f
#endif

#ifdef cl_khr_fp64
#define FCL_DIGAMMA_MIN_ARG_FOR_POLYNOMIAL_APPROX   8.5
#define FCL_DIGAMMA_MAX_ARG_FOR_LINEAR_APPROX       0.00001
#define FCL_DIGAMMA_D1                             -0.5772156649
#define FCL_DIGAMMA_S3                              0.08333333333
#define FCL_DIGAMMA_S4                              0.0083333333333
#define FCL_DIGAMMA_S5                              0.003968253968
#else
#define FCL_DIGAMMA_MIN_ARG_FOR_POLYNOMIAL_APPROX   8.5f
#define FCL_DIGAMMA_MAX_ARG_FOR_LINEAR_APPROX       0.00001f
#define FCL_DIGAMMA_D1                             -0.5772156649f
#define FCL_DIGAMMA_S3                              0.08333333333f
#define FCL_DIGAMMA_S4                              0.0083333333333f
#define FCL_DIGAMMA_S5                              0.003968253968f
#endif

#ifdef cl_khr_fp64
#define FCL_TRIGAMMA_MIN_ARG_FOR_POLYNOMIAL_APPROX  5.
#define FCL_TRIGAMMA_MAX_ARG_FOR_LINEAR_APPROX      0.0001
#define FCL_TRIGAMMA_B2                             0.1666666667
#define FCL_TRIGAMMA_B4                            -0.03333333333
#define FCL_TRIGAMMA_B6                             0.02380952381
#define FCL_TRIGAMMA_B8                            -0.03333333333
#else
#define FCL_TRIGAMMA_MIN_ARG_FOR_POLYNOMIAL_APPROX  5.f
#define FCL_TRIGAMMA_MAX_ARG_FOR_LINEAR_APPROX      0.0001f
#define FCL_TRIGAMMA_B2                             0.1666666667f
#define FCL_TRIGAMMA_B4                            -0.03333333333f
#define FCL_TRIGAMMA_B6                             0.02380952381f
#define FCL_TRIGAMMA_B8                            -0.03333333333f
#endif

#define fcl_slice_to_sliver(i) (4*i)

fcl_fp sum4(fcl_fp4 vec);
fcl_fp sum16(fcl_fp16 vec);
fcl_fp dot16(fcl_fp16 vec1, fcl_fp16 vec2);

fcl_fp4 select4(fcl_fp4 c, fcl_fp4 b, int4 a);
fcl_fp16 select16(fcl_fp16 c, fcl_fp16 b, int16 a);

fcl_fp fcl_logistic_cdf(fcl_fp alpha);
fcl_fp fcl_normal_cdf(fcl_fp alpha);

fcl_fp fcl_normal_pdf(fcl_fp alpha);
fcl_fp fcl_inverse_mills_ratio(fcl_fp alpha);

fcl_fp digamma(fcl_fp alpha);
fcl_fp16 digamma16(fcl_fp16 alpha);

fcl_fp trigamma(fcl_fp alpha);
fcl_fp16 trigamma16(fcl_fp16 alpha);

int4 convert_fp4_to_int4(fcl_fp4 res_fp);
int16 convert_fp16_to_int16(fcl_fp16 res_fp);

fcl_fp4 convert_int4_to_fp4(int4 res_i);
fcl_fp16 convert_int16_to_fp16(int16 res_i);

int4 int_array_elements4(__global int *array, int4 idx);
int16 int_array_elements16(__global int *array, int16 idx);
fcl_fp4 array_elements4(__global fcl_fp *array, int4 idx);
fcl_fp16 array_elements16(__global fcl_fp *array, int16 idx);

fcl_fp dot16(fcl_fp16 vec1, fcl_fp16 vec2) {
    fcl_fp result = 0.f;
    result += dot(vec1.lo.lo, vec2.lo.lo);
    result += dot(vec1.lo.hi, vec2.lo.hi);
    result += dot(vec1.hi.lo, vec2.hi.lo);
    result += dot(vec1.hi.hi, vec2.hi.hi);
    return result;
}

fcl_fp4 select4(fcl_fp4 c, fcl_fp4 b, int4 a) {
    fcl_fp4 res_fp;
    res_fp.s0 = a.s0 ? b.s0 : c.s0;
    res_fp.s1 = a.s1 ? b.s1 : c.s1;
    res_fp.s2 = a.s2 ? b.s2 : c.s2;
    res_fp.s3 = a.s3 ? b.s3 : c.s3;
    return res_fp;
}

fcl_fp16 select16(fcl_fp16 c, fcl_fp16 b, int16 a) {
    fcl_fp16 res_fp;
    res_fp.lo.lo = select4(c.lo.lo, b.lo.lo, a.lo.lo);
    res_fp.lo.hi = select4(c.lo.hi, b.lo.hi, a.lo.hi);
    res_fp.hi.lo = select4(c.hi.lo, b.hi.lo, a.hi.lo);
    res_fp.hi.hi = select4(c.hi.hi, b.hi.hi, a.hi.hi);
    return res_fp;
}

fcl_fp sum4(fcl_fp4 vec) {
    return vec.x + vec.y + vec.z + vec.w;
}

fcl_fp sum16(fcl_fp16 vec) {
    return sum4(vec.lo.lo) + sum4(vec.lo.hi) + sum4(vec.hi.lo) + sum4(vec.hi.hi);
}

fcl_fp fcl_logistic_cdf(fcl_fp alpha) {
    return 1.f / (1.f + exp(-alpha));
}

fcl_fp fcl_normal_cdf(fcl_fp alpha) {
    return 0.5f * erfc(-alpha * M_SQRT1_2);
}

fcl_fp fcl_normal_pdf(fcl_fp alpha) {
    return FCL_1_OVER_SQRT_2_PI * exp(-0.5f * alpha * alpha);
}

fcl_fp fcl_inverse_mills_ratio(fcl_fp alpha) {
    return fcl_normal_pdf(alpha) / fcl_normal_cdf(alpha);
}

fcl_fp digamma(fcl_fp alpha) {
    fcl_fp r;
    fcl_fp value = 0.f;
    fcl_fp y = alpha;

    if (alpha <= FCL_DIGAMMA_MAX_ARG_FOR_LINEAR_APPROX) {
        return FCL_DIGAMMA_D1 - 1.f / y;
    }
    while (y < FCL_DIGAMMA_MIN_ARG_FOR_POLYNOMIAL_APPROX) {
        value -= 1.f / y;
        y += 1.f;
    }
    value += log(y) - 0.5f / y;
    r = 1.f / (y * y);

    return value - r*(FCL_DIGAMMA_S3 - r*(FCL_DIGAMMA_S4 - r*FCL_DIGAMMA_S5));
}

fcl_fp16 digamma16(fcl_fp16 alpha) {
    fcl_fp16 r;
    fcl_fp16 y = alpha;
    fcl_fp16 value = (fcl_fp16)(0.f);
    while (any(y < FCL_DIGAMMA_MIN_ARG_FOR_POLYNOMIAL_APPROX)) {
        value -= 1.f / y;
        y += 1.f;
    }
    value += log(y) - 0.5f / y;
    r = 1.f / (y * y);

    return select(value - r*(FCL_DIGAMMA_S3 - r*(FCL_DIGAMMA_S4 - r*FCL_DIGAMMA_S5)), 
            FCL_DIGAMMA_D1 - 1.f / alpha, alpha <= FCL_DIGAMMA_MAX_ARG_FOR_LINEAR_APPROX);
}

fcl_fp trigamma(fcl_fp alpha) {
    fcl_fp value = 0.f;
    fcl_fp y;
    fcl_fp z = alpha;

    if (alpha <= FCL_TRIGAMMA_MAX_ARG_FOR_LINEAR_APPROX) {
        return 1.f / (alpha * alpha);
    }

    while (z < FCL_TRIGAMMA_MIN_ARG_FOR_POLYNOMIAL_APPROX) {
        value += 1.f / (z * z);
        z += 1.f;
    }

    y = 1.f / (z * z);

    return value + 0.5f * y + ( 1.f
            + y * ( FCL_TRIGAMMA_B2
                + y * ( FCL_TRIGAMMA_B4
                    + y * ( FCL_TRIGAMMA_B6
                        + y * FCL_TRIGAMMA_B8 )))) / z;
}

fcl_fp16 trigamma16(fcl_fp16 alpha) {
    fcl_fp16 value = 0.f;
    fcl_fp16 y;
    fcl_fp16 z = alpha;

    while (any(z < FCL_TRIGAMMA_MIN_ARG_FOR_POLYNOMIAL_APPROX)) {
        value += 1.f / (z * z);
        z += 1.f;
    }

    y = 1.f / (z * z);

    return select(value + 0.5f * y + ( 1.f
            + y * ( FCL_TRIGAMMA_B2
                + y * ( FCL_TRIGAMMA_B4
                    + y * ( FCL_TRIGAMMA_B6
                        + y * FCL_TRIGAMMA_B8 )))) / z,
                        1.f / (alpha * alpha), alpha <= FCL_TRIGAMMA_MAX_ARG_FOR_LINEAR_APPROX);
}

int4 convert_fp4_to_int4(fcl_fp4 res_fp) {
    int4 res_i;
    res_i.s0 = res_fp.s0;
    res_i.s1 = res_fp.s1;
    res_i.s2 = res_fp.s2;
    res_i.s3 = res_fp.s3;
    return res_i;
}

int16 convert_fp16_to_int16(fcl_fp16 res_fp) {
    int16 res_i;
    res_i.lo.lo = convert_fp4_to_int4(res_fp.lo.lo);
    res_i.lo.hi = convert_fp4_to_int4(res_fp.lo.hi);
    res_i.hi.lo = convert_fp4_to_int4(res_fp.hi.lo);
    res_i.hi.hi = convert_fp4_to_int4(res_fp.hi.hi);
    return res_i;
}

fcl_fp4 convert_int4_to_fp4(int4 res_i) {
    fcl_fp4 res_fp;
    res_fp.s0 = res_i.s0;
    res_fp.s1 = res_i.s1;
    res_fp.s2 = res_i.s2;
    res_fp.s3 = res_i.s3;
    return res_fp;
}

fcl_fp16 convert_int16_to_fp16(int16 res_i) {
    fcl_fp16 res_fp;
    res_fp.lo.lo = convert_int4_to_fp4(res_i.lo.lo);
    res_fp.lo.hi = convert_int4_to_fp4(res_i.lo.hi);
    res_fp.hi.lo = convert_int4_to_fp4(res_i.hi.lo);
    res_fp.hi.hi = convert_int4_to_fp4(res_i.hi.hi);
    return res_fp;
}

int4 int_array_elements4(__global int *array, int4 idx) {
    int4 res_i;
    res_i.s0 = array[idx.s0];
    res_i.s1 = array[idx.s1];
    res_i.s2 = array[idx.s2];
    res_i.s3 = array[idx.s3];
    return res_i;
}

int16 int_array_elements16(__global int *array, int16 idx) {
    int16 res_i;
    res_i.lo.lo = int_array_elements4(array, idx.lo.lo);
    res_i.lo.hi = int_array_elements4(array, idx.lo.hi);
    res_i.hi.lo = int_array_elements4(array, idx.hi.lo);
    res_i.hi.hi = int_array_elements4(array, idx.hi.hi);
    return res_i;
}

__kernel void fcl_compute_cox_proportional_hazards_work(
        __global fcl_fp *data,
        __global fcl_fp *rank2weight,
        __global fcl_fp *rank2y,
        __global int *independent_cols,
        int slices,
        __global int *id2col,
        __global fcl_fp *theta,
        __global int *rank2row,
        __global int *stratum_offsets,
        __global int *stratum_sizes,
        int cox_work_dim,
        __global fcl_fp *cox_work,
        __global fcl_fp *cox_work_sum
) {
    int s = get_global_id(0);
    int j = get_global_id(1);

    int rows = stratum_sizes[s];
    int offset = stratum_offsets[s];

    int idx = id2col[j];
    int i;

    __global fcl_fp *col_work = &cox_work[cox_work_dim*j];
    __global fcl_fp *col_work_sum = &cox_work_sum[cox_work_dim*j];

    __global fcl_fp *col = &data[independent_cols[idx]*16*slices];

    for (i=0; i<rows; i++) {
        col_work[offset+i] = rank2weight[offset+i] * col[rank2row[offset+i]] * theta[offset+i];
    }

    if (rows > 0) {
        col_work_sum[offset+0] = col_work[offset+0];
    }

    for (i=1; i<rows; i++) {
        col_work_sum[offset+i] = (col_work[offset+i] + col_work_sum[offset+i-1]);
    }

    // Fix up ties
    for (i=rows-1; i>=1; i--) {
        if (rank2y[offset+i] == rank2y[offset+i-1]) {
            col_work_sum[offset+i-1] = col_work_sum[offset+i];
        }
    }
}

__kernel void fcl_compute_jacobian_cox_prortional_hazards(
		__global fcl_fp *rank2weight,
        __global fcl_fp *theta,
        __global fcl_fp *cumulative_theta,
        __global int *censoring,
        int rows,
        __global fcl_fp *cox_work,
        __global fcl_fp *cox_work_sum,
        __global fcl_fp *jacobian_out
) {
    int k = get_global_id(0);
    fcl_fp result = 0.f;
    __global fcl_fp *col_work = &cox_work[rows*k];
    __global fcl_fp *col_work_sum = &cox_work_sum[rows*k];
    int i;

    for (i=0; i<rows; i++) {
        result += (censoring[i] == 0) ? 
            col_work[i] / theta[i] - rank2weight[i] * col_work_sum[i] / cumulative_theta[i] :
            0.f;
    }

    jacobian_out[k] = result;
}

__kernel void fcl_compute_hessian_cox_proportional_hazards(
		__global fcl_fp *rank2weight,
		__global fcl_fp *rank2y,
        __global fcl_fp *theta,
        __global fcl_fp *cumulative_theta,
        __global int *stratum_offsets,
        __global int *stratum_sizes,
        int stratum_count,
        int cox_work_dim,
        __global fcl_fp *cox_work,
        __global fcl_fp *cox_work_sum,
        __global fcl_fp *hessian_out,
        int hessian_cols_count
) {
    int k1 = get_global_id(0);
    int k2 = get_global_id(1);

    __global fcl_fp *col_work1 = &cox_work[cox_work_dim*k1];
    __global fcl_fp *col_work2 = &cox_work[cox_work_dim*k2];

    __global fcl_fp *col_work_sum1 = &cox_work_sum[cox_work_dim*k1];
    __global fcl_fp *col_work_sum2 = &cox_work_sum[cox_work_dim*k2];

    fcl_fp result = 0.f;

    int i, s;

    if (k1 > k2)
        return;

    for (s=0; s<stratum_count; s++) {
        int rows = stratum_sizes[s];
        int offset = stratum_offsets[s];

        fcl_fp cumulative_cross_product = 0.f;
        fcl_fp cumulative_cohort_cross_product = 0.f;

        for (i=0; i<rows; i++) {
            cumulative_cross_product += col_work1[offset+i] * col_work2[offset+i] / (rank2weight[offset+i] * theta[offset+i]);

            result += (rank2weight[offset+i] * 
                    col_work_sum1[offset+i] / cumulative_theta[offset+i] * 
                    col_work_sum2[offset+i] / cumulative_theta[offset+i]);
        }

        // Fix up ties 
        for (i=rows-1; i>=1; i--) {
            result -= rank2weight[offset+i] * cumulative_cross_product / cumulative_theta[offset+i];

            cumulative_cohort_cross_product += col_work1[offset+i] * col_work2[offset+i] / (rank2weight[offset+i] * theta[offset+i]);

            if (rank2y[offset+i] != rank2y[offset+i-1]) {
                cumulative_cross_product -= cumulative_cohort_cross_product;
                cumulative_cohort_cross_product = 0.f;
            }
        }

        result -= rank2weight[offset+0] * cumulative_cross_product / cumulative_theta[offset+0];
    }

	hessian_out[k1*hessian_cols_count+k2] = result;
	hessian_out[k2*hessian_cols_count+k1] = result;
}
__kernel void fcl_compute_fixed_effects_nuisance_jacobian_cox_proportional_hazards(
        __global fcl_fp *global_rank2weight,
        __global fcl_fp *global_rank2y,
        __global fcl_fp *cumulative_theta,
        __global int *stratum_offsets,
        __global int *stratum_rows,
        int stratum_count,

        __global fcl_fp *panel_rank2weight,
        __global fcl_fp *panel_rank2y,
        __global fcl_fp *panel_cumulative_theta,
        __global int *panel_stratum_offsets,
        __global int *panel_stratum_rows,

        __global int *censoring,
        int hessian_cols_count,
        __global fcl_fp *jacobian_out
) {
    int j = get_global_id(0);
    int i, s;

    int core_cols_count = hessian_cols_count - get_global_size(0);

    fcl_fp result = 0.f;

    for (s=0; s<stratum_count; s++) {
        int global_offset = stratum_offsets[s];
        int global_rows = stratum_rows[s];

        int panel_offset = panel_stratum_offsets[j*stratum_count+s];
        int panel_rows = panel_stratum_rows[j*stratum_count+s];

        int panel_i = 0;
        int row_i = 0, panel_row_i = 0;

        for (i=0; i<panel_rows; i++) {
            result += (censoring[panel_offset+i] == 0) ?
                panel_rank2weight[panel_offset+i] :
                0.f;
        }

        for (i=0; i<global_rows; i++) {
            while (panel_i < panel_rows - 1 && 
                panel_rank2y[panel_offset+panel_i+1] >= global_rank2y[global_offset+i]) {
                panel_i++;
            }
            if (panel_rank2y[panel_offset+panel_i] < global_rank2y[global_offset+i])
                continue;

            fcl_fp p = (panel_cumulative_theta[panel_offset+panel_i] / cumulative_theta[global_offset+i]);
            result -= global_rank2weight[global_offset+i] * p;
        }
    }

    jacobian_out[core_cols_count+j] = result;
}

__kernel void fcl_compute_fixed_effects_cross_hessian_cox_proportional_hazards(
        __global fcl_fp *global_rank2weight,
        __global fcl_fp *global_rank2y,
        __global fcl_fp *cumulative_theta,
        __global int *stratum_offsets,
        __global int *stratum_rows,
        int stratum_count,

        __global fcl_fp *panel_rank2y,
        __global fcl_fp *panel_cumulative_theta,
        __global int *panel_stratum_offsets,
        __global int *panel_stratum_rows,

        int cox_work_dim,
        __global fcl_fp *cox_work_sum,
        __global fcl_fp *fe_cox_work_sum,
        __global fcl_fp *hessian_out
) {
    int k = get_global_id(0);
    int j = get_global_id(1);
    int i, s;

    __global fcl_fp *col_work_sum = &cox_work_sum[cox_work_dim*k];
    __global fcl_fp *fe_col_work_sum = &fe_cox_work_sum[cox_work_dim*k];

    int core_cols_count = get_global_size(0);
    int hessian_cols_count = core_cols_count + get_global_size(1);

    fcl_fp result = 0.f;

    for (s=0; s<stratum_count; s++) {
        int global_offset = stratum_offsets[s];
        int global_rows = stratum_rows[s];

        int panel_offset = panel_stratum_offsets[j*stratum_count+s];
        int panel_rows = panel_stratum_rows[j*stratum_count+s];

        int panel_i = 0;

        for (i=0; i<global_rows; i++) {
            while (panel_i < panel_rows - 1 && 
                panel_rank2y[panel_offset+panel_i+1] >= global_rank2y[global_offset+i]) {
                panel_i++;
            }
            if (panel_rank2y[panel_offset+panel_i] < global_rank2y[global_offset+i])
                continue;

            result += -(global_rank2weight[global_offset+i] * 
                    col_work_sum[global_offset+i] / cumulative_theta[global_offset+i] * 
                    panel_cumulative_theta[panel_offset+panel_i] / cumulative_theta[global_offset+i]);

            result += (global_rank2weight[global_offset+i] * 
                    fe_col_work_sum[panel_offset+panel_i] / cumulative_theta[global_offset+i]);
        }
    }

    hessian_out[hessian_cols_count*(core_cols_count+j)+k] = -result;
    hessian_out[hessian_cols_count*k+(core_cols_count+j)] = -result;
}

__kernel void fcl_compute_fixed_effects_offdiagonal_hessian_cox_proportional_hazards(
        __global fcl_fp *global_rank2weight,
        __global fcl_fp *global_rank2y,
        __global fcl_fp *cumulative_theta,
        __global int *stratum_offsets,
        __global int *stratum_rows,
        int stratum_count,

        __global fcl_fp *panel_rank2y,
        __global fcl_fp *panel_cumulative_theta,
        __global int *panel_stratum_offsets,
        __global int *panel_stratum_rows,

        int hessian_cols_count,
        __global fcl_fp *hessian_out
) {
    int j1 = get_global_id(0);
    int j2 = get_global_id(1);
    int i, s;

    int core_cols_count = hessian_cols_count - get_global_size(0);

    fcl_fp result = 0.f;

    if (j1 >= j2)
        return;

    for (s=0; s<stratum_count; s++) {
        int global_offset = stratum_offsets[s];
        int global_rows = stratum_rows[s];

        int panel1_offset = panel_stratum_offsets[j1*stratum_count+s];
        int panel1_rows = panel_stratum_rows[j1*stratum_count+s];

        int panel2_offset = panel_stratum_offsets[j2*stratum_count+s];
        int panel2_rows = panel_stratum_rows[j2*stratum_count+s];

        int panel1_i = 0;
        int panel2_i = 0;

        for (i=0; i<global_rows; i++) {
            while (panel1_i < panel1_rows - 1 && 
                panel_rank2y[panel1_offset+panel1_i+1] >= global_rank2y[global_offset+i]) {
                panel1_i++;
            }
            if (panel_rank2y[panel1_offset+panel1_i] < global_rank2y[global_offset+i])
                continue;

            while (panel2_i < panel2_rows - 1 && 
                panel_rank2y[panel2_offset+panel2_i+1] >= global_rank2y[global_offset+i]) {
                panel2_i++;
            }
            if (panel_rank2y[panel2_offset+panel2_i] < global_rank2y[global_offset+i])
                continue;

            result += (global_rank2weight[global_offset+i] * 
                    panel_cumulative_theta[panel1_offset+panel1_i] / cumulative_theta[global_offset+i] *
                    panel_cumulative_theta[panel2_offset+panel2_i] / cumulative_theta[global_offset+i]);
        }
    }

    hessian_out[hessian_cols_count*(core_cols_count+j1)+core_cols_count+j2] = result;
    hessian_out[hessian_cols_count*(core_cols_count+j2)+core_cols_count+j1] = result;
}

__kernel void fcl_compute_fixed_effects_diagonal_hessian_cox_proportional_hazards(
        __global fcl_fp *global_rank2weight,
        __global fcl_fp *global_rank2y,
        __global fcl_fp *cumulative_theta,
        __global int *stratum_offsets,
        __global int *stratum_rows,
        int stratum_count,

        __global fcl_fp *panel_rank2y,
        __global fcl_fp *panel_cumulative_theta,
        __global int *panel_stratum_offsets,
        __global int *panel_stratum_rows,

        int hessian_cols_count,
        __global fcl_fp *hessian_out
) {
    int j = get_global_id(0);
    int i, s, panel_i;
    int core_cols_count = hessian_cols_count - get_global_size(0);
    fcl_fp result = 0.f;

    for (s=0; s<stratum_count; s++) {
        int global_offset = stratum_offsets[s];
        int global_rows = stratum_rows[s];

        int panel_offset = panel_stratum_offsets[j*stratum_count+s];
        int panel_rows = panel_stratum_rows[j*stratum_count+s];

        panel_i = 0;

        for (i=0; i<global_rows; i++) {
            while (panel_i < panel_rows - 1 && 
                panel_rank2y[panel_offset+panel_i+1] >= global_rank2y[global_offset+i]) {
                panel_i++;
            }
            if (panel_rank2y[panel_offset+panel_i] < global_rank2y[global_offset+i])
                continue;

            fcl_fp p = (panel_cumulative_theta[panel_offset+panel_i] / cumulative_theta[global_offset+i]);
            result += global_rank2weight[global_offset+i] * p * (1.f - p);
        }
    }

    hessian_out[hessian_cols_count*(core_cols_count+j)+core_cols_count+j] = -result;
}
