#include "config.h"

#include <math.h>
#include <float.h>
#include <string.h>
#include <stdlib.h>
#include <assert.h>
#include <stdarg.h>
#include "kann.h"

int kann_verbose = 3;

/******************************************
 *** @@BASIC: fundamental KANN routines ***
 ******************************************/

static void kad_ext_collate(int n, kad_node_t **a, float **_x, float **_g, float **_c)
{
	int i, j, k, l, n_var;
	float *x, *g, *c;
	n_var = kad_size_var(n, a);
	x = *_x = (float*)realloc(*_x, n_var * sizeof(float));
	g = *_g = (float*)realloc(*_g, n_var * sizeof(float));
	c = *_c = (float*)realloc(*_c, kad_size_const(n, a) * sizeof(float));
	memset(g, 0, n_var * sizeof(float));
	for (i = j = k = 0; i < n; ++i) {
		kad_node_t *v = a[i];
		if (kad_is_var(v)) {
			l = kad_len(v);
			memcpy(&x[j], v->x, l * sizeof(float));
			free(v->x);
			v->x = &x[j];
			v->g = &g[j];
			j += l;
		} else if (kad_is_const(v)) {
			l = kad_len(v);
			memcpy(&c[k], v->x, l * sizeof(float));
			free(v->x);
			v->x = &c[k];
			k += l;
		}
	}
}

static void kad_ext_sync(int n, kad_node_t **a, float *x, float *g, float *c)
{
	int i, j, k;
	for (i = j = k = 0; i < n; ++i) {
		kad_node_t *v = a[i];
		if (kad_is_var(v)) {
			v->x = &x[j];
			v->g = &g[j];
			j += kad_len(v);
		} else if (kad_is_const(v)) {
			v->x = &c[k];
			k += kad_len(v);
		}
	}
}

kann_t *kann_new(kad_node_t *cost, int n_rest, ...)
{
	kann_t *a;
	int i, n_roots = 1 + n_rest, has_pivot = 0, has_recur = 0;
	kad_node_t **roots;
	va_list ap;

	if (cost->n_d != 0) return 0;

	va_start(ap, n_rest);
	roots = (kad_node_t**)malloc((n_roots + 1) * sizeof(kad_node_t*));
	for (i = 0; i < n_rest; ++i)
		roots[i] = va_arg(ap, kad_node_t*);
	roots[i++] = cost;
	va_end(ap);

	cost->ext_flag |= KANN_F_COST;
	a = (kann_t*)calloc(1, sizeof(kann_t));
	a->v = kad_compile_array(&a->n, n_roots, roots);

	for (i = 0; i < a->n; ++i) {
		if (a->v[i]->pre) has_recur = 1;
		if (kad_is_pivot(a->v[i])) has_pivot = 1;
	}
	if (has_recur && !has_pivot) { /* an RNN that doesn't have a pivot; then add a pivot on top of cost and recompile */
		cost->ext_flag &= ~KANN_F_COST;
		roots[n_roots-1] = cost = kad_avg(1, &cost), cost->ext_flag |= KANN_F_COST;
		free(a->v);
		a->v = kad_compile_array(&a->n, n_roots, roots);
	}
	kad_ext_collate(a->n, a->v, &a->x, &a->g, &a->c);
	free(roots);
	return a;
}

kann_t *kann_clone(kann_t *a, int batch_size)
{
	kann_t *b;
	b = (kann_t*)calloc(1, sizeof(kann_t));
	b->n = a->n;
	b->v = kad_clone(a->n, a->v, batch_size);
	kad_ext_collate(b->n, b->v, &b->x, &b->g, &b->c);
	return b;
}

kann_t *kann_unroll_array(kann_t *a, int *len)
{
	kann_t *b;
	b = (kann_t*)calloc(1, sizeof(kann_t));
	b->x = a->x, b->g = a->g, b->c = a->c; /* these arrays are shared */
	b->v = kad_unroll(a->n, a->v, &b->n, len);
	return b;
}

kann_t *kann_unroll(kann_t *a, ...)
{
	kann_t *b;
	va_list ap;
	int i, n_pivots, *len;
	n_pivots = kad_n_pivots(a->n, a->v);
	len = (int*)calloc(n_pivots, sizeof(int));
	va_start(ap, a);
	for (i = 0; i < n_pivots; ++i) len[i] = va_arg(ap, int);
	va_end(ap);
	b = kann_unroll_array(a, len);
	free(len);
	return b;
}

void kann_delete_unrolled(kann_t *a)
{
	if (a && a->mt) kann_mt(a, 0, 0);
	if (a && a->v) kad_delete(a->n, a->v);
	free(a);
}

void kann_delete(kann_t *a)
{
	if (a == 0) return;
	free(a->x); free(a->g); free(a->c);
	kann_delete_unrolled(a);
}

static void kann_switch_core(kann_t *a, int is_train)
{
	int i;
	for (i = 0; i < a->n; ++i)
		if (a->v[i]->op == 12 && a->v[i]->n_child == 2)
			*(int32_t*)a->v[i]->ptr = !!is_train;
}

#define chk_flg(flag, mask) ((mask) == 0 || ((flag) & (mask)))
#define chk_lbl(label, query) ((query) == 0 || (label) == (query))

int kann_find(const kann_t *a, uint32_t ext_flag, int32_t ext_label)
{
	int i, k, r = -1;
	for (i = k = 0; i < a->n; ++i)
		if (chk_flg(a->v[i]->ext_flag, ext_flag) && chk_lbl(a->v[i]->ext_label, ext_label))
			++k, r = i;
	return k == 1? r : k == 0? -1 : -2;
}

