alpha = new double[alpha_.Length];
			alpha_.CopyTo(alpha, 0);
			this.Cp = Cp;
			this.Cn = Cn;
			this.eps = eps;
			this.unshrinked = false;
			
			// initialize alpha_status
			{
				alpha_status = new sbyte[l];
				for (int i = 0; i < l; i++)
					update_alpha_status(i);
			}
			
			// initialize active set (for shrinking)
			{
				active_set = new int[l];
				for (int i = 0; i < l; i++)
					active_set[i] = i;
				active_size = l;
			}
			
			// initialize gradient
			{
				G = new double[l];
				G_bar = new double[l];
				int i;
				for (i = 0; i < l; i++)
				{
					G[i] = b[i];
					G_bar[i] = 0;
				}
				for (i = 0; i < l; i++)
					if (!is_lower_bound(i))
					{
						float[] Q_i = Q.get_Q(i, l);
						double alpha_i = alpha[i];
						int j;
						for (j = 0; j < l; j++)
							G[j] += alpha_i * Q_i[j];
						if (is_upper_bound(i))
							for (j = 0; j < l; j++)
								G_bar[j] += get_C(i) * Q_i[j];
					}
			}
			
			// optimization step
			
			int iter = 0;
			int counter = System.Math.Min(l, 1000) + 1;
			int[] working_set = new int[2];
			
			while (true)
			{
				// show progress and do shrinking
				
				if (--counter == 0)
				{
					counter = System.Math.Min(l, 1000);
					if (shrinking != 0)
						do_shrinking();
					System.Console.Error.Write(".");
				}
				
				if (select_working_set(working_set) != 0)
				{
					// reconstruct the whole gradient
					reconstruct_gradient();
					// reset active set size and check
					active_size = l;
					System.Console.Error.Write("*");
					if (select_working_set(working_set) != 0)
						break;
					else
						counter = 1; // do shrinking next iteration
				}
				
				int i = working_set[0];
				int j = working_set[1];
				
				++iter;
				
				// update alpha[i] and alpha[j], handle bounds carefully
				
				float[] Q_i = Q.get_Q(i, active_size);
				float[] Q_j = Q.get_Q(j, active_size);
				
				double C_i = get_C(i);
				double C_j = get_C(j);
				
				double old_alpha_i = alpha[i];
				double old_alpha_j = alpha[j];
				
				if (y[i] != y[j])
				{
					double delta = (- G[i] - G[j]) / System.Math.Max(Q_i[i] + Q_j[j] + 2 * Q_i[j], (float) 0);
					double diff = alpha[i] - alpha[j];
					alpha[i] += delta;
					alpha[j] += delta;
					
					if (diff > 0)
					{
						if (alpha[j] < 0)
						{
							alpha[j] = 0;
							alpha[i] = diff;
						}
					}
					else
					{
						if (alpha[i] < 0)
						{
							alpha[i] = 0;
							alpha[j] = - diff;
						}
					}
					if (diff > C_i - C_j)
					{
						if (alpha[i] > C_i)
						{
							alpha[i] = C_i;
							alpha[j] = C_i - diff;
						}
					}
					else
					{
						if (alpha[j] > C_j)
						{
							alpha[j] = C_j;
							alpha[i] = C_j + diff;
						}
					}
				}
				else
				{
					double delta = (G[i] - G[j]) / System.Math.Max(Q_i[i] + Q_j[j] - 2 * Q_i[j], (float) 0);
					double sum = alpha[i] + alpha[j];
					alpha[i] -= delta;
					alpha[j] += delta;
					if (sum > C_i)
					{
						if (alpha[i] > C_i)
						{
							alpha[i] = C_i;
							alpha[j] = sum - C_i;
						}
					}
					else
					{
						if (alpha[j] < 0)
						{
							alpha[j] = 0;
							alpha[i] = sum;
						}
					}
					if (sum > C_j)
					{
						if (alpha[j] > C_j)
						{
							alpha[j] = C_j;
							alpha[i] = sum - C_j;
						}
					}
					else
					{
						if (alpha[i] < 0)
						{
							alpha[i] = 0;
							alpha[j] = sum;
						}
					}
				}
				
				// update G
				
				double delta_alpha_i = alpha[i] - old_alpha_i;
				double delta_alpha_j = alpha[j] - old_alpha_j;
				
				for (int k = 0; k < active_size; k++)
				{
					G[k] += Q_i[k] * delta_alpha_i + Q_j[k] * delta_alpha_j;
				}
				
				// update alpha_status and G_bar
				
				{
					bool ui = is_upper_bound(i);
					bool uj = is_upper_bound(j);
					update_alpha_status(i);
					update_alpha_status(j);
					int k;
					if (ui != is_upper_bound(i))
					{
						Q_i = Q.get_Q(i, l);
						if (ui)
							for (k = 0; k < l; k++)
								G_bar[k] -= C_i * Q_i[k];
						else
							for (k = 0; k < l; k++)
								G_bar[k] += C_i * Q_i[k];
					}
					
					if (uj != is_upper_bound(j))
					{
						Q_j = Q.get_Q(j, l);
						if (uj)
							for (k = 0; k < l; k++)
								G_bar[k] -= C_j * Q_j[k];
						else
							for (k = 0; k < l; k++)
								G_bar[k] += C_j * Q_j[k];
					}
				}
			}
			
			// calculate rho
			
			si.rho = calculate_rho();
			
			// calculate objective value
			{
				double v = 0;
				int i;
				for (i = 0; i < l; i++)
					v += alpha[i] * (G[i] + b[i]);
				
				si.obj = v / 2;
			}
			
			// put back the solution
			{
				for (int i = 0; i < l; i++)
					alpha_[active_set[i]] = alpha[i];
			}
			
			si.upper_bound_p = Cp;
			si.upper_bound_n = Cn;
			
			System.Console.Out.Write("\noptimization finished, #iter = " + iter + "\n");
		}
		
		// return 1 if already optimal, return 0 otherwise
		internal virtual int select_working_set(int[] working_set)
		{
			// return i,j which maximize -grad(f)^T d , under constraint
			// if alpha_i == C, d != +1
			// if alpha_i == 0, d != -1
			
			double Gmax1 = - INF; // max { -grad(f)_i * d | y_i*d = +1 }
			int Gmax1_idx = - 1;
			
			double Gmax2 = - INF; // max { -grad(f)_i * d | y_i*d = -1 }
			int Gmax2_idx = - 1;
			
			for (int i = 0; i < active_size; i++)
			{
				if (y[i] == + 1)
				// y = +1
				{
					if (!is_upper_bound(i))
					// d = +1
					{
						if (- G[i] > Gmax1)
						{
							Gmax1 = - G[i];
							Gmax1_idx = i;
						}
					}
					if (!is_lower_bound(i))
					// d = -1
					{
						if (G[i] > Gmax2)
						{
							Gmax2 = G[i];
							Gmax2_idx = i;
						}
					}
				}
				// y = -1
				else
				{
					if (!is_upper_bound(i))
					// d = +1
					{
						if (- G[i] > Gmax2)
						{
							Gmax2 = - G[i];
							Gmax2_idx = i;
						}
					}
					if (!is_lower_bound(i))
					// d = -1
					{
						if (G[i] > Gmax1)
						{
							Gmax1 = G[i];
							Gmax1_idx = i;
						}
					}
				}
			}
			
			if (Gmax1 + Gmax2 < eps)
				return 1;
			
			working_set[0] = Gmax1_idx;
			working_set[1] = Gmax2_idx;
			return 0;
		}
		
		internal virtual void  do_shrinking()
		{
			int i, j, k;
			int[] working_set = new int[2];
			if (select_working_set(working_set) != 0)
				return ;
			i = working_set[0];
			j = working_set[1];
			double Gm1 = (- y[j]) * G[j];
			double Gm2 = y[i] * G[i];
			
			// shrink
			
			for (k = 0; k < active_size; k++)
			{
				if (is_lower_bound(k))
				{
					if (y[k] == + 1)
					{
						if (- G[k] >= Gm1)
							continue;
					}
					else if (- G[k] >= Gm2)
						continue;
				}
				else if (is_upper_bound(k))
				{
					if (y[k] == + 1)
					{
						if (G[k] >= Gm2)
							continue;
					}
					else if (G[k] >= Gm1)
						continue;
				}
				else
					continue;
				
				--active_size;
				swap_index(k, active_size);
				--k; // look at the newcomer
			}
			
			// unshrink, check all variables again before final iterations
			
			if (unshrinked || - (Gm1 + Gm2) > eps * 10)
				return ;
			
			unshrinked = true;
			reconstruct_gradient();
			
			for (k = l - 1; k >= active_size; k--)
			{
				if (is_lower_bound(k))
				{
					if (y[k] == + 1)
					{
						if (- G[k] < Gm1)
							continue;
					}
					else if (- G[k] < Gm2)
						continue;
				}
				else if (is_upper_bound(k))
				{
					if (y[k] == + 1)
					{
						if (G[k] < Gm2)
							continue;
					}
					else if (G[k] < Gm1)
						continue;
				}
				else
					continue;
				
				swap_index(k, active_size);
				active_size++;
				++k; // look at the newcomer
			}
		}
		
		internal virtual double calculate_rho()
		{
			double r;
			int nr_free = 0;
			double ub = INF, lb = - INF, sum_free = 0;
			for (int i = 0; i < active_size; i++)
			{
				double yG = y[i] * G[i];
				
				if (is_lower_bound(i))
				{
					if (y[i] > 0)
						ub = System.Math.Min(ub, yG);
					else
						lb = System.Math.Max(lb, yG);
				}
				else if (is_upper_bound(i))
				{
					if (y[i] < 0)
						ub = System.Math.Min(ub, yG);
					else
						lb = System.Math.Max(lb, yG);
				}
				else
				{
					++nr_free;
					sum_free += yG;
				}
			}
			
			if (nr_free > 0)
				r = sum_free / nr_free;
			else
				r = (ub + lb) / 2;
			
			return r;
		}
	}
	
	//
	// Solver for nu-svm classification and regression
	//
	// additional constraint: e^T \alpha = constant
	//
	sealed class Solver_NU:Solver
	{
		private SolutionInfo si;
		
		internal override void  Solve(int l, Kernel Q, double[] b, sbyte[] y, double[] alpha, double Cp, double Cn, double eps, SolutionInfo si, int shrinking)
		{
			this.si = si;
			base.Solve(l, Q, b, y, alpha, Cp, Cn, eps, si, shrinking);
		}
		
		internal override int select_working_set(int[] working_set)
		{
			// return i,j which maximize -grad(f)^T d , under constraint
			// if alpha_i == C, d != +1
			// if alpha_i == 0, d != -1
			
			double Gmax1 = - INF; // max { -grad(f)_i * d | y_i = +1, d = +1 }
			int Gmax1_idx = - 1;
			
			double Gmax2 = - INF; // max { -grad(f)_i * d | y_i = +1, d = -1 }
			int Gmax2_idx = - 1;
			
			double Gmax3 = - INF; // max { -grad(f)_i * d | y_i = -1, d = +1 }
			int Gmax3_idx = - 1;
			
			double Gmax4 = - INF; // max { -grad(f)_i * d | y_i = -1, d = -1 }
			int Gmax4_idx = - 1;
			
			for (int i = 0; i < active_size; i++)
			{
				if (y[i] == + 1)
				// y == +1
				{
					if (!is_upper_bound(i))
					// d = +1
					{
						if (- G[i] > Gmax1)
						{
							Gmax1 = - G[i];
							Gmax1_idx = i;
						}
					}
					if (!is_lower_bound(i))
					// d = -1
					{
						if (G[i] > Gmax2)
						{
							Gmax2 = G[i];
							Gmax2_idx = i;
						}
					}
				}
				// y == -1
				else
				{
					if (!is_upper_bound(i))
					// d = +1