/************************************************************************* AP library Copyright (c) 2003-2009 Sergey Bochkanov (ALGLIB project). >>> SOURCE LICENSE >>> This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation (www.fsf.org); either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. A copy of the GNU General Public License is available at http://www.fsf.org/licensing/licenses >>> END OF LICENSE >>> *************************************************************************/ using System; public partial class alglib { /******************************************************************** Callback definitions for optimizers/fitters/solvers. Callbacks for unparameterized (general) functions: * ndimensional_func calculates f(arg), stores result to func * ndimensional_grad calculates func = f(arg), grad[i] = df(arg)/d(arg[i]) * ndimensional_hess calculates func = f(arg), grad[i] = df(arg)/d(arg[i]), hess[i,j] = d2f(arg)/(d(arg[i])*d(arg[j])) Callbacks for systems of functions: * ndimensional_fvec calculates vector function f(arg), stores result to fi * ndimensional_jac calculates f[i] = fi(arg) jac[i,j] = df[i](arg)/d(arg[j]) Callbacks for parameterized functions, i.e. for functions which depend on two vectors: P and Q. Gradient and Hessian are calculated with respect to P only. * ndimensional_pfunc calculates f(p,q), stores result to func * ndimensional_pgrad calculates func = f(p,q), grad[i] = df(p,q)/d(p[i]) * ndimensional_phess calculates func = f(p,q), grad[i] = df(p,q)/d(p[i]), hess[i,j] = d2f(p,q)/(d(p[i])*d(p[j])) Callbacks for progress reports: * ndimensional_rep reports current position of optimization algo Callbacks for ODE solvers: * ndimensional_ode_rp calculates dy/dx for given y[] and x Callbacks for integrators: * integrator1_func calculates f(x) for given x (additional parameters xminusa and bminusx contain x-a and b-x) ********************************************************************/ public delegate void ndimensional_func (double[] arg, ref double func, object obj); public delegate void ndimensional_grad (double[] arg, ref double func, double[] grad, object obj); public delegate void ndimensional_hess (double[] arg, ref double func, double[] grad, double[,] hess, object obj); public delegate void ndimensional_fvec (double[] arg, double[] fi, object obj); public delegate void ndimensional_jac (double[] arg, double[] fi, double[,] jac, object obj); public delegate void ndimensional_pfunc(double[] p, double[] q, ref double func, object obj); public delegate void ndimensional_pgrad(double[] p, double[] q, ref double func, double[] grad, object obj); public delegate void ndimensional_phess(double[] p, double[] q, ref double func, double[] grad, double[,] hess, object obj); public delegate void ndimensional_rep(double[] arg, double func, object obj); public delegate void ndimensional_ode_rp (double[] y, double x, double[] dy, object obj); public delegate void integrator1_func (double x, double xminusa, double bminusx, ref double f, object obj); /******************************************************************** Class defining a complex number with double precision. ********************************************************************/ public struct complex { public double x; public double y; public complex(double _x) { x = _x; y = 0; } public complex(double _x, double _y) { x = _x; y = _y; } public static implicit operator complex(double _x) { return new complex(_x); } public static bool operator==(complex lhs, complex rhs) { return ((double)lhs.x==(double)rhs.x) & ((double)lhs.y==(double)rhs.y); } public static bool operator!=(complex lhs, complex rhs) { return ((double)lhs.x!=(double)rhs.x) | ((double)lhs.y!=(double)rhs.y); } public static complex operator+(complex lhs) { return lhs; } public static complex operator-(complex lhs) { return new complex(-lhs.x,-lhs.y); } public static complex operator+(complex lhs, complex rhs) { return new complex(lhs.x+rhs.x,lhs.y+rhs.y); } public static complex operator-(complex lhs, complex rhs) { return new complex(lhs.x-rhs.x,lhs.y-rhs.y); } public static complex operator*(complex lhs, complex rhs) { return new complex(lhs.x*rhs.x-lhs.y*rhs.y, lhs.x*rhs.y+lhs.y*rhs.x); } public static complex operator/(complex lhs, complex rhs) { complex result; double e; double f; if( System.Math.Abs(rhs.y)(T[] a) { return a.Length; } public static int rows(T[,] a) { return a.GetLength(0); } public static int cols(T[,] a) { return a.GetLength(1); } public static void swap(ref T a, ref T b) { T t = a; a = b; b = t; } public static void assert(bool cond, string s) { if( !cond ) throw new alglibexception(s); } public static void assert(bool cond) { assert(cond, "ALGLIB: assertion failed"); } /**************************************************************** returns dps (digits-of-precision) value corresponding to threshold. dps(0.9) = dps(0.5) = dps(0.1) = 0 dps(0.09) = dps(0.05) = dps(0.01) = 1 and so on ****************************************************************/ public static int threshold2dps(double threshold) { int result = 0; double t; for (result = 0, t = 1; t / 10 > threshold*(1+1E-10); result++, t /= 10) ; return result; } /**************************************************************** prints formatted complex ****************************************************************/ public static string format(complex a, int _dps) { int dps = Math.Abs(_dps); string fmt = _dps>=0 ? "F" : "E"; string fmtx = String.Format("{{0:"+fmt+"{0}}}", dps); string fmty = String.Format("{{0:"+fmt+"{0}}}", dps); string result = String.Format(fmtx, a.x) + (a.y >= 0 ? "+" : "-") + String.Format(fmty, Math.Abs(a.y)) + "i"; result = result.Replace(',', '.'); return result; } /**************************************************************** prints formatted array ****************************************************************/ public static string format(bool[] a) { string[] result = new string[len(a)]; int i; for(i=0; i= 0 ? "F" : "E"; string fmt = String.Format("{{0:" + sfmt + "{0}}}", dps); string[] result = new string[len(a)]; int i; for (i = 0; i < len(a); i++) { result[i] = String.Format(fmt, a[i]); result[i] = result[i].Replace(',', '.'); } return "{" + String.Join(",", result) + "}"; } /**************************************************************** prints formatted array ****************************************************************/ public static string format(complex[] a, int _dps) { int dps = Math.Abs(_dps); string fmt = _dps >= 0 ? "F" : "E"; string fmtx = String.Format("{{0:"+fmt+"{0}}}", dps); string fmty = String.Format("{{0:"+fmt+"{0}}}", dps); string[] result = new string[len(a)]; int i; for (i = 0; i < len(a); i++) { result[i] = String.Format(fmtx, a[i].x) + (a[i].y >= 0 ? "+" : "-") + String.Format(fmty, Math.Abs(a[i].y)) + "i"; result[i] = result[i].Replace(',', '.'); } return "{" + String.Join(",", result) + "}"; } /**************************************************************** prints formatted matrix ****************************************************************/ public static string format(bool[,] a) { int i, j, m, n; n = cols(a); m = rows(a); bool[] line = new bool[n]; string[] result = new string[m]; for (i = 0; i < m; i++) { for (j = 0; j < n; j++) line[j] = a[i, j]; result[i] = format(line); } return "{" + String.Join(",", result) + "}"; } /**************************************************************** prints formatted matrix ****************************************************************/ public static string format(int[,] a) { int i, j, m, n; n = cols(a); m = rows(a); int[] line = new int[n]; string[] result = new string[m]; for (i = 0; i < m; i++) { for (j = 0; j < n; j++) line[j] = a[i, j]; result[i] = format(line); } return "{" + String.Join(",", result) + "}"; } /**************************************************************** prints formatted matrix ****************************************************************/ public static string format(double[,] a, int dps) { int i, j, m, n; n = cols(a); m = rows(a); double[] line = new double[n]; string[] result = new string[m]; for (i = 0; i < m; i++) { for (j = 0; j < n; j++) line[j] = a[i, j]; result[i] = format(line, dps); } return "{" + String.Join(",", result) + "}"; } /**************************************************************** prints formatted matrix ****************************************************************/ public static string format(complex[,] a, int dps) { int i, j, m, n; n = cols(a); m = rows(a); complex[] line = new complex[n]; string[] result = new string[m]; for (i = 0; i < m; i++) { for (j = 0; j < n; j++) line[j] = a[i, j]; result[i] = format(line, dps); } return "{" + String.Join(",", result) + "}"; } /**************************************************************** checks that matrix is symmetric. max|A-A^T| is calculated; if it is within 1.0E-14 of max|A|, matrix is considered symmetric ****************************************************************/ public static bool issymmetric(double[,] a) { int i, j, n; double err, mx, v1, v2; if( rows(a)!=cols(a) ) return false; n = rows(a); if( n==0 ) return true; mx = 0; err = 0; for( i=0; iyabs ? xabs : yabs; v = xabs63 ) return '?'; return _sixbits2char_tbl[v]; } /************************************************************************ This function converts character to six-bit value (from 0 to 63). This function is inverse of ae_sixbits2char() If c is not correct character, this function returns -1. ************************************************************************/ private static int[] _char2sixbits_tbl = new int[128] { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 62, -1, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, -1, -1, -1, -1, -1, -1, -1, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, -1, -1, -1, -1, 63, -1, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, -1, -1, -1, -1, -1 }; private static int char2sixbits(char c) { return (c>=0 && c<127) ? _char2sixbits_tbl[c] : -1; } /************************************************************************ This function converts three bytes (24 bits) to four six-bit values (24 bits again). src array src_offs offset of three-bytes chunk dst array for ints dst_offs offset of four-ints chunk ************************************************************************/ private static void threebytes2foursixbits(byte[] src, int src_offs, int[] dst, int dst_offs) { dst[dst_offs+0] = src[src_offs+0] & 0x3F; dst[dst_offs+1] = (src[src_offs+0]>>6) | ((src[src_offs+1]&0x0F)<<2); dst[dst_offs+2] = (src[src_offs+1]>>4) | ((src[src_offs+2]&0x03)<<4); dst[dst_offs+3] = src[src_offs+2]>>2; } /************************************************************************ This function converts four six-bit values (24 bits) to three bytes (24 bits again). src pointer to four ints src_offs offset of the chunk dst pointer to three bytes dst_offs offset of the chunk ************************************************************************/ private static void foursixbits2threebytes(int[] src, int src_offs, byte[] dst, int dst_offs) { dst[dst_offs+0] = (byte)(src[src_offs+0] | ((src[src_offs+1]&0x03)<<6)); dst[dst_offs+1] = (byte)((src[src_offs+1]>>2) | ((src[src_offs+2]&0x0F)<<4)); dst[dst_offs+2] = (byte)((src[src_offs+2]>>4) | (src[src_offs+3]<<2)); } /************************************************************************ This function serializes boolean value into buffer v boolean value to be serialized buf buffer, at least 11 characters wide offs offset in the buffer after return from this function, offs points to the char's past the value being read. ************************************************************************/ private static void bool2str(bool v, char[] buf, ref int offs) { char c = v ? '1' : '0'; int i; for(i=0; i=SER_ENTRY_LENGTH ) throw new alglib.alglibexception(emsg); sixbits[sixbitsread] = d; sixbitsread++; offs++; } if( sixbitsread==0 ) throw new alglib.alglibexception(emsg); for(i=sixbitsread; i<12; i++) sixbits[i] = 0; foursixbits2threebytes(sixbits, 0, bytes, 0); foursixbits2threebytes(sixbits, 4, bytes, 3); foursixbits2threebytes(sixbits, 8, bytes, 6); c = (bytes[sizeof(int)-1] & 0x80)!=0 ? (byte)0xFF : (byte)0x00; for(i=sizeof(int); i<8; i++) if( bytes[i]!=c ) throw new alglib.alglibexception(emsg3264); for(i=0; i=SER_ENTRY_LENGTH ) throw new alglib.alglibexception(emsg); sixbits[sixbitsread] = d; sixbitsread++; offs++; } if( sixbitsread!=SER_ENTRY_LENGTH ) throw new alglib.alglibexception(emsg); sixbits[SER_ENTRY_LENGTH] = 0; foursixbits2threebytes(sixbits, 0, bytes, 0); foursixbits2threebytes(sixbits, 4, bytes, 3); foursixbits2threebytes(sixbits, 8, bytes, 6); for(i=0; i(shared_pool pool, ref T obj) where T : alglib.apobject { alglib.apobject new_obj; /* assert that pool was seeded */ alglib.ap.assert(pool.seed_object!=null, "ALGLIB: shared pool is not seeded, PoolRetrieve() failed"); /* acquire lock */ ae_acquire_lock(pool.pool_lock); /* try to reuse recycled objects */ if( pool.recycled_objects!=null ) { /* retrieve entry/object from list of recycled objects */ sharedpoolentry result = pool.recycled_objects; pool.recycled_objects = pool.recycled_objects.next_entry; new_obj = result.obj; result.obj = null; /* move entry to list of recycled entries */ result.next_entry = pool.recycled_entries; pool.recycled_entries = result; /* release lock */ ae_release_lock(pool.pool_lock); /* assign object to smart pointer */ obj = (T)new_obj; return; } /* * release lock; we do not need it anymore because * copy constructor does not modify source variable. */ ae_release_lock(pool.pool_lock); /* create new object from seed */ new_obj = pool.seed_object.make_copy(); /* assign object to pointer and return */ obj = (T)new_obj; } /************************************************************************ This function recycles object owned by the source variable by moving it to internal storage of the shared pool. Source variable must own the object, i.e. be the only place where reference to object is stored. After call to this function source variable becomes NULL. pool pool obj source variable NOTE: this function IS thread-safe. It acquires pool lock during its operation and can be used simultaneously from several threads. ************************************************************************/ public static void ae_shared_pool_recycle(shared_pool pool, ref T obj) where T : alglib.apobject { sharedpoolentry new_entry; /* assert that pool was seeded */ alglib.ap.assert(pool.seed_object!=null, "ALGLIB: shared pool is not seeded, PoolRecycle() failed"); /* assert that pointer non-null */ alglib.ap.assert(obj!=null, "ALGLIB: obj in ae_shared_pool_recycle() is NULL"); /* acquire lock */ ae_acquire_lock(pool.pool_lock); /* acquire shared pool entry (reuse one from recycled_entries or malloc new one) */ if( pool.recycled_entries!=null ) { /* reuse previously allocated entry */ new_entry = pool.recycled_entries; pool.recycled_entries = new_entry.next_entry; } else { /* * Allocate memory for new entry. * * NOTE: we release pool lock during allocation because new() may raise * exception and we do not want our pool to be left in the locked state. */ ae_release_lock(pool.pool_lock); new_entry = new sharedpoolentry(); ae_acquire_lock(pool.pool_lock); } /* add object to the list of recycled objects */ new_entry.obj = obj; new_entry.next_entry = pool.recycled_objects; pool.recycled_objects = new_entry; /* release lock object */ ae_release_lock(pool.pool_lock); /* release source pointer */ obj = null; } /************************************************************************ This function clears internal list of recycled objects, but does not change seed object managed by the pool. pool pool NOTE: this function is NOT thread-safe. It does not acquire pool lock, so you should NOT call it when lock can be used by another thread. ************************************************************************/ public static void ae_shared_pool_clear_recycled(shared_pool pool) { pool.recycled_objects = null; } /************************************************************************ This function allows to enumerate recycled elements of the shared pool. It stores reference to the first recycled object in the smart pointer. IMPORTANT: * in case target variable owns non-NULL value, it is rewritten * recycled object IS NOT removed from pool * target variable DOES NOT become owner of the new value; you can use reference to recycled object, but you do not own it. * this function IS NOT thread-safe * you SHOULD NOT modify shared pool during enumeration (although you can modify state of the objects retrieved from pool) * in case there is no recycled objects in the pool, NULL is stored to obj * in case pool is not seeded, NULL is stored to obj pool pool obj reference ************************************************************************/ public static void ae_shared_pool_first_recycled(shared_pool pool, ref T obj) where T : alglib.apobject { /* modify internal enumeration counter */ pool.enumeration_counter = pool.recycled_objects; /* exit on empty list */ if( pool.enumeration_counter==null ) { obj = null; return; } /* assign object to smart pointer */ obj = (T)pool.enumeration_counter.obj; } /************************************************************************ This function allows to enumerate recycled elements of the shared pool. It stores pointer to the next recycled object in the smart pointer. IMPORTANT: * in case target variable owns non-NULL value, it is rewritten * recycled object IS NOT removed from pool * target pointer DOES NOT become owner of the new value * this function IS NOT thread-safe * you SHOULD NOT modify shared pool during enumeration (although you can modify state of the objects retrieved from pool) * in case there is no recycled objects left in the pool, NULL is stored. * in case pool is not seeded, NULL is stored. pool pool obj target variable ************************************************************************/ public static void ae_shared_pool_next_recycled(shared_pool pool, ref T obj) where T : alglib.apobject { /* exit on end of list */ if( pool.enumeration_counter==null ) { obj = null; return; } /* modify internal enumeration counter */ pool.enumeration_counter = pool.enumeration_counter.next_entry; /* exit on empty list */ if( pool.enumeration_counter==null ) { obj = null; return; } /* assign object to smart pointer */ obj = (T)pool.enumeration_counter.obj; } /************************************************************************ This function clears internal list of recycled objects and seed object. However, pool still can be used (after initialization with another seed). pool pool state ALGLIB environment state NOTE: this function is NOT thread-safe. It does not acquire pool lock, so you should NOT call it when lock can be used by another thread. ************************************************************************/ public static void ae_shared_pool_reset(shared_pool pool) { pool.seed_object = null; pool.recycled_objects = null; pool.enumeration_counter = null; } } } public partial class alglib { public partial class smp { public static int cores_count = 1; } public class smpselftests { public static bool runtests() { return true; } } }