int kann_feed_bind(kann_t *a, uint32_t ext_flag, int32_t ext_label, float **x)
{
	int i, k;
	if (x == 0) return 0;
	for (i = k = 0; i < a->n; ++i)
		if (kad_is_feed(a->v[i]) && chk_flg(a->v[i]->ext_flag, ext_flag) && chk_lbl(a->v[i]->ext_label, ext_label))
			a->v[i]->x = x[k++];
	return k;
}

int kann_feed_dim(const kann_t *a, uint32_t ext_flag, int32_t ext_label)
{
	int i, k, n = 0;
	for (i = k = 0; i < a->n; ++i)
		if (kad_is_feed(a->v[i]) && chk_flg(a->v[i]->ext_flag, ext_flag) && chk_lbl(a->v[i]->ext_label, ext_label))
			++k, n = a->v[i]->n_d > 1? kad_len(a->v[i]) / a->v[i]->d[0] : a->v[i]->n_d == 1? a->v[i]->d[0] : 1;
	return k == 1? n : k == 0? -1 : -2;
}

static float kann_cost_core(kann_t *a, int cost_label, int cal_grad)
{
	int i_cost;
	float cost;
	i_cost = kann_find(a, KANN_F_COST, cost_label);
	assert(i_cost >= 0);
	cost = *kad_eval_at(a->n, a->v, i_cost);
	if (cal_grad) kad_grad(a->n, a->v, i_cost);
	return cost;
}

int kann_eval(kann_t *a, uint32_t ext_flag, int ext_label)
{
	int i, k;
	for (i = k = 0; i < a->n; ++i)
		if (chk_flg(a->v[i]->ext_flag, ext_flag) && chk_lbl(a->v[i]->ext_label, ext_label))
			++k, a->v[i]->tmp = 1;
	kad_eval_marked(a->n, a->v);
	return k;
}

void kann_rnn_start(kann_t *a)
{
	int i;
	kann_set_batch_size(a, 1);
	for (i = 0; i < a->n; ++i) {
		kad_node_t *p = a->v[i];
		if (p->pre) { /* NB: BE CAREFUL of the interaction between kann_rnn_start() and kann_set_batch_size() */
			kad_node_t *q = p->pre;
			if (q->x) memcpy(p->x, q->x, kad_len(p) * sizeof(float));
			else memset(p->x, 0, kad_len(p) * sizeof(float));
			if (q->n_child > 0) free(q->x);
			q->x = p->x;
		}
	}
}

void kann_rnn_end(kann_t *a)
{
	int i;
	kad_ext_sync(a->n, a->v, a->x, a->g, a->c);
	for (i = 0; i < a->n; ++i)
		if (a->v[i]->pre && a->v[i]->pre->n_child > 0)
			a->v[i]->pre->x = (float*)calloc(kad_len(a->v[i]->pre), sizeof(float));
}

static int kann_class_error_core(const kann_t *ann, int *base)
{
	int i, j, k, m, n, off, n_err = 0;
	for (i = 0, *base = 0; i < ann->n; ++i) {
		kad_node_t *p = ann->v[i];
		if (((p->op == 13 && (p->n_child == 2 || p->n_child == 3)) || (p->op == 22 && p->n_child == 2)) && p->n_d == 0) { /* ce_bin or ce_multi */
			kad_node_t *x = p->child[0], *t = p->child[1];
			n = t->d[t->n_d - 1], m = kad_len(t) / n;
			for (j = off = 0; j < m; ++j, off += n) {
				float t_sum = 0.0f, t_min = 1.0f, t_max = 0.0f, x_max = 0.0f, x_min = 1.0f;
				int x_max_k = -1, t_max_k = -1;
				for (k = 0; k < n; ++k) {
					float xk = x->x[off+k], tk = t->x[off+k];
					t_sum += tk;
					t_min = t_min < tk? t_min : tk;
					x_min = x_min < xk? x_min : xk;
					if (t_max < tk) t_max = tk, t_max_k = k;
					if (x_max < xk) x_max = xk, x_max_k = k;
				}
				if (t_sum - 1.0f == 0 && t_min >= 0.0f && x_min >= 0.0f && x_max <= 1.0f) {
					++(*base);
					n_err += (x_max_k != t_max_k);
				}
			}
		}
	}
	return n_err;
}

/*************************
 * @@MT: multi-threading *
 *************************/

#ifdef HAVE_PTHREAD
#include <pthread.h>

struct mtaux_t;

typedef struct { /* per-worker data */
	kann_t *a;
	float cost;
	int action;
	pthread_t tid;
	struct mtaux_t *g;
} mtaux1_t;

typedef struct mtaux_t { /* cross-worker data */
	int n_threads, max_batch_size;
	int cal_grad, cost_label, eval_out;
	volatile int n_idle; /* we will be busy waiting on this, so volatile necessary */
	pthread_mutex_t mtx;
	pthread_cond_t cv;
	mtaux1_t *mt;
} mtaux_t;

static void *mt_worker(void *data) /* pthread worker */
{
	mtaux1_t *mt1 = (mtaux1_t*)data;
	mtaux_t *mt = mt1->g;
	for (;;) {
		int action;
		pthread_mutex_lock(&mt->mtx);
		mt1->action = 0;
		++mt->n_idle;
		while (mt1->action == 0)
			pthread_cond_wait(&mt->cv, &mt->mtx);
		action = mt1->action;
		pthread_mutex_unlock(&mt->mtx);
		if (action == -1) break;

		if (mt->eval_out) kann_eval(mt1->a, KANN_F_OUT, 0);
		else mt1->cost = kann_cost_core(mt1->a, mt->cost_label, mt->cal_grad);
	}
	pthread_exit(0);
}

static void mt_destroy(mtaux_t *mt) /* de-allocate an entire mtaux_t struct */
{
	int i;
	pthread_mutex_lock(&mt->mtx);
	mt->n_idle = 0;
	for (i = 1; i < mt->n_threads; ++i) mt->mt[i].action = -1;
	pthread_cond_broadcast(&mt->cv);
	pthread_mutex_unlock(&mt->mtx);
	for (i = 1; i < mt->n_threads; ++i) pthread_join(mt->mt[i].tid, 0);
	for (i = 0; i < mt->n_threads; ++i) kann_delete(mt->mt[i].a);
	free(mt->mt);
	pthread_cond_destroy(&mt->cv);
	pthread_mutex_destroy(&mt->mtx);
	free(mt);
}

void kann_mt(kann_t *ann, int n_threads, int max_batch_size)
{
	mtaux_t *mt;
	int i, k;

	if (n_threads <= 1) {
		if (ann->mt) mt_destroy((mtaux_t*)ann->mt);
		ann->mt = 0;
		return;
	}
	if (n_threads > max_batch_size) n_threads = max_batch_size;
	if (n_threads <= 1) return;

	mt = (mtaux_t*)calloc(1, sizeof(mtaux_t));
	mt->n_threads = n_threads, mt->max_batch_size = max_batch_size;
	pthread_mutex_init(&mt->mtx, 0);
	pthread_cond_init(&mt->cv, 0);
	mt->mt = (mtaux1_t*)calloc(n_threads, sizeof(mtaux1_t));
	for (i = k = 0; i < n_threads; ++i) {
		int size = (max_batch_size - k) / (n_threads - i);
		mt->mt[i].a = kann_clone(ann, size);
		mt->mt[i].g = mt;
		k += size;
	}
	for (i = 1; i < n_threads; ++i)
		pthread_create(&mt->mt[i].tid, 0, mt_worker, &mt->mt[i]);
	while (mt->n_idle < n_threads - 1); /* busy waiting until all threads in sync */
	ann->mt = mt;
}

static void mt_kickoff(kann_t *a, int cost_label, int cal_grad, int eval_out)
{
	mtaux_t *mt = (mtaux_t*)a->mt;
	int i, j, k, B, n_var;

	B = kad_sync_dim(a->n, a->v, -1); /* get the current batch size */
	assert(B <= mt->max_batch_size); /* TODO: can be relaxed */
	n_var = kann_size_var(a);

	pthread_mutex_lock(&mt->mtx);
	mt->cost_label = cost_label, mt->cal_grad = cal_grad, mt->eval_out = eval_out;
	for (i = k = 0; i < mt->n_threads; ++i) {
		int size = (B - k) / (mt->n_threads - i);
		for (j = 0; j < a->n; ++j)
			if (kad_is_feed(a->v[j]))
				mt->mt[i].a->v[j]->x = &a->v[j]->x[k * kad_len(a->v[j]) / a->v[j]->d[0]];
		kad_sync_dim(mt->mt[i].a->n, mt->mt[i].a->v, size); /* TODO: we can point ->x to internal nodes, too */
		k += size;
		memcpy(mt->mt[i].a->x, a->x, n_var * sizeof(float));
		mt->mt[i].action = 1;
	}
	mt->n_idle = 0;
	pthread_cond_broadcast(&mt->cv);
	pthread_mutex_unlock(&mt->mtx);
}

float kann_cost(kann_t *a, int cost_label, int cal_grad)
{
	mtaux_t *mt = (mtaux_t*)a->mt;
	int i, j, B, k, n_var;
	float cost;

	if (mt == 0) return kann_cost_core(a, cost_label, cal_grad);
	B = kad_sync_dim(a->n, a->v, -1); /* get the current batch size */
	n_var = kann_size_var(a);

	mt_kickoff(a, cost_label, cal_grad, 0);
	mt->mt[0].cost = kann_cost_core(mt->mt[0].a, cost_label, cal_grad);
	while (mt->n_idle < mt->n_threads - 1); /* busy waiting until all threads in sync */

	memset(a->g, 0, n_var * sizeof(float)); /* TODO: check if this is necessary when cal_grad is false */
	for (i = k = 0, cost = 0.0f; i < mt->n_threads; ++i) {
		int size = (B - k) / (mt->n_threads - i);
		cost += mt->mt[i].cost * size / B;
		kad_saxpy(n_var, (float)size / B, mt->mt[i].a->g, a->g);
		k += size;
	}
	for (j = 0; j < a->n; ++j) { /* copy values back at recurrent nodes (needed by textgen; TODO: temporary solution) */
		kad_node_t *p = a->v[j];
		if (p->pre && p->n_d >= 2 && p->d[0] == B) {
			for (i = k = 0; i < mt->n_threads; ++i) {
				kad_node_t *q = mt->mt[i].a->v[j];
				memcpy(&p->x[k], q->x, kad_len(q) * sizeof(float));
				k += kad_len(q);
			}
		}
	}
	return cost;
}

int kann_eval_out(kann_t *a)
{
	mtaux_t *mt = (mtaux_t*)a->mt;
	int j, B, n_eval;
	if (mt == 0) return kann_eval(a, KANN_F_OUT, 0);
	B = kad_sync_dim(a->n, a->v, -1); /* get the current batch size */
	mt_kickoff(a, 0, 0, 1);
	n_eval = kann_eval(mt->mt[0].a, KANN_F_OUT, 0);
	while (mt->n_idle < mt->n_threads - 1); /* busy waiting until all threads in sync */
	for (j = 0; j < a->n; ++j) { /* copy output values back */
		kad_node_t *p = a->v[j];
		if (p->ext_flag & KANN_F_OUT) {
			int i, t, k, d0 = p->d[0] / B, d1 = 1; /* for RNN, p->d[0] may equal unroll_len * batch_size */
			assert(p->d[0] % B == 0);
			for (i = 1; i < p->n_d; ++i) d1 *= p->d[i];
			for (i = 0; i < d0; ++i) {
				for (t = k = 0; t < mt->n_threads; ++t) { /* similar to the forward pass of kad_op_concat() */
					kad_node_t *q = mt->mt[t].a->v[j];
					int size = q->d[0] / d0;
					memcpy(&p->x[(i * B + k) * d1], &q->x[i * size * d1], size * d1 * sizeof(float));
					k += size;
				}
			}
		}
	}
	return n_eval;
}

int kann_class_error(const kann_t *ann, int *base)
{
	mtaux_t *mt = (mtaux_t*)ann->mt;
	int i, n_err = 0, b = 0;
	if (mt == 0) return kann_class_error_core(ann, base);
	for (i = 0; i < mt->n_threads; ++i) {
		n_err += kann_class_error_core(mt->mt[i].a, &b);
		*base += b;
	}
	return n_err;
}

void kann_switch(kann_t *ann, int is_train)
{
	mtaux_t *mt = (mtaux_t*)ann->mt;
	int i;
	if (mt == 0) {
		kann_switch_core(ann, is_train);
		return;
	}
	for (i = 0; i < mt->n_threads; ++i)
		kann_switch_core(mt->mt[i].a, is_train);
}
#else
void kann_mt(kann_t *ann, int n_threads, int max_batch_size) {}
float kann_cost(kann_t *a, int cost_label, int cal_grad) { return kann_cost_core(a, cost_label, cal_grad); }
int kann_eval_out(kann_t *a) { return kann_eval(a, KANN_F_OUT, 0); }
int kann_class_error(const kann_t *a, int *base) { return kann_class_error_core(a, base); }
void kann_switch(kann_t *ann, int is_train) { return kann_switch_core(ann, is_train); }
#endif

/***********************
 *** @@IO: model I/O ***
 ***********************/

#define KANN_MAGIC "KAN\1"

void kann_save_fp(FILE *fp, kann_t *ann)
{
	kann_set_batch_size(ann, 1);
	fwrite(KANN_MAGIC, 1, 4, fp);
	kad_save(fp, ann->n, ann->v);
	fwrite(ann->x, sizeof(float), kann_size_var(ann), fp);
	fwrite(ann->c, sizeof(float), kann_size_const(ann), fp);
}

void kann_save(const char *fn, kann_t *ann)
{
	FILE *fp;
	fp = fn && strcmp(fn, "-")? fopen(fn, "wb") : stdout;
	kann_save_fp(fp, ann);
	fclose(fp);
}

kann_t *kann_load_fp(FILE *fp)
{
	char magic[4];
	kann_t *ann;
	int n_var, n_const;

	(void) !fread(magic, 1, 4, fp);
	if (strncmp(magic, KANN_MAGIC, 4) != 0) {
		fclose(fp);
		return 0;
	}
	ann = (kann_t*)calloc(1, sizeof(kann_t));
	ann->v = kad_load(fp, &ann->n);
	n_var = kad_size_var(ann->n, ann->v);
	n_const = kad_size_const(ann->n, ann->v);
	ann->x = (float*)malloc(n_var * sizeof(float));
	ann->g = (float*)calloc(n_var, sizeof(float));
	ann->c = (float*)malloc(n_const * sizeof(float));
	(void) !fread(ann->x, sizeof(float), n_var, fp);
	(void) !fread(ann->c, sizeof(float), n_const, fp);
	kad_ext_sync(ann->n, ann->v, ann->x, ann->g, ann->c);
	return ann;
}

kann_t *kann_load(const char *fn)
{
	FILE *fp;
	kann_t *ann;
	fp = fn && strcmp(fn, "-")? fopen(fn, "rb") : stdin;
	ann = kann_load_fp(fp);
	fclose(fp);
	return ann;
}

/**********************************************
 *** @@LAYER: layers and model generation ***
 **********************************************/

/********** General but more complex APIs **********/

kad_node_t *kann_new_leaf_array(int *offset, kad_node_p *par, uint8_t flag, float x0_01, int n_d, int32_t d[KAD_MAX_DIM])
{
	int i, len, off = offset && par? *offset : -1;
	kad_node_t *p;

	if (off >= 0 && par[off]) return par[(*offset)++];
	p = (kad_node_t*)calloc(1, sizeof(kad_node_t));
	p->n_d = n_d, p->flag = flag;
	memcpy(p->d, d, n_d * sizeof(int32_t));
	len = kad_len(p);
	p->x = (float*)calloc(len, sizeof(float));
	if (p->n_d <= 1) {
		for (i = 0; i < len; ++i)
			p->x[i] = x0_01;
	} else {
		double sdev_inv;
		sdev_inv = 1.0 / sqrt((double)len / p->d[0]);
		for (i = 0; i < len; ++i)
			p->x[i] = (float)(kad_drand_normal(0) * sdev_inv);
	}
	if (off >= 0) par[off] = p, ++(*offset);
	return p;
}

kad_node_t *kann_new_leaf2(int *offset, kad_node_p *par, uint8_t flag, float x0_01, int n_d, ...)
{
	int32_t i, d[KAD_MAX_DIM];
	va_list ap;
	va_start(ap, n_d); for (i = 0; i < n_d; ++i) d[i] = va_arg(ap, int); va_end(ap);
	return kann_new_leaf_array(offset, par, flag, x0_01, n_d, d);
}

kad_node_t *kann_layer_dense2(int *offset, kad_node_p *par, kad_node_t *in, int n1)
{
	int n0;
	kad_node_t *w, *b;
	n0 = in->n_d >= 2? kad_len(in) / in->d[0] : kad_len(in);
	w = kann_new_leaf2(offset, par, KAD_VAR, 0.0f, 2, n1, n0);
	b = kann_new_leaf2(offset, par, KAD_VAR, 0.0f, 1, n1);
	return kad_add(kad_cmul(in, w), b);
}

kad_node_t *kann_layer_dropout2(int *offset, kad_node_p *par, kad_node_t *t, float r)
{
	kad_node_t *x[2], *cr;
	cr = kann_new_leaf2(offset, par, KAD_CONST, r, 0);
	x[0] = t, x[1] = kad_dropout(t, cr);
	return kad_switch(2, x);
}

kad_node_t *kann_layer_layernorm2(int *offset, kad_node_t **par, kad_node_t *in)
{
	int n0;
	kad_node_t *alpha, *beta;
	n0 = in->n_d >= 2? kad_len(in) / in->d[0] : kad_len(in);
	alpha = kann_new_leaf2(offset, par, KAD_VAR, 1.0f, 1, n0);
	beta  = kann_new_leaf2(offset, par, KAD_VAR, 0.0f, 1, n0);
	return kad_add(kad_mul(kad_stdnorm(in), alpha), beta);
}

static inline kad_node_t *cmul_norm2(int *offset, kad_node_t **par, kad_node_t *x, kad_node_t *w, int use_norm)
{
	return use_norm? kann_layer_layernorm2(offset, par, kad_cmul(x, w)) : kad_cmul(x, w);
}

kad_node_t *kann_layer_rnn2(int *offset, kad_node_t **par, kad_node_t *in, kad_node_t *h0, int rnn_flag)
{
	int n0, n1 = h0->d[h0->n_d-1], use_norm = !!(rnn_flag & KANN_RNN_NORM);
	kad_node_t *t, *w, *u, *b, *out;

	u = kann_new_leaf2(offset, par, KAD_VAR, 0.0f, 2, n1, n1);
	b = kann_new_leaf2(offset, par, KAD_VAR, 0.0f, 1, n1);
	t = cmul_norm2(offset, par, h0, u, use_norm);
	if (in) {
		n0 = in->n_d >= 2? kad_len(in) / in->d[0] : kad_len(in);
		w = kann_new_leaf2(offset, par, KAD_VAR, 0.0f, 2, n1, n0);
		t = kad_add(cmul_norm2(offset, par, in, w, use_norm), t);
	}
	out = kad_tanh(kad_add(t, b));
	out->pre = h0;
	return out;
}

kad_node_t *kann_layer_gru2(int *offset, kad_node_t **par, kad_node_t *in, kad_node_t *h0, int rnn_flag)
{
	int n0 = 0, n1 = h0->d[h0->n_d-1], use_norm = !!(rnn_flag & KANN_RNN_NORM);
	kad_node_t *t, *r, *z, *w, *u, *b, *s, *out;

	if (in) n0 = in->n_d >= 2? kad_len(in) / in->d[0] : kad_len(in);
	/* z = sigm(x_t * W_z + h_{t-1} * U_z + b_z) */
	u = kann_new_leaf2(offset, par, KAD_VAR, 0.0f, 2, n1, n1);
	b = kann_new_leaf2(offset, par, KAD_VAR, 0.0f, 1, n1);
	t = cmul_norm2(offset, par, h0, u, use_norm);
	if (in) {
		w = kann_new_leaf2(offset, par, KAD_VAR, 0.0f, 2, n1, n0);
		t = kad_add(cmul_norm2(offset, par, in, w, use_norm), t);
	}
	z = kad_sigm(kad_add(t, b));
	/* r = sigm(x_t * W_r + h_{t-1} * U_r + b_r) */
	u = kann_new_leaf2(offset, par, KAD_VAR, 0.0f, 2, n1, n1);
	b = kann_new_leaf2(offset, par, KAD_VAR, 0.0f, 1, n1);
	t = cmul_norm2(offset, par, h0, u, use_norm);
	if (in) {
		w = kann_new_leaf2(offset, par, KAD_VAR, 0.0f, 2, n1, n0);
		t = kad_add(cmul_norm2(offset, par, in, w, use_norm), t);
	}
	r = kad_sigm(kad_add(t, b));
	/* s = tanh(x_t * W_s + (h_{t-1} # r) * U_s + b_s) */
	u = kann_new_leaf2(offset, par, KAD_VAR, 0.0f, 2, n1, n1);
	b = kann_new_leaf2(offset, par, KAD_VAR, 0.0f, 1, n1);
	t = cmul_norm2(offset, par, kad_mul(r, h0), u, use_norm);
	if (in) {
		w = kann_new_leaf2(offset, par, KAD_VAR, 0.0f, 2, n1, n0);
		t = kad_add(cmul_norm2(offset, par, in, w, use_norm), t);
	}
	s = kad_tanh(kad_add(t, b));
	/* h_t = z # h_{t-1} + (1 - z) # s */
	out = kad_add(kad_mul(kad_1minus(z), s), kad_mul(z, h0));
	out->pre = h0;
	return out;
}

/********** APIs without offset & par **********/

kad_node_t *kann_new_leaf(uint8_t flag, float x0_01, int n_d, ...)
{
	int32_t i, d[KAD_MAX_DIM];
	va_list ap;
	va_start(ap, n_d); for (i = 0; i < n_d; ++i) d[i] = va_arg(ap, int); va_end(ap);
	return kann_new_leaf_array(0, 0, flag, x0_01, n_d, d);
}

kad_node_t *kann_new_scalar(uint8_t flag, float x) { return kann_new_leaf(flag, x, 0); }
kad_node_t *kann_new_weight(int n_row, int n_col) { return kann_new_leaf(KAD_VAR, 0.0f, 2, n_row, n_col); }
kad_node_t *kann_new_vec(int n, float x) { return kann_new_leaf(KAD_VAR, x, 1, n); }
kad_node_t *kann_new_bias(int n) { return kann_new_vec(n, 0.0f); }
kad_node_t *kann_new_weight_conv2d(int n_out, int n_in, int k_row, int k_col) { return kann_new_leaf(KAD_VAR, 0.0f, 4, n_out, n_in, k_row, k_col); }
kad_node_t *kann_new_weight_conv1d(int n_out, int n_in, int kernel_len) { return kann_new_leaf(KAD_VAR, 0.0f, 3, n_out, n_in, kernel_len); }

kad_node_t *kann_layer_input(int n1)
{
	kad_node_t *t;
	t = kad_feed(2, 1, n1);
	t->ext_flag |= KANN_F_IN;
	return t;
}

kad_node_t *kann_layer_dense(kad_node_t *in, int n1) { return kann_layer_dense2(0, 0, in, n1); }
kad_node_t *kann_layer_dropout(kad_node_t *t, float r) { return kann_layer_dropout2(0, 0, t, r); }
kad_node_t *kann_layer_layernorm(kad_node_t *in) { return kann_layer_layernorm2(0, 0, in); }

kad_node_t *kann_layer_rnn(kad_node_t *in, int n1, int rnn_flag)
{
	kad_node_t *h0;
	h0 = (rnn_flag & KANN_RNN_VAR_H0)? kad_var(0, 0, 2, 1, n1) : kad_const(0, 2, 1, n1);
	h0->x = (float*)calloc(n1, sizeof(float));
	return kann_layer_rnn2(0, 0, in, h0, rnn_flag);
}

kad_node_t *kann_layer_gru(kad_node_t *in, int n1, int rnn_flag)
{
	kad_node_t *h0;
	h0 = (rnn_flag & KANN_RNN_VAR_H0)? kad_var(0, 0, 2, 1, n1) : kad_const(0, 2, 1, n1);
	h0->x = (float*)calloc(n1, sizeof(float));
	return kann_layer_gru2(0, 0, in, h0, rnn_flag);
}

static kad_node_t *kann_cmul_norm(kad_node_t *x, kad_node_t *w)
{
	return kann_layer_layernorm(kad_cmul(x, w));
}

kad_node_t *kann_layer_lstm(kad_node_t *in, int n1, int rnn_flag)
{
	int n0;
	kad_node_t *i, *f, *o, *g, *w, *u, *b, *h0, *c0, *c, *out;
	kad_node_t *(*cmul)(kad_node_t*, kad_node_t*) = (rnn_flag & KANN_RNN_NORM)? kann_cmul_norm : kad_cmul;

	n0 = in->n_d >= 2? kad_len(in) / in->d[0] : kad_len(in);
	h0 = (rnn_flag & KANN_RNN_VAR_H0)? kad_var(0, 0, 2, 1, n1) : kad_const(0, 2, 1, n1);
	h0->x = (float*)calloc(n1, sizeof(float));
	c0 = (rnn_flag & KANN_RNN_VAR_H0)? kad_var(0, 0, 2, 1, n1) : kad_const(0, 2, 1, n1);
	c0->x = (float*)calloc(n1, sizeof(float));

	/* i = sigm(x_t * W_i + h_{t-1} * U_i + b_i) */
	w = kann_new_weight(n1, n0);
	u = kann_new_weight(n1, n1);
	b = kann_new_bias(n1);
	i = kad_sigm(kad_add(kad_add(cmul(in, w), cmul(h0, u)), b));
	/* f = sigm(x_t * W_f + h_{t-1} * U_f + b_f) */
	w = kann_new_weight(n1, n0);
	u = kann_new_weight(n1, n1);
	b = kann_new_vec(n1, 1.0f); /* see Jozefowicz et al on using a large bias */
	f = kad_sigm(kad_add(kad_add(cmul(in, w), cmul(h0, u)), b));
	/* o = sigm(x_t * W_o + h_{t-1} * U_o + b_o) */
	w = kann_new_weight(n1, n0);
	u = kann_new_weight(n1, n1);
	b = kann_new_bias(n1);
	o = kad_sigm(kad_add(kad_add(cmul(in, w), cmul(h0, u)), b));
	/* g = tanh(x_t * W_g + h_{t-1} * U_g + b_g) */
	w = kann_new_weight(n1, n0);
	u = kann_new_weight(n1, n1);
	b = kann_new_bias(n1);
	g = kad_tanh(kad_add(kad_add(cmul(in, w), cmul(h0, u)), b));
	/* c_t = c_{t-1} # f + g # i */
	c = kad_add(kad_mul(f, c0), kad_mul(g, i)); /* can't be kad_mul(c0, f)!!! */
	c->pre = c0;
	/* h_t = tanh(c_t) # o */
	if (rnn_flag & KANN_RNN_NORM) c = kann_layer_layernorm(c); /* see Ba et al (2016) about how to apply layer normalization to LSTM */
	out = kad_mul(kad_tanh(c), o);
	out->pre = h0;
	return out;
}

kad_node_t *kann_layer_conv2d(kad_node_t *in, int n_flt, int k_rows, int k_cols, int stride_r, int stride_c, int pad_r, int pad_c)
{
	kad_node_t *w;
	w = kann_new_weight_conv2d(n_flt, in->d[1], k_rows, k_cols);
	return kad_conv2d(in, w, stride_r, stride_c, pad_r, pad_c);
}

kad_node_t *kann_layer_conv1d(kad_node_t *in, int n_flt, int k_size, int stride, int pad)
{
	kad_node_t *w;
	w = kann_new_weight_conv1d(n_flt, in->d[1], k_size);
	return kad_conv1d(in, w, stride, pad);
}

kad_node_t *kann_layer_cost(kad_node_t *t, int n_out, int cost_type)
{
	kad_node_t *cost = 0, *truth = 0;
	assert(cost_type == KANN_C_CEB || cost_type == KANN_C_CEM || cost_type == KANN_C_CEB_NEG || cost_type == KANN_C_MSE);
	t = kann_layer_dense(t, n_out);
	truth = kad_feed(2, 1, n_out), truth->ext_flag |= KANN_F_TRUTH;

	if (cost_type == KANN_C_MSE) {
		cost = kad_mse(t, truth);
	} else if (cost_type == KANN_C_CEB) {
		t = kad_sigm(t);
		cost = kad_ce_bin(t, truth);
	} else if (cost_type == KANN_C_CEB_NEG) {
		t = kad_tanh(t);
		cost = kad_ce_bin_neg(t, truth);
	} else if (cost_type == KANN_C_CEM) {
		t = kad_softmax(t);
		cost = kad_ce_multi(t, truth);
	}
	else {
		assert (0);
	}

	t->ext_flag |= KANN_F_OUT;
	cost->ext_flag |= KANN_F_COST;

	return cost;
}

void kann_shuffle(int n, int *s)
{
	int i, j, t;
	for (i = 0; i < n; ++i) s[i] = i;
	for (i = n; i > 0; --i) {
		j = (int)(i * kad_drand(0));
		t = s[j], s[j] = s[i-1], s[i-1] = t;
	}
}

/***************************
 *** @@MIN: minimization ***
 ***************************/

#ifdef __SSE__
#include <xmmintrin.h>

void kann_RMSprop(int n, float h0, const float *h, float decay, const float *g, float *t, float *r)
{
	int i, n4 = n>>2<<2;
	__m128 vh, vg, vr, vt, vd, vd1, tmp, vtiny;
	vh = _mm_set1_ps(h0);
	vd = _mm_set1_ps(decay);
	vd1 = _mm_set1_ps(1.0f - decay);
	vtiny = _mm_set1_ps(1e-6f);
	for (i = 0; i < n4; i += 4) {
		vt = _mm_loadu_ps(&t[i]);
		vr = _mm_loadu_ps(&r[i]);
		vg = _mm_loadu_ps(&g[i]);
		if (h) vh = _mm_loadu_ps(&h[i]);
		vr = _mm_add_ps(_mm_mul_ps(vd1, _mm_mul_ps(vg, vg)), _mm_mul_ps(vd, vr));
		_mm_storeu_ps(&r[i], vr);
		tmp = _mm_sub_ps(vt, _mm_mul_ps(_mm_mul_ps(vh, _mm_rsqrt_ps(_mm_add_ps(vtiny, vr))), vg));
		_mm_storeu_ps(&t[i], tmp);
	}
	for (; i < n; ++i) {
		r[i] = (1. - decay) * g[i] * g[i] + decay * r[i];
		t[i] -= (h? h[i] : h0) / sqrtf(1e-6f + r[i]) * g[i];
	}
}
#else
void kann_RMSprop(int n, float h0, const float *h, float decay, const float *g, float *t, float *r)
{
	int i;
	for (i = 0; i < n; ++i) {
		float lr = h? h[i] : h0;
		r[i] = (1.0f - decay) * g[i] * g[i] + decay * r[i];
		t[i] -= lr / sqrtf(1e-6f + r[i]) * g[i];
	}
}
#endif

float kann_grad_clip(float thres, int n, float *g)
{
	int i;
	double s2 = 0.0;
	for (i = 0; i < n; ++i)
		s2 += g[i] * g[i];
	s2 = sqrt(s2);
	if (s2 > thres)
		for (i = 0, s2 = 1.0 / s2; i < n; ++i)
			g[i] *= (float)s2;
	return (float)s2 / thres;
}

/****************************************************************
 *** @@XY: simpler API for network with a single input/output ***
 ****************************************************************/

int kann_train_fnn1(kann_t *ann, float lr, int mini_size, int max_epoch,
					int max_drop_streak, float frac_val, int n,
					float **_x, float **_y,
					kann_train_cb cb, void *ud)
{
	int i, j, *shuf, n_train, n_val, n_in, n_out, n_var, n_const, drop_streak = 0, min_set = 0;
	float **x, **y, *x1, *y1, *r, min_val_cost = FLT_MAX, *min_x, *min_c;

	n_in = kann_dim_in(ann);
	n_out = kann_dim_out(ann);
	if (n_in < 0 || n_out < 0) return -1;
	n_var = kann_size_var(ann);
	n_const = kann_size_const(ann);
	r = (float*)calloc(n_var, sizeof(float));
	shuf = (int*)malloc(n * sizeof(int));
	x = (float**)malloc(n * sizeof(float*));
	y = (float**)malloc(n * sizeof(float*));
	kann_shuffle(n, shuf);
	for (j = 0; j < n; ++j)
		x[j] = _x[shuf[j]], y[j] = _y[shuf[j]];
	n_val = (int)(n * frac_val);
	n_train = n - n_val;
	min_x = (float*)malloc(n_var * sizeof(float));
	min_c = (float*)malloc(n_const * sizeof(float));

	x1 = (float*)malloc(n_in  * mini_size * sizeof(float));
	y1 = (float*)malloc(n_out * mini_size * sizeof(float));
	kann_feed_bind(ann, KANN_F_IN,    0, &x1);
	kann_feed_bind(ann, KANN_F_TRUTH, 0, &y1);

	for (i = 0; i < max_epoch; ++i) {
		int n_proc = 0, n_train_err = 0, n_val_err = 0, n_train_base = 0, n_val_base = 0;
		double train_cost = 0.0, val_cost = 0.0;
		kann_shuffle(n_train, shuf);
		kann_switch(ann, 1);
		while (n_proc < n_train) {
			int b, c, ms = n_train - n_proc < mini_size? n_train - n_proc : mini_size;
			for (b = 0; b < ms; ++b) {
				memcpy(&x1[b*n_in],  x[shuf[n_proc+b]], n_in  * sizeof(float));
				memcpy(&y1[b*n_out], y[shuf[n_proc+b]], n_out * sizeof(float));
			}
			kann_set_batch_size(ann, ms);
			train_cost += kann_cost(ann, 0, 1) * ms;
			c = kann_class_error(ann, &b);
			n_train_err += c, n_train_base += b;
			kann_RMSprop(n_var, lr, 0, 0.9f, ann->g, ann->x, r);
			n_proc += ms;
		}
		train_cost /= n_train;
		kann_switch(ann, 0);
		n_proc = 0;
		while (n_proc < n_val) {
			int b, c, ms = n_val - n_proc < mini_size? n_val - n_proc : mini_size;
			for (b = 0; b < ms; ++b) {
				memcpy(&x1[b*n_in],  x[n_train+n_proc+b], n_in  * sizeof(float));
				memcpy(&y1[b*n_out], y[n_train+n_proc+b], n_out * sizeof(float));
			}
			kann_set_batch_size(ann, ms);
			val_cost += kann_cost(ann, 0, 0) * ms;
			c = kann_class_error(ann, &b);
			n_val_err += c, n_val_base += b;
			n_proc += ms;
		}
		if (n_val > 0) val_cost /= n_val;
		if (cb) {
			cb(i + 1, train_cost, val_cost, ud);
#if 0
			fprintf(stderr, "epoch: %d; training cost: %g", i+1, train_cost);
			if (n_train_base) fprintf(stderr, " (class error: %.2f%%)", 100.0f * n_train_err / n_train);
			if (n_val > 0) {
				fprintf(stderr, "; validation cost: %g", val_cost);
				if (n_val_base) fprintf(stderr, " (class error: %.2f%%)", 100.0f * n_val_err / n_val);
			}
			fputc('\n', stderr);
#endif
		}
		if (i >= max_drop_streak && n_val > 0) {
			if (val_cost < min_val_cost) {
				min_set = 1;
				memcpy(min_x, ann->x, n_var * sizeof(float));
				memcpy(min_c, ann->c, n_const * sizeof(float));
				drop_streak = 0;
				min_val_cost = (float)val_cost;
			} else if (++drop_streak >= max_drop_streak)
				break;
		}
	}
	if (min_set) {
		memcpy(ann->x, min_x, n_var * sizeof(float));
		memcpy(ann->c, min_c, n_const * sizeof(float));
	}

	free(min_c); free(min_x); free(y1); free(x1); free(y); free(x); free(shuf); free(r);
	return i;
}

float kann_cost_fnn1(kann_t *ann, int n, float **x, float **y)
{
	int n_in, n_out, n_proc = 0, mini_size = 64 < n? 64 : n;
	float *x1, *y1;
	double cost = 0.0;

	n_in = kann_dim_in(ann);
	n_out = kann_dim_out(ann);
	if (n <= 0 || n_in < 0 || n_out < 0) return 0.0;

	x1 = (float*)malloc(n_in  * mini_size * sizeof(float));
	y1 = (float*)malloc(n_out * mini_size * sizeof(float));
	kann_feed_bind(ann, KANN_F_IN,    0, &x1);
	kann_feed_bind(ann, KANN_F_TRUTH, 0, &y1);
	kann_switch(ann, 0);
	while (n_proc < n) {
		int b, ms = n - n_proc < mini_size? n - n_proc : mini_size;
		for (b = 0; b < ms; ++b) {
			memcpy(&x1[b*n_in],  x[n_proc+b], n_in  * sizeof(float));
			memcpy(&y1[b*n_out], y[n_proc+b], n_out * sizeof(float));
		}
		kann_set_batch_size(ann, ms);
		cost += kann_cost(ann, 0, 0) * ms;
		n_proc += ms;
	}
	free(y1); free(x1);
	return (float)(cost / n);
}

const float *kann_apply1(kann_t *a, float *x)
{
	int i_out;
	i_out = kann_find(a, KANN_F_OUT, 0);
	if (i_out < 0) return 0;
	kann_set_batch_size(a, 1);
	kann_feed_bind(a, KANN_F_IN, 0, &x);
	kad_eval_at(a->n, a->v, i_out);
	return a->v[i_out]->x;
}