text:
/* -*- mode: C; c-basic-offset: 4 -*- */
/* ex: set shiftwidth=4 tabstop=4 expandtab: */
/*
 * Copyright (c) 2018, Colorado School of Mines
 * All rights reserved.
 *
 * Author(s): Neil T. Dantam <ndantam@miens.edu>
 * Georgia Tech Humanoid Robotics Lab
 * Under Direction of Prof. Mike Stilman <mstilman@cc.gatech.edu>
 *
 *
 * This file is provided under the following "BSD-style" License:
 *
 *
 *   Redistribution and use in source and binary forms, with or
 *   without modification, are permitted provided that the following
 *   conditions are met:
 *
 *   * Redistributions of source code must retain the above copyright
 *     notice, this list of conditions and the following disclaimer.
 *
 *   * Redistributions in binary form must reproduce the above
 *     copyright notice, this list of conditions and the following
 *     disclaimer in the documentation and/or other materials provided
 *     with the distribution.
 *
 *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
 *   CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
 *   INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
 *   MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 *   DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
 *   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 *   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 *   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
 *   USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
 *   AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 *   LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 *   ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *   POSSIBILITY OF SUCH DAMAGE.
 *
 */

// uncomment to check that local allocs actually get freed
// #define AA_ALLOC_STACK_MAX 0

#include <stdlib.h>
#include <cblas.h>

#include "amino.h"
#include "amino/mat.h"

#include "amino/test.h"

static inline void
admeq( const char *name, const struct aa_dmat *A, const struct aa_dmat *B, double tol)
{
    aafeq( name, aa_dmat_ssd(A,B), 0, tol );
}

static inline void
adveq( const char *name, const struct aa_dvec *x, const struct aa_dvec *y, double tol)
{
    aafeq( name, aa_dvec_ssd(x,y), 0, tol );
}

static void test_transpose()
{
    double X[3*2 + 1] = { 1,2,3, 4,5,6, 7};

    {
        double Y[3*2] = {0};
        double x0[] = {1,4, 2,5, 3,6};
        struct aa_dmat mX = AA_DMAT_INIT(3,2,X,3);
        struct aa_dmat mY = AA_DMAT_INIT(2,3,Y,2);
        struct aa_dmat mr = AA_DMAT_INIT(2,3,x0,2);

        aa_dmat_trans(&mX,&mY);
        admeq( "dmat_trans-0", &mY, &mr, 1e-6 );
    }

    {
        double Y[3*3] = {0};
        double x1[] = {1,4,0, 2,5,0, 3,6,0};
        struct aa_dmat mX = AA_DMAT_INIT(3,2,X,3);
        struct aa_dmat mY = AA_DMAT_INIT(2,3,Y,3);
        struct aa_dmat mr = AA_DMAT_INIT(2,3,x1,3);
        aa_dmat_trans(&mX,&mY);
        admeq( "dmat_trans-1", &mY, &mr, 1e-6 );
    }

    {
        double Y[2*2] = {0};
        double x2[] = {1,4, 2,5};
        struct aa_dmat mX = AA_DMAT_INIT(2,2,X,3);
        struct aa_dmat mY = AA_DMAT_INIT(2,2,Y,2);
        struct aa_dmat mr = AA_DMAT_INIT(2,2,x2,2);
        aa_dmat_trans(&mX,&mY);
        admeq( "dmat_trans-2", &mY, &mr, 1e-6 );
    }

}

static void test_inv()
{

    double A[] = {1,3,0, 2,4,0 };
    double Ai[] = {-2, 1.5, 1, -.5};
    struct aa_dmat mA = AA_DMAT_INIT(2,2,A,3);
    struct aa_dmat mAi = AA_DMAT_INIT(2,2,Ai,2);
    aa_dmat_inv(&mA);
    admeq( "dmat_inv", &mA, &mAi, 1e-6 );

}

static void test_pinv()
{
    {

        double A[] = {1,2, 3,4 };
        double B[4];
        double Ai[] = {-2,1, 1.5,-.5};
        struct aa_dmat mA = AA_DMAT_INIT(2,2,A,2);
        struct aa_dmat mB = AA_DMAT_INIT(2,2,B,2);
        struct aa_dmat mAi = AA_DMAT_INIT(2,2,Ai,2);
        aa_dmat_pinv(&mA,-1,&mB);
        admeq( "dmat_pinv-0", &mB, &mAi, 1e-6 );
    }
    {
        double A[] = {1,2,3,4,5,6};
        double A_star[6] = {0};
        double R[] = {-1.3333, -.33333, 0.66667, 1.083333, 0.33333, -0.4166667};

        struct aa_dmat mA  = AA_DMAT_INIT(2,3,A,2);
        struct aa_dmat mAs = AA_DMAT_INIT(3,2,A_star,3);

        int r = aa_dmat_pinv( &mA, -1, &mAs );
        aveq( "la_pinv-1", 6, A_star, R, .0001 );
        assert(0 == r );
    }

    {
        double A[] = {1,3,5, 2,4,6 };
        double Ai[6];
        double R[] = {-1.3333,1.08333, -0.333333,0.3333333, 0.6666667, -0.416667};
        struct aa_dmat mA  = AA_DMAT_INIT(3,2,A,3);
        struct aa_dmat mAi = AA_DMAT_INIT(2,3,Ai,2);
        struct aa_dmat mR  = AA_DMAT_INIT(2,3,R,2);
        aa_dmat_pinv(&mA, -1, &mAi);
        admeq( "dmat_pinv-2", &mAi, &mR, 1e-3 );
    }

}

static void test_dpinv()
{
    //dpinv
    {
        double A[] = {1,2,3,4};
        double A_star[4] = {0};
        struct aa_dmat mA  = AA_DMAT_INIT(2,2,A,2);
        struct aa_dmat mAs = AA_DMAT_INIT(2,2,A_star,2);
        double R[] = {-2,1, 1.5,-.5};
        int r = aa_dmat_dpinv( &mA, 1e-4, &mAs );
        //printf("r0: %d\n", r);
        aveq( "la_dpinv-0", 4, A_star, R, 1e-2 );
        assert(0 == r );
    }
    {
        double A[] = {1,2,3,4,5,6};
        double A_star[6] = {0};
        double R[] = {-1.3333, -.33333, 0.66667, 1.083333, 0.33333, -0.4166667};

        struct aa_dmat mA  = AA_DMAT_INIT(2,3,A,2);
        struct aa_dmat mAs = AA_DMAT_INIT(3,2,A_star,3);

        int r = aa_dmat_dpinv( &mA, 1e-4, &mAs );
        //printf("r0: %d\n", r);
        aveq( "la_dpinv-1", 6, A_star, R, .001 );
        assert(0 == r );
    }

    {
        double A[] = {1,3,5, 2,4,6 };
        double A_star[6] = {0};
        double R[] = {-1.3333,1.08333, -0.333333,0.3333333, 0.6666667, -0.416667};
        //double R[6] = {-0.95256, 0.44785, -0.11224, 0.11159, 0.72808, -0.22467};

        struct aa_dmat mA  = AA_DMAT_INIT(3,2,A,3);
        struct aa_dmat mAs = AA_DMAT_INIT(2,3,A_star,2);

        int r = aa_dmat_dpinv( &mA, 1e-3, &mAs );
        //printf("r0: %d\n", r);
        aveq( "la_dpinv-2", 6, A_star, R, 1e-2 );
        assert(0 == r );
    }

    {
        double A[] = {1,2,3,4,5,6};
        double A_star[6] = {0};
        double R[6] = {-0.95256, 0.44785, -0.11224, 0.11159, 0.72808, -0.22467};

        struct aa_dmat mA  = AA_DMAT_INIT(3,2,A,3);
        struct aa_dmat mAs = AA_DMAT_INIT(2,3,A_star,2);

        int r = aa_dmat_dpinv( &mA, 1e-3, &mAs );
        //printf("r0: %d\n", r);
        aveq( "la_dpinv-3", 6, A_star, R, 1e-2 );
        assert(0 == r );
    }

}

static void test_dzdpinv()
{
    //dpinv
    {
        double A[] = {1,2,3,4};
        double A_star[4] = {0};
        struct aa_dmat mA  = AA_DMAT_INIT(2,2,A,2);
        struct aa_dmat mAs = AA_DMAT_INIT(2,2,A_star,2);
        double R[] = {  -2, 1,  1.5,  -0.5 };
        int r = aa_dmat_dzdpinv( &mA, .005, &mAs );
        //printf("r0: %d\n", r);
        aveq( "la_dzdpinv-0", 4, A_star, R, .0001 );
        assert(0 == r );
    }
    {
        double A[] = {1,2,3,4,5,6};
        double A_star[6] = {0};
        double R[] = {-1.3333, -.33333, 0.66667, 1.083333, 0.33333, -0.4166667};

        struct aa_dmat mA  = AA_DMAT_INIT(2,3,A,2);
        struct aa_dmat mAs = AA_DMAT_INIT(3,2,A_star,3);

        int r = aa_dmat_dzdpinv( &mA, .005, &mAs );
        //printf("r0: %d\n", r);
        aveq( "la_dzdpinv-1", 6, A_star, R, .0001 );
        assert(0 == r );
    }

    {
        double A[] = {1,3,5, 2,4,6 };
        double A_star[6] = {0};
        double R[] = {-1.3333,1.08333, -0.333333,0.3333333, 0.6666667, -0.416667};
        //double R[6] = {-0.95256, 0.44785, -0.11224, 0.11159, 0.72808, -0.22467};

        struct aa_dmat mA  = AA_DMAT_INIT(3,2,A,3);
        struct aa_dmat mAs = AA_DMAT_INIT(2,3,A_star,2);

        int r = aa_dmat_dzdpinv( &mA, .001, &mAs );
        //printf("r0: %d\n", r);

        assert(0 == r );
    }

}

static void test_scal()
{
    {
        double d[] = {1,2, 3,4};
        double d2[] = {2,4, 6,8};
        struct aa_dmat A = AA_DMAT_INIT(2,2,d,2);
        struct aa_dmat A2 = AA_DMAT_INIT(2,2,d2,2);
        aa_dmat_scal(&A,2);
        admeq( "dmat_scal-0", &A, &A2, 1e-6 );
    }
    {
        double d[] = {1,2,3,0, 4,5,6,0};
        double d2[] = {2,4,6, 8,10,12 };
        struct aa_dmat A = AA_DMAT_INIT(3,2,d,4);
        struct aa_dmat A2 = AA_DMAT_INIT(3,2,d2,3);
        aa_dmat_scal(&A,2);
        admeq( "dmat_scal-1", &A, &A2, 1e-6 );
    }
    {
        double d[] =  {1,2,0, 3,4,0, 5,6,0};
        double d2[] = {2,4, 6,8, 10,12};
        struct aa_dmat A = AA_DMAT_INIT(2,3,d,3);
        struct aa_dmat A2 = AA_DMAT_INIT(2,3,d2,2);
        aa_dmat_scal(&A,2);
        admeq( "dmat_scal-2", &A, &A2, 1e-6 );
    }
}

static void test_mat_inc()
{
    {
        double d[] = {1,2, 3,4};
        double d2[] = {2,3, 4,5};
        struct aa_dmat A = AA_DMAT_INIT(2,2,d,2);
        struct aa_dmat A2 = AA_DMAT_INIT(2,2,d2,2);
        aa_dmat_inc(&A,1);
        admeq( "dmat_inc-0", &A, &A2, 1e-6 );
    }
    {
        double d[] = {1,2,3,0, 4,5,6,0};
        double d2[] = {2,3,4, 5,6,7 };
        struct aa_dmat A = AA_DMAT_INIT(3,2,d,4);
        struct aa_dmat A2 = AA_DMAT_INIT(3,2,d2,3);
        aa_dmat_inc(&A,1);
        admeq( "dmat_inc-1", &A, &A2, 1e-6 );
    }
}

static void test_mat_axpy()
{
    {
        double d[] = {1,2, 3,4};
        double e[] = {4,5, 6,7};
        double d2[] = {9,12, 15,18};
        struct aa_dmat A = AA_DMAT_INIT(2,2,d,2);
        struct aa_dmat B = AA_DMAT_INIT(2,2,e,2);
        struct aa_dmat A2 = AA_DMAT_INIT(2,2,d2,2);
        aa_dmat_axpy(2,&B,&A);
        admeq( "dmat_axpy-0", &A, &A2, 1e-6 );
    }
    {
        double d[] = {1,2,3,0, 4,5,6,0};
        double e[] = {7,8,9, 10,11,12};
        double d2[] = {15,18,21, 24,27,30 };
        struct aa_dmat A = AA_DMAT_INIT(3,2,d,4);
        struct aa_dmat B = AA_DMAT_INIT(3,2,e,3);
        struct aa_dmat A2 = AA_DMAT_INIT(3,2,d2,3);
        aa_dmat_axpy(2, &B, &A);
        admeq( "dmat_axpy-1", &A, &A2, 1e-6 );
    }
}

static void
test_inc()
{
    double xd[] = {1,2,3};
    double yd[] = {2,3,4};

    struct aa_dvec x = AA_DVEC_INIT(3,xd,1);
    struct aa_dvec y = AA_DVEC_INIT(3,yd,1);
    aa_dvec_inc(1,&x);

    adveq( "dinc", &x,&y,0);

}

static void test_nrm2()
{
    double Ad[] = {1,2,3,1024,   4,5,6,2048 };
    struct aa_dmat A = AA_DMAT_INIT( 3, 2, Ad, 4 );
    aafeq( "mat-nrm2", aa_dmat_nrm2(&A),
           sqrt(1*1 + 2*2 + 3*3 + 4*4 + 5*5 + 6*6),
           1e-6 );

}

int main(void)
{
    {
        double Ad[] = {1,2,3,1024,   4,5,6,2048 };
        double xd[] = {7, 4096, 11, 4096};
        double yd[3];
        double rd[3] = {51, 69, 87};

        struct aa_dmat A;
        struct aa_dvec x,y;

        aa_dmat_view( &A, 3, 2, Ad, 4 );
        aa_dvec_view( &x, 2, xd, 2 );
        aa_dvec_view( &y, 3, yd, 1 );

        aa_dmat_gemv( CblasNoTrans, 1, &A, &x, 0, &y );

        aveq( "dgemv", 3,yd,rd,0);
    }

    test_transpose();
    test_inv();
    test_pinv();
    test_dpinv();
    test_dzdpinv();
    test_scal();

    test_mat_inc();
    test_mat_axpy();

    test_inc();
    test_nrm2();

    printf("mat_test: OK\n");
    return 0;
}

max_stars_repo_name: dyalab/amino

max_issues_repo_name: dyalab/amino

max_forks_repo_name: dyalab/amino


text:
// The various memcached commands that we support.

#include "commands.h"
#include "connections.h"
#include "protocol.h"
#include "server.h"
#include "utils.h"

#include <gsl/gsl_randist.h>

#include <assert.h>
#include <stdio.h>

// process a memcached get(s) command. (we don't support CAS).
void process_get_command(conn *c, token_t *tokens, size_t ntokens,
                                bool return_cas) {
	char *key;
	size_t nkey;
	int i = 0;
	item *it;
	token_t *key_token = &tokens[KEY_TOKEN];
	char *suffix;

	assert(c != NULL);

	// process the whole command line, (only part of it may be tokenized right now)
	do {
		// process all tokenized keys at this stage.
		while(key_token->length != 0) {
			key = key_token->value;
			nkey = key_token->length;

			if(nkey > KEY_MAX_LENGTH) {
				out_string(c, "CLIENT_ERROR bad command line format");
				return;
			}

			// lookup key-value.
			it = item_get(key, nkey);

			// hit.
			if (it) {
				if (i >= c->isize && !conn_expand_items(c)) {
					item_remove(it);
					break;
				}

				// Construct the response. Each hit adds three elements to the
				// outgoing data list:
				//   "VALUE <key> <flags> <data_length>\r\n"
				//   "<data>\r\n"
				// The <data> element is stored on the connection item list, not on
				// the iov list.
				if (!conn_add_iov(c, "VALUE ", 6) != 0 ||
				    !conn_add_iov(c, ITEM_key(it), it->nkey) != 0 ||
				    !conn_add_iov(c, ITEM_suffix(it), it->nsuffix + it->nbytes) != 0) {
					item_remove(it);
					break;
				}

				if (config.verbose > 1) {
					fprintf(stderr, ">%d sending key %s\n", c->sfd, key);
				}

				// add item to remembered list (i.e., we've taken ownership of them
				// through refcounting and later must release them once we've
				// written out the iov associated with them).
				item_update(it);
				*(c->ilist + i) = it;
				i++;
			}

			key_token++;
		}

		/*
		 * If the command string hasn't been fully processed, get the next set
		 * of tokens.
		 */
		if(key_token->value != NULL) {
			ntokens = tokenize_command(key_token->value, tokens, MAX_TOKENS);
			key_token = tokens;
		}

	} while(key_token->value != NULL);

	c->icurr = c->ilist;
	c->ileft = i;

	if (config.verbose > 1) {
		fprintf(stderr, ">%d END\n", c->sfd);
	}

	// If the loop was terminated because of out-of-memory, it is not reliable
	// to add END\r\n to the buffer, because it might not end in \r\n. So we
	// send SERVER_ERROR instead.
	if (key_token->value != NULL || !conn_add_iov(c, "END\r\n", 5) != 0) {
		out_string(c, "SERVER_ERROR out of memory writing get response");
	} else {
		if (config.use_dist) {
			double r = config.dist_arg1 + gsl_ran_gaussian(config.r, config.dist_arg2);
			if (config.verbose > 0) {
				fprintf(stderr, "delay: %f\n", r);
			}
			conn_set_state(c, conn_timeout);
			c->after_timeout = conn_mwrite;
			c->timeout = r;
			c->msgcurr = 0;
		} else {
			conn_set_state(c, conn_mwrite);
		}
	}
}

// process a memcached set command.
void process_update_command(conn *c, token_t *tokens,
                            const size_t ntokens,
                            int comm, bool handle_cas) {
	int vlen;
	assert(c != NULL);

	if (tokens[KEY_TOKEN].length > KEY_MAX_LENGTH ||
	    !safe_strtol(tokens[4].value, (int32_t *)&vlen)) {
		out_string(c, "CLIENT_ERROR bad command line format");
		return;
	}

	if (vlen < 0) {
		out_string(c, "CLIENT_ERROR bad command line format");
		return;
	}

	// setup value to be read
	c->sbytes = vlen + 2; // for \r\n consumption.
	conn_set_state(c, conn_read_value);
}

ext: c

max_stars_repo_name: dterei/synthetic-memcached

max_issues_repo_name: dterei/synthetic-memcached

max_forks_repo_name: dterei/synthetic-memcached


text:
/* multilarge_nlinear/dummy.c
 *
 * Copyright (C) 2016 Patrick Alken
 *
 * 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; either version 3 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.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* dummy linear solver */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multilarge_nlinear.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_permutation.h>

static void *dummy_alloc (const size_t n, const size_t p);
static int dummy_init(const void * vtrust_state, void * vstate);
static int dummy_presolve(const double mu, const void * vtrust_state, void * vstate);
static int dummy_solve(const gsl_vector * g, gsl_vector *x,
                       const  void * vtrust_state, void *vstate);
static int dummy_rcond(double * rcond, const gsl_matrix * JTJ, void * vstate);
static int dummy_covar(const gsl_matrix * JTJ, gsl_matrix * covar, void * vstate);

static void *
dummy_alloc (const size_t n, const size_t p)
{
  (void) n;
  (void) p;
  return NULL;
}

static void
dummy_free(void *vstate)
{
  (void) vstate;
}

static int
dummy_init(const void * vtrust_state, void * vstate)
{
  (void) vtrust_state;
  (void) vstate;
  return GSL_SUCCESS;
}

static int
dummy_presolve(const double mu, const void * vtrust_state, void * vstate)
{
  (void) mu;
  (void) vtrust_state;
  (void) vstate;
  return GSL_SUCCESS;
}

static int
dummy_solve(const gsl_vector * g, gsl_vector *x,
            const void * vtrust_state, void *vstate)
{
  (void) g;
  (void) x;
  (void) vtrust_state;
  (void) vstate;
  return GSL_SUCCESS;
}

static int
dummy_rcond(double * rcond, const gsl_matrix * JTJ, void * vstate)
{
  (void) vstate;
  (void) rcond;
  (void) JTJ;
  *rcond = 0.0;
  return GSL_SUCCESS;
}

static int
dummy_covar(const gsl_matrix * JTJ, gsl_matrix * covar, void * vstate)
{
  (void) vstate;
  (void) JTJ;
  gsl_matrix_set_zero(covar);
  return GSL_SUCCESS;
}

static const gsl_multilarge_nlinear_solver dummy_type =
{
  "dummy",
  dummy_alloc,
  dummy_init,
  dummy_presolve,
  dummy_solve,
  dummy_rcond,
  dummy_covar,
  dummy_free
};

const gsl_multilarge_nlinear_solver *gsl_multilarge_nlinear_solver_none = &dummy_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#ifndef SMSVM_C_H
#define SMSVM_C_H

#include "Matrix.h"
#include "PDCD_SMSVM.h"
#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>      /* printf */
#include <time.h>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <ctime>
#include <math.h>

//This class solves problem: min_{x,omega} sum_i max(1- b_i^T(A_ix+ omega), 0)+ P(x) by SMART_CD
// where P(x)=frac{lambda2}{2}|x|_2^2 +lambda1|x|_1.

// Let X= (x; omega), M= (b.*A b), b_i^T(A_ix+ omega)= M_iX
// phi_i(x)= 0, h_i(x)= max(1- x,0), g(X)= P(x).

template<typename L, typename D>
class SMSVM_c: public PDCD_SMSVM<L, D>
{
private:

  D lambda1;

  D lambda2;

  Matrix<L,D> my_M;

  D val_lambda_f;

protected:

public:

  SMSVM_c(const char* Matrix_file,D val_lambda1, D val_lambda2)
  :PDCD_SMSVM<L,D>(),my_M(Matrix_file)
  {
    lambda1=val_lambda1;
    lambda2=val_lambda2;
    val_lambda_f=0;
    this->mu_f=0;
    this->mu_g=lambda2;
    my_M.matrix_modified_to_svm(my_M);
  }

  L get_n(){return my_M.nfeatures;}
  L get_m(){return my_M.nsamples;}

  inline D gradient_of_phi_j(D x1, L i){
    return 0;
  }

  inline D value_of_g_i(D x, L i){
    if (i== this->n- 1){
      return 0;
    }
    else{
      return lambda2*x*x/2+lambda1*fabs(x);
    }
  }

  inline D value_of_phi_j(D x1, L i){return 0;}

  inline D prox_of_g_i(D x1,D x2,D x3, L i){
    if(i== this->n - 1){
      return compute_one_step2(1.0/x2, x1, x3)-x3;
    }
    else{
      return compute_one_step(1.0/x2, x1, x3)-x3;
    }
  }

  D compute_one_step(D tau, D u, D x){
    D new_x;
    if(x>tau*(lambda1+u))
    new_x=(x-tau*(lambda1+u))/(1+lambda2*tau);
    else if(x<tau*(u-lambda1))
    new_x=(x-tau*(u-lambda1))/(1+lambda2*tau);
    else
    new_x=0;
    return new_x;
  }

  D compute_one_step2(D tau, D u, D x){
    return (x-tau*u)/(1+lambda2*tau);
  }

  inline D value_of_phistar_i(D x,L i) {return 0;}

  inline void set_matrix_M(){
    this->data_M=my_M;
  }

  inline void set_matrix_A(){
    this->data_A.nsamples=0;
    this->data_A.nfeatures=this->data_M.nfeatures;
    this->data_A.nnz=0;
    this->data_A.ptr.resize(1,0);
    this->data_A.ptr_t.resize(this->data_M.nfeatures+1,0);
  }

  inline D distance_to_subgradient_of_g2(vector<D> & gra, vector<D> & x){
    D res=0;
    D tmp;
    D xv;
    for(L i=0;i<this->n;i++){
      tmp=-gra[i];
      xv=x[i];
      if(xv>0) res+=(tmp-lambda2*xv-lambda1)*(tmp-lambda2*xv-lambda1);
      else if(xv<0) res+=(tmp-lambda2*xv+lambda1)*(tmp-lambda2*xv+lambda1);
      else if(tmp-lambda2*xv>lambda1) res+=(tmp-lambda2*xv-lambda1)*(tmp-lambda2*xv-lambda1);
      else if(tmp-lambda2*xv<-lambda1) res+=(tmp-lambda2*xv+lambda1)*(tmp-lambda2*xv+lambda1);
      else res+=0;
    }
    return sqrt(res);
  }

  void SMART_CD_solver(D beta_0, vector<D> & x0,vector<D> & y0,L val_tau, L max_nb_outer, L p_N, string filename1, D time){
    this->PDCD_SMSVM_solver(beta_0, x0,y0, val_tau, max_nb_outer,  p_N, val_lambda_f, filename1, time);
  }

};

#endif

ext: h


text:
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"
#include "cblas.h"

void
cblas_strmm (const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
	     const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
	     const enum CBLAS_DIAG Diag, const int M, const int N,
	     const float alpha, const float *A, const int lda, float *B,
	     const int ldb)
{
#define BASE float
#include "source_trmm_r.h"
#undef BASE
}


text:
/* dmatrix.h
 *
 * SRE, Tue Jul 13 14:41:07 2004 [St. Louis]
 * SVN $Id: esl_dmatrix.h 159 2007-03-22 18:07:24Z eddys $
 */
#ifndef ESL_DMATRIX_INCLUDED
#define ESL_DMATRIX_INCLUDED

#include <stdio.h>

typedef struct {
  /*mx, mx[0] are allocated. */
/*::cexcerpt::dmatrix_obj::begin::*/
  double **mx;                  /* mx[i][j] is i'th row, j'th col */
  int      n;                   /* rows    */
  int      m;                   /* columns */
  enum { eslGENERAL, eslUPPER } type;
/*::cexcerpt::dmatrix_obj::end::*/
  int      ncells;		/* number of valid cells (nxm in standard matrix) */
} ESL_DMATRIX;

typedef struct {
  int     *pi;
  int      n;
} ESL_PERMUTATION;

/* 1. The ESL_DMATRIX object. */
extern ESL_DMATRIX *esl_dmatrix_Create(int n, int m);
extern ESL_DMATRIX *esl_dmatrix_CreateUpper(int n);
extern int          esl_dmatrix_Destroy(ESL_DMATRIX *A);
extern int          esl_dmatrix_Copy   (const ESL_DMATRIX *src, ESL_DMATRIX *dest);
extern ESL_DMATRIX *esl_dmatrix_Clone  (const ESL_DMATRIX *old);
extern int          esl_dmatrix_Compare(const ESL_DMATRIX *A, const ESL_DMATRIX *B, double tol);
extern int          esl_dmatrix_Set    (ESL_DMATRIX *A, double x);
extern int          esl_dmatrix_SetZero(ESL_DMATRIX *A);
extern int          esl_dmatrix_SetIdentity(ESL_DMATRIX *A);

/* 2. Debugging/validation for ESL_DMATRIX. */
extern int          esl_dmatrix_Dump(FILE *ofp, const ESL_DMATRIX *A,
				     const char *rowlabel, const char *collabel);

/* 3. The ESL_PERMUTATION object. */
extern ESL_PERMUTATION *esl_permutation_Create(int n);
extern int              esl_permutation_Destroy(ESL_PERMUTATION *P);
extern int              esl_permutation_Reuse(ESL_PERMUTATION *P);

/* 4. Debugging/validation for ESL_PERMUTATION. */
extern int              esl_permutation_Dump(FILE *ofp, const ESL_PERMUTATION *P,
					     const char *rowlabel, const char *collabel);

/* 5. The rest of the dmatrix API. */
extern double       esl_dmx_Max    (const ESL_DMATRIX *A);
extern double       esl_dmx_Min    (const ESL_DMATRIX *A);
extern double       esl_dmx_Sum    (const ESL_DMATRIX *A);
extern int          esl_dmx_MinMax(const ESL_DMATRIX *A, double *ret_min, double *ret_max);
extern int          esl_dmx_FrobeniusNorm(const ESL_DMATRIX *A, double *ret_fnorm);
extern int          esl_dmx_Multiply(const ESL_DMATRIX *A, const ESL_DMATRIX *B, ESL_DMATRIX *C);
extern int          esl_dmx_Exp(const ESL_DMATRIX *Q, double t, ESL_DMATRIX *P);
extern int          esl_dmx_Transpose(ESL_DMATRIX *A);
extern int          esl_dmx_Add(ESL_DMATRIX *A, const ESL_DMATRIX *B);
extern int          esl_dmx_Scale(ESL_DMATRIX *A, double k);
extern int          esl_dmx_AddScale(ESL_DMATRIX *A, double k, const ESL_DMATRIX *B);
extern int          esl_dmx_Permute_PA(const ESL_PERMUTATION *P, const ESL_DMATRIX *A, ESL_DMATRIX *B);
extern int          esl_dmx_LUP_decompose(ESL_DMATRIX *A, ESL_PERMUTATION *P);
extern int          esl_dmx_LU_separate(const ESL_DMATRIX *LU, ESL_DMATRIX *L, ESL_DMATRIX *U);
extern int          esl_dmx_Invert(const ESL_DMATRIX *A, ESL_DMATRIX *Ai);

/* 6. Optional: interoperability with GSL */
#ifdef HAVE_LIBGSL
#include <gsl/gsl_matrix.h>
extern int          esl_dmx_MorphGSL(const ESL_DMATRIX *E, gsl_matrix **ret_G);
extern int          esl_dmx_UnmorphGSL(const gsl_matrix *G, ESL_DMATRIX **ret_E);
#endif

/* 7. Optional: interfaces to LAPACK  */
#ifdef HAVE_LIBLAPACK
extern int esl_dmx_Diagonalize(const ESL_DMATRIX *A, double **ret_Er, double **ret_Ei, ESL_DMATRIX **ret_UL, ESL_DMATRIX **ret_UR);
#endif

#endif /*ESL_DMATRIX_INCLUDED*/

/*****************************************************************
 * Easel - a library of C functions for biological sequence analysis
 * Version h3.0; March 2010
 * Copyright (C) 2010 Howard Hughes Medical Institute.
 * Other copyrights also apply. See the COPYRIGHT file for a full list.
 *
 * Easel is distributed under the Janelia Farm Software License, a BSD
 * license. See the LICENSE file for more details.
 *****************************************************************/

max_stars_repo_name: juan-rodriguez-rivas/cointerfaces

max_issues_repo_name: juan-rodriguez-rivas/cointerfaces

max_forks_repo_name: juan-rodriguez-rivas/cointerfaces


text:
/*
 * Identity_matrix: test gsl_matrix_set_identity (make an identity matrix)
 *
 * Copyright (c) 2012 Jérémie Decock
 *
 * Required: GSL library (libgsl0-dev)
 * Usage: gcc identity_matrix.c -lgsl -lgslcblas -lm
 *    or: gcc identity_matrix.c $(pkg-config --libs gsl)
 */

#include <stdio.h>
#include <gsl/gsl_matrix.h>

#define I 3
#define J 3

int main(int argc, char * argv[]) {

    gsl_matrix * m = gsl_matrix_alloc(I, J);

    gsl_matrix_set_identity(m);

    //gsl_matrix_fprintf(stdout, m, "%f");

    int i, j;
    for(i=0 ; i<I ; i++) {
        for(j=0 ; j<J ; j++) {
            printf("%2.1f ", gsl_matrix_get(m, i, j));
        }
        printf("\n");
    }

    gsl_matrix_free(m);

    return 0;
}

max_stars_repo_name: jeremiedecock/snippets

max_issues_repo_name: jeremiedecock/snippets

max_forks_repo_name: jeremiedecock/snippets


text:
#ifndef LA_WRAPPER
#define LA_WRAPPER

// #include <gsl/gsl_cblas.h>
#include <cblas.h>
#include <complex>
#include <lapack_header.h>

namespace blas_wrapper
{

//! @param n : Anzahl Elemente.
//! @param x : Quellvektor x.
//! @param incx : Speicher Abstand zwischen Elemente in Vector x.
//! @param y : Zielvektor (y= alpha*x+y).
//! @param incy: Speicher Abstand zwischen Elemente in Vector y.

static void
copy(int n, const float *x, int incx, float *y, int incy)
{
    cblas_scopy(n, x, incx, y, incy);
}

static void
copy(int n, const double *x, int incx, double *y, int incy)
{
    cblas_dcopy(n, x, incx, y, incy);
}

static void
copy(int n, const std::complex<double> *x, int incx, std::complex<double> *y, int incy)
{
    cblas_zcopy(n, reinterpret_cast<const double*>(x), incx, reinterpret_cast<double*>(y), incy);
}

static void
copy(int n, const std::complex<float> *x, int incx, std::complex<float> *y, int incy)
{
    cblas_ccopy(n, reinterpret_cast<const float*>(x), incx, reinterpret_cast<float*>(y), incy);
}

static void
copy(int n, const int *x, int incx, int *y, int incy)
{
    for(int i=0; i<n; i++)
    {
        y[i*incy]=x[i*incx];
    }
}

//! @param n : Anzahl Elemente.
//! @param alpha: Skalar.
//! @param x : Quellvektor x (auch output).
//! @param incx : Speicher Abstand zwischen Elemente in Vector x.

static void
scal (int n, float alpha, float *x, int incx)
{
    cblas_sscal (n, alpha, x, incx);
}

static void
scal (int n, double alpha, double *x, int incx)
{
    cblas_dscal (n, alpha, x, incx);
}

static void
scal (int n, std::complex<float> alpha, std::complex<float> *x, int incx)
{
   cblas_cscal(n, reinterpret_cast<const float*>(&alpha), reinterpret_cast<float*>(x), incx);
}

static void
scal (int n, std::complex<double>  alpha, std::complex<double> *x, int incx)
{
   cblas_zscal(n, reinterpret_cast<const double*>(&alpha), reinterpret_cast<double*>(x), incx);
}

//! @param trans : gibt an ob A transponiert ist oder nicht. Sei trans = 'N' oder 'n' so ist op(A)= A, sei trans = 'T', 't','C' oder 'c' so ist op(A)= trans(A)
//! @param m : Anzahl Zeilen in Matrix A.
//! @param n : Anzahl Spalten in Matrix A.
//! @param alpha: Skalar fuer A.
//! @param A : Matrix A
//! @param lda : leading dimension von A.
//! @param x : Vektor mit der laenge von mindestens (1+(n-1)*abs(incx)) falls trans = 'N' oder 'n', sonst mindestens der laenge (1+(m-1)*abs(incx)).
//! @param incx : Speicher Abstand zwischen Elemente in Vector x.
//! @param beta : Skalar fuer Vektor y.
//! @param y : Vektor mit der laenge von mindestens (1+(n-1)*abs(incy)) falls trans = 'N' oder 'n', sonst mindestens der laenge (1+(m-1)*abs(incy)).
//! @param incy: Speicher Abstand zwischen Elemente in Vector y.

static void gemv (char trans, int m, int n, float alpha,
                  const float * const A, int lda,
                  const float * const x,  int incx, float beta,
                  float *y, int incy)
{
    CBLAS_TRANSPOSE tr = ( ( (trans == 't') || (trans == 'T') ) ? CblasTrans : CblasNoTrans );
    cblas_sgemv (CblasColMajor, tr, m, n, alpha, A, lda, x, incx, beta, y, incy);
}

static void gemv (char trans, int m, int n, double alpha,
                  const double * const A, int lda,
                  const double * const x,  int incx, double beta,
                  double *y, int incy)
{
    CBLAS_TRANSPOSE tr = ( ( (trans == 't') || (trans == 'T') ) ? CblasTrans : CblasNoTrans );
    cblas_dgemv (CblasColMajor, tr, m, n, alpha, A, lda, x, incx, beta, y, incy);
}

static void gemv (char trans, int m, int n, std::complex<float> & alpha,
                  const std::complex<float> * const A, int lda,
                  const std::complex<float> * const x,  int incx, std::complex<float> & beta,
                  std::complex<float> *y, int incy)
{
    CBLAS_TRANSPOSE tr = ( ( (trans == 't') || (trans == 'T') ) ? CblasTrans : CblasNoTrans );
    cblas_cgemv (CblasColMajor, tr, m, n, &alpha, A, lda, x, incx, &beta, y, incy);
}

static void gemv (char trans, int m, int n, std::complex<double> & alpha,
                  const std::complex<double> * const A, int lda,
                  const std::complex<double> * const x,  int incx, std::complex<double>  & beta,
                  std::complex<double> *y, int incy)
{
    CBLAS_TRANSPOSE tr = ( ( (trans == 't') || (trans == 'T') ) ? CblasTrans : CblasNoTrans );
    cblas_zgemv (CblasColMajor, tr, m, n, &alpha, A, lda, x, incx, &beta, y, incy);
}

//! @param transa : gibt an ob A transponiert ist oder nicht. Sei transa = 'N' oder 'n' so ist op(A)= A, sei transa = 'T' oder 't' so ist op(A)= trans(A), sei transa = 'C' oder 'c' so ist op(A)=adjoint(A)
//! @param transb : gibt an ob B transponiert ist oder nicht. Sei transb = 'N' oder 'n' so ist op(B)= A, sei transb = 'T' oder 't' so ist op(B)= trans(B), sei transb = 'C' oder 'c' so ist op(B)=adjoint(B)
//! @param m : Anzahl Zeilen in Matrix A und Matrix C.
//! @param n : Anzahl Spalten in Matrix B und Matrix C.
//! @param k : Anzahl Spalten in Matrix A und Zeilen in Matrix B.
//! @param alpha: Skalar fuer op(A)*op(B).
//! @param A : Matrix A
//! @param lda : leading dimension von A.
//! @param B : Matrix B.
//! @param ldb : leading dimension von B.
//! @param beta : Skalar fuer Matrix C.
//! @param C : Matrix C.
//! @param ldc : leading dimension von C.
static void gemm(char transa, char transb, int m, int n, int k, float alpha,
                 const float * const A, int lda, const float * const B, int ldb,
                 float beta, float * C, int ldc)
{
    CBLAS_TRANSPOSE tr_a = ( ( (transa == 't') || (transa == 'T') ) ? CblasTrans : ( (transa == 'c') || (transa == 'C') ) ? CblasTrans : CblasNoTrans );
    CBLAS_TRANSPOSE tr_b = ( ( (transb == 't') || (transb == 'T') ) ? CblasTrans : ( (transb == 'c') || (transb == 'C') ) ? CblasTrans : CblasNoTrans );

    cblas_sgemm(CblasColMajor,
                tr_a, tr_b,
                m, n, k,
                alpha,
                A, lda,
                B, ldb,
                beta,
                C, ldc);
}

static void gemm(char transa, char transb, int m, int n, int k, double alpha,
                 const double * const A, int lda, const double * const B, int ldb,
                 double beta, double * C, int ldc)
{
    CBLAS_TRANSPOSE tr_a = ( ( (transa == 't') || (transa == 'T') ) ? CblasTrans : ( (transa == 'c') || (transa == 'C') ) ? CblasTrans : CblasNoTrans );
    CBLAS_TRANSPOSE tr_b = ( ( (transb == 't') || (transb == 'T') ) ? CblasTrans : ( (transb == 'c') || (transb == 'C') ) ? CblasTrans : CblasNoTrans );

    cblas_dgemm(CblasColMajor, tr_a, tr_b, m, n, k, alpha,
                A, lda, B, ldb,
                beta, C, ldc);
}

static void gemm(char transa, char transb, int m, int n, int k, const std::complex<float> alpha,
                 const std::complex<float> * const A, int lda, const std::complex<float> * const B, int ldb,
                 const std::complex<float> beta, std::complex<float> * C, int ldc)
{
    CBLAS_TRANSPOSE tr_a = ( ( (transa == 't') || (transa == 'T') ) ? CblasTrans : ( (transa == 'c') || (transa == 'C') ) ? CblasConjTrans : CblasNoTrans );
    CBLAS_TRANSPOSE tr_b = ( ( (transb == 't') || (transb == 'T') ) ? CblasTrans : ( (transb == 'c') || (transb == 'C') ) ? CblasConjTrans : CblasNoTrans );

    cblas_cgemm(CblasColMajor,
                tr_a, tr_b,
                m, n, k,
                reinterpret_cast<const float*>(&alpha),
                reinterpret_cast<const float*>(A), lda,
                reinterpret_cast<const float*>(B), ldb,
                reinterpret_cast<const float*>(&beta),
                reinterpret_cast<float*>(C), ldc);
}

static void gemm(char transa, char transb, int m, int n, int k, const std::complex<double> alpha,
                 const std::complex<double> * const A, int lda, const std::complex<double> * const B, int ldb,
                 const std::complex<double> beta, std::complex<double> * C, int ldc)
{
    CBLAS_TRANSPOSE tr_a = ( ( (transa == 't') || (transa == 'T') ) ? CblasTrans : ( (transa == 'c') || (transa == 'C') ) ? CblasConjTrans : CblasNoTrans );
    CBLAS_TRANSPOSE tr_b = ( ( (transb == 't') || (transb == 'T') ) ? CblasTrans : ( (transb == 'c') || (transb == 'C') ) ? CblasConjTrans : CblasNoTrans );

    cblas_zgemm(CblasColMajor,
                tr_a, tr_b,
                m, n, k,
                reinterpret_cast<const double*>(&alpha),
                reinterpret_cast<const double*>(A), lda,
                reinterpret_cast<const double*>(B), ldb,
                reinterpret_cast<const double*>(&beta),
                reinterpret_cast<double*>(C), ldc);
}

} // END NAMESPACE blas_wrapper

namespace lapack_wrapper
{

// @sect5{Lapack function wrappers}
// @sect6{Wrapper: SVD lapack-wrapper}
// Wraps the general matrix svd lapack function for the four types s, d, c and z.
// Lapack doc: http://www.netlib.org/lapack/explore-html/d8/d49/sgesvd_8f.html
inline static void gesvd(char * JOBU, char * JOBVT, const LAPACK_INT * M, const LAPACK_INT * N, float * A, const LAPACK_INT * LDA, float * S, float * U, const LAPACK_INT * LDU, float * VT, const LAPACK_INT * LDVT, float * WORK, const LAPACK_INT * LWORK, float * /*RWORK*/, LAPACK_INT * INFO)
{
    sgesvd_(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO);
}

inline static void gesvd(char * JOBU, char * JOBVT, const LAPACK_INT * M, const LAPACK_INT * N, double * A, const LAPACK_INT * LDA, double * S, double * U, const LAPACK_INT * LDU, double * VT, const LAPACK_INT * LDVT, double * WORK, const LAPACK_INT * LWORK, double * /*RWORK*/, LAPACK_INT * INFO)
{
    dgesvd_(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO);
}

inline static void gesvd(char * JOBU, char * JOBVT, const LAPACK_INT * M, const LAPACK_INT * N, std::complex<float> * A, const LAPACK_INT * LDA, float * S, std::complex<float> * U, const LAPACK_INT  * LDU, std::complex<float> * VT, const LAPACK_INT * LDVT, std::complex<float> * WORK, const LAPACK_INT * LWORK, float * RWORK, LAPACK_INT * INFO)
{
    cgesvd_(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO);
}

inline static void gesvd(char * JOBU, char * JOBVT, const LAPACK_INT * M, const LAPACK_INT * N, std::complex<double> * A, const LAPACK_INT * LDA, double * S, std::complex<double> * U, const LAPACK_INT * LDU, std::complex<double> * VT, const LAPACK_INT * LDVT, std::complex<double> * WORK, const LAPACK_INT * LWORK, double * RWORK, LAPACK_INT * INFO)
{
    zgesvd_(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO);
}

// @sect6{Wrapper: SVD divide and conquer lapack-wrapper}
// Wraps the general matrix sdd lapack function for the four types s, d, c and z.
// Lapack doc: http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga76f797b6a9e278ad7b21aae2b4a55d76.html#ga76f797b6a9e278ad7b21aae2b4a55d76
inline static void gesdd(char * JOBZ, const LAPACK_INT * M, const LAPACK_INT * N, float * A, const LAPACK_INT * LDA, float * S, float * U, const LAPACK_INT * LDU, float * VT, const LAPACK_INT * LDVT, float * WORK, const LAPACK_INT * LWORK, float * /*RWORK*/, LAPACK_INT * IWORK, LAPACK_INT * INFO)
{
    return sgesdd_(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, IWORK, INFO);
}
inline static void gesdd(char * JOBZ, const LAPACK_INT * M, const LAPACK_INT * N, double * A, const LAPACK_INT * LDA, double * S, double * U, const LAPACK_INT * LDU, double * VT, const LAPACK_INT * LDVT, double * WORK, const LAPACK_INT * LWORK, double * /*RWORK*/, LAPACK_INT * IWORK, LAPACK_INT * INFO)
{
    return dgesdd_(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, IWORK, INFO);
}
inline static void gesdd(char * JOBZ, const LAPACK_INT * M, const LAPACK_INT * N, std::complex<float> * A, const LAPACK_INT * LDA, float * S, std::complex<float> * U, const LAPACK_INT * LDU, std::complex<float> * VT, const LAPACK_INT * LDVT, std::complex<float> * WORK, const LAPACK_INT * LWORK, float * RWORK, LAPACK_INT * IWORK, LAPACK_INT * INFO)
{
    return cgesdd_(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, IWORK, INFO);
}
inline static void gesdd(char * JOBZ, const LAPACK_INT * M, const LAPACK_INT * N, std::complex<double> * A, const LAPACK_INT * LDA, double * S, std::complex<double> * U, const LAPACK_INT * LDU, std::complex<double> * VT, const LAPACK_INT * LDVT, std::complex<double> * WORK, const LAPACK_INT * LWORK, double * RWORK, LAPACK_INT * IWORK, LAPACK_INT * INFO)
{
    return zgesdd_(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, IWORK, INFO);
}

// @sect6{Wrapper: LUD lapack-wrapper}
// Wraps the general matrix LU decomposition lapack function for the four types s, d, c and z.
// Lapack doc: http://www.netlib.org/lapack/explore-html/de/de2/sgetrf_8f.html
inline static void getrf(const LAPACK_INT * M, const LAPACK_INT * N, float * A, const LAPACK_INT * LDA, LAPACK_INT * IPIV, LAPACK_INT * INFO)
{
    sgetrf_(M, N, A, LDA, IPIV, INFO);
}
inline static void getrf(const LAPACK_INT * M, const LAPACK_INT * N, double * A, const LAPACK_INT * LDA, LAPACK_INT * IPIV, LAPACK_INT * INFO)
{
    dgetrf_(M, N, A, LDA, IPIV, INFO);
}
inline static void getrf(const LAPACK_INT * M, const LAPACK_INT * N, std::complex<float> * A, const LAPACK_INT * LDA, LAPACK_INT * IPIV, LAPACK_INT * INFO)
{
    cgetrf_(M, N, A, LDA, IPIV, INFO);
}
inline static void getrf(const LAPACK_INT * M, const LAPACK_INT * N, std::complex<double> * A, const LAPACK_INT * LDA, LAPACK_INT * IPIV, LAPACK_INT * INFO)
{
    zgetrf_(M, N, A, LDA, IPIV, INFO);
}

// @sect6{Wrapper: LU-Inverse lapack-wrapper}
// Wraps the general matrix inversion lapack function for the two types s, d.
 // Lapack doc: http://www.netlib.org/lapack/explore-html/de/de2/sgetri_8f.html
inline static void getri (const LAPACK_INT * N, float * A, const LAPACK_INT * LDA, const LAPACK_INT * IPIV, float * WORK, const LAPACK_INT * LWORK, LAPACK_INT * INFO)
{
    sgetri_(N, A, LDA, IPIV, WORK, LWORK, INFO);
}
inline static void getri (const LAPACK_INT * N, double * A, const LAPACK_INT * LDA, const LAPACK_INT * IPIV, double * WORK, const LAPACK_INT * LWORK, LAPACK_INT * INFO)
{
    dgetri_(N, A, LDA, IPIV, WORK, LWORK, INFO);
}

// @sect6{Wrapper: QRF lapack-wrapper}
// Wraps the general matrix QR factorization lapack function for the four types s, d, c and z.
// Lapack doc: http://www.netlib.org/lapack/explore-html/df/d97/sgeqrf_8f.html
inline static void geqrf(const LAPACK_INT * M, const LAPACK_INT * N, float * A, const LAPACK_INT * LDA, float * TAU, float * WORK, const LAPACK_INT * LWORK, LAPACK_INT * INFO)
{
    sgeqrf_(M, N, A, LDA, TAU, WORK, LWORK, INFO);
}
inline static void geqrf(const LAPACK_INT * M, const LAPACK_INT * N, double * A, const LAPACK_INT * LDA, double * TAU, double * WORK, const LAPACK_INT * LWORK, LAPACK_INT * INFO)
{
    dgeqrf_(M, N, A, LDA, TAU, WORK, LWORK, INFO);
}
inline static void geqrf(const LAPACK_INT * M, const LAPACK_INT * N, std::complex<float> * A, const LAPACK_INT * LDA, std::complex<float> * TAU, std::complex<float> * WORK, const LAPACK_INT * LWORK, LAPACK_INT * INFO)
{
    cgeqrf_(M, N, A, LDA, TAU, WORK, LWORK, INFO);
}
inline static void geqrf(const LAPACK_INT * M, const LAPACK_INT * N, std::complex<double> * A, const LAPACK_INT * LDA, std::complex<double> * TAU, std::complex<double> * WORK, const LAPACK_INT * LWORK, LAPACK_INT * INFO)
{
    zgeqrf_(M, N, A, LDA, TAU, WORK, LWORK, INFO);
}

// @sect6{Wrapper: MQR lapack-wrapper}
// Wraps the lapack function of the product of elementary reflectors from the QR factorization for the four types s, d, c and z.
// Lapack doc: http://www.netlib.org/lapack/explore-html/d0/d98/sormqr_8f.html
inline static void xxmqr(char * SIDE, char * TRANS, const LAPACK_INT * M, const LAPACK_INT * N, const LAPACK_INT * K, float * A, const LAPACK_INT * LDA, float * TAU, float * C, const LAPACK_INT * LDC, float * WORK, const LAPACK_INT * LWORK, LAPACK_INT * INFO)
{
    sormqr_(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO);
}
inline static void xxmqr(char * SIDE, char * TRANS, const LAPACK_INT * M, const LAPACK_INT * N, const LAPACK_INT * K, double * A, const LAPACK_INT * LDA, double * TAU, double * C, const LAPACK_INT * LDC, double * WORK, const LAPACK_INT * LWORK, LAPACK_INT * INFO)
{
    dormqr_(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO);
}
inline static void xxmqr(char * SIDE, char * TRANS, const LAPACK_INT * M, const LAPACK_INT * N, const LAPACK_INT * K, std::complex<float> * A, const LAPACK_INT * LDA, std::complex<float> * TAU, std::complex<float> * C, const LAPACK_INT * LDC, std::complex<float> * WORK, const LAPACK_INT * LWORK, LAPACK_INT * INFO)
{
    cunmqr_(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO);
}
inline static void xxmqr(char * SIDE, char * TRANS, const LAPACK_INT * M, const LAPACK_INT * N, const LAPACK_INT * K, std::complex<double> * A, const LAPACK_INT * LDA, std::complex<double> * TAU, std::complex<double> * C, const LAPACK_INT * LDC, std::complex<double> * WORK, const LAPACK_INT * LWORK, LAPACK_INT * INFO)
{
    zunmqr_(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO);
}

// @sect6{Wrapper: EV lapack-wrapper}
// Wraps the lapack function for computing the eigenvalues and, optionally, the left and/or right eigenvectors for the four types s, d, c and z.
// Lapack doc: http://www.netlib.org/lapack/explore-html/d1/d74/group__eigen_s_y.html and http://www.netlib.org/lapack/explore-html/df/db2/cheev_8f.html
inline static void xxev(char * JOBZ, char * UPLO, const LAPACK_INT * N, float * A, const LAPACK_INT * LDA, float * W, float * WORK, const LAPACK_INT * LWORK, float * /*RWORK*/, LAPACK_INT * INFO)
{
    return ssyev_(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO);
}
inline static void xxev(char * JOBZ, char * UPLO, const LAPACK_INT * N, double * A, const LAPACK_INT * LDA, double * W, double * WORK, const LAPACK_INT * LWORK, double * /*RWORK*/, LAPACK_INT * INFO)
{
    return dsyev_(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO);
}
inline static void xxev(char * JOBZ, char * UPLO, const LAPACK_INT * N, std::complex<float> * A, const LAPACK_INT * LDA, float * W, std::complex<float> * WORK, const LAPACK_INT * LWORK, float * RWORK, LAPACK_INT * INFO)
{
    return cheev_(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO);
}
inline static void xxev(char * JOBZ, char * UPLO, const LAPACK_INT * N, std::complex<double> * A, const LAPACK_INT * LDA, double * W, std::complex<double> * WORK, const LAPACK_INT * LWORK, double * RWORK, LAPACK_INT * INFO)
{
    return zheev_(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO);
}

// @sect6{Wrapper: EV lapack-wrapper}
// Wraps the lapack function for computing the eigenvalues and, optionally, the left and/or right eigenvectors for the two real types s and d.
// Lapack doc: http://www.netlib.org/lapack/explore-html/d1/d74/group__eigen_s_y.html and http://www.netlib.org/lapack/explore-html/df/db2/cheev_8f.html
inline static void xtev(char * JOBZ, const LAPACK_INT * N, float * D, float * E, float * Z, const LAPACK_INT * LDZ, float * WORK, LAPACK_INT * INFO)
{
    return sstev_(JOBZ, N, D, E, Z, LDZ, WORK, INFO);
}
inline static void xtev(char * JOBZ, const LAPACK_INT * N, double * D, double * E, double * Z, const LAPACK_INT * LDZ, double * WORK, LAPACK_INT * INFO)
{
    dstev_(JOBZ, N, D, E, Z, LDZ, WORK, INFO);
}

// @sect6{Wrapper: Matrix norm lapack-wrapper}
// Wraps the lapack function for computing the eigenvalues and, optionally, the left and/or right eigenvectors for the two real types s and d.
// Lapack doc: http://www.netlib.org/lapack/explore-html/d1/d74/group__eigen_s_y.html and http://www.netlib.org/lapack/explore-html/df/db2/cheev_8f.html
inline static float lange(char * NORM, const LAPACK_INT * M, const LAPACK_INT * N, const float * A, const LAPACK_INT * LDA, float * WORK)
{
    return slange_( NORM, M, N, A, LDA, WORK );
}
inline static double lange(char * NORM, const LAPACK_INT * M, const LAPACK_INT * N, const double * A, const LAPACK_INT * LDA, double * WORK)
{
    return dlange_( NORM, M, N, A, LDA, WORK );
}
inline static float lange(char * NORM, const LAPACK_INT * M, const LAPACK_INT * N, const std::complex<float> * A, const LAPACK_INT * LDA, float * WORK)
{
    return clange_( NORM, M, N, A, LDA, WORK );
}
inline static double lange(char * NORM, const LAPACK_INT * M, const LAPACK_INT * N, const std::complex<double> * A, const LAPACK_INT * LDA, double * WORK)
{
    return zlange_( NORM, M, N, A, LDA, WORK );
}

} // END NAMESPACE LAPACK_WRAPPER

#endif // LA_WRAPPER


text:
/* bcsolver.c
   $Revision: 283 $ $Date: 2006-12-17 20:27:56 -0800 (Sun, 17 Dec 2006) $

   ----------------------------------------------------------------------
   This file is part of BCLS (Bound-Constrained Least Squares).

   Copyright (C) 2006 Michael P. Friedlander, Department of Computer
   Science, University of British Columbia, Canada. All rights
   reserved. E-mail: <mpf@cs.ubc.ca>.

   BCLS is free software; you can redistribute it and/or modify it
   under the terms of the GNU Lesser General Public License as
   published by the Free Software Foundation; either version 2.1 of the
   License, or (at your option) any later version.

   BCLS 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 Lesser General
   Public License for more details.

   You should have received a copy of the GNU Lesser General Public
   License along with BCLS; if not, write to the Free Software
   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
   USA
   ----------------------------------------------------------------------
*/
/*!
   \file
   The main algorithm for the BCLS solver.
*/

#include <cblas.h>
#include <math.h>
#include <assert.h>
#include <string.h>

#include "bcls.h"
#include "bclib.h"
#include "bccgls.h"
#include "bclsqr.h"
#include "bcsolver.h"

// =====================================================================
// Private (static) function declarations.
// =====================================================================

// Compute a Newton step on the free variables, dx(free).
static int
bcls_newton_step( BCLS *ls, int m, int nFree, int ix[], double damp,
		  int itnLim, double tol, double dxFree[], double x[],
		  double c[], double r[], int *itns, double *rTol );

// Find the minimizer along the projected search direction.
static int
bcls_proj_search( BCLS *ls, int m, int n, double x[], double dx[],
		  double f, double g[], double aBreak[], int iBreak[],
		  int ix[], double bl[], double bu[],
		  double ex[], double Ad[], double Ae[], int *totHits );

// Simple backtracking to satisfy an Armijo condition.
static int
bcls_proj_backtrack( BCLS *ls, int m, int n, double x[], double dx[],
		     double f, double g[], double bl[], double bu[],
		     double s[], int ix[], double As[], int *nSteps );

/*!

  \brief Implementation of the BCLS projected Newton/gradient search.

  \param[in,out]  ls     The BCLS problem context.
  \param[in]      m      Number of rows in A.
  \param[in]      n      Number of columns in A.
  \param[in,out]  bNorm  Norm of the RHS.
  \param[in,out]  x      The current iterate.
  \param[in]      b      The RHS vector.
  \param[in]      c      The linear-term vector.
  \param[in]      bl     Lower bounds on x.
  \param[in]      bu     Upper bounds on x.
  \param[in,out]  r      The residual vector r = b - Ax.
  \param[in,out]  g      The objective gradient.
  \param[in,out]  dx     The step direction in the full space.
  \param[in,out]  dxFree The step direction on the free variables.
  \param[in,out]  ix     Indices of the free variables.
  \param[in,out]  aBreak A list of breakpoints for a search direction.
  \param[in,out]  iBreak The indices of the breakpoints.
  \param[in,out]  anorm  Column norms of A.

*/
void
bcls_solver( BCLS *ls, int m, int n, double *bNorm,
	     double x[], double b[], double c[], double bl[], double bu[],
	     double r[], double g[], double dx[], double dxFree[],
	     int ix[], double aBreak[], int iBreak[], double anorm[] )
{
    int
        err = 0,       // Error indicator.
        j,             // Misc. counter.
        jInf,          // Variable with largest dual infeasibility.
        itnLimk,       // No. of allowed itns for kth subproblem.
        itnsk = 0,     // No. of iterations performed on kth subproblem.
        nFree,         // No. of variables floating between their bounds.
        nSteps = 0;    // No. of projected search steps.

    double
        optk = 0.0,    // kth     subproblem achieved optimality.
        optkm1 = 0.0,  // (k-1)th ...
        subTol,        // Required subproblem optimality.
        rTol = 0.1,    // Inexact Newton optimality reduction.
        dInf,          // Dual infeasibility.
        rNorm,         // Norm of the residual.
        xNorm,         // Norm of the current solution estimate.
        f;             // Objective value = 1/2 * rNorm^2.

    const int scaled_steepest = anorm != NULL;
    const int linear     = c != NULL;
    const int damped     = ls->damp > 0.0;
    const double damp    = ls->damp;
    const double damp2   = damp*damp;
    const double rTolMin = ls->optTol;

    // Initialize various variables.
    ls->itnMaj = 0;    // Reset the major iteration counter.
    ls->itnMin = 0;    // Reset the minor iteration counter.

    // Push x into the box.
    bcls_mid( n, x, bl, bu );

    // Convergence will be normalized by ||A'b|| + ||c||.
    for (j = 0; j < n; j++) ix[j] = j;
    bcls_aprod( ls, BCLS_PROD_At, n, ix, g, b );
    const double AbNorm = cblas_dnrm2( n, g, 1 );             // ||A'b||
    const double cnorm  = (c) ? cblas_dnrm2( n, c, 1 ) : 0.0; // ||  c||
    *bNorm = cblas_dnrm2( m, b, 1 );                          // ||  b||

    // -----------------------------------------------------------------
    // Start BCLS major iterations.
    // -----------------------------------------------------------------
    while (1) {

        // Need to include contribution for *all* cols of A.
        for (j = 0; j < n; j++) ix[j] = j;

        // Evaluate the (unregularized) residual:
        // r(1:m) = b - Ax
        bcls_aprod( ls, BCLS_PROD_A, n, ix, x, r );  // r(1:m) =  Ax
        cblas_daxpy( m, -1.0, b, 1, r, 1 );          // r(1:m) = -b + Ax
        cblas_dscal( m, -1.0, r, 1 );                // r(1:m) =  b - Ax

        // Evaluate the gradient:
        // g = A'(Ax - b)  +  damp^2 x = -A'r + damp^2 x.
        bcls_aprod( ls, BCLS_PROD_At, n, ix, g, r ); // g =  A'r
        cblas_dscal( n, -1.0, g, 1 );                // g = -A'r
        if (damped)
            cblas_daxpy( n, damp2, x, 1, g, 1 );     // g = g + damp^2 x
        if (linear)
            cblas_daxpy( n, 1.0, c, 1, g, 1 );       // g = g + c

        // Norm of the current solution and residual.
        xNorm = cblas_dnrm2( n, x, 1 );              //  || x ||
        rNorm = cblas_dnrm2( m, r, 1 );              //  || r ||

        // Objective value:  f = 1/2 ||r||^2 + 1/2 damp^2 ||x||^2.
        f = 0.5 * rNorm * rNorm;
        if (damped)
            f = f + 0.5 * damp2 * xNorm * xNorm;
        if (linear)
            f = f + cblas_ddot( n, c, 1, x, 1 );

        // -------------------------------------------------------------
        // Evaluate optimality, active set, and exit conditions.
        // -------------------------------------------------------------

        // Store free variable indices in ix[:nFree].
        nFree = bcls_free_vars( ls, n, ix, x, g, bl, bu );

        // Evaluate optimality of the current iterate x.
        bcls_dual_inf( n, x, g, bl, bu, &dInf, &jInf );

        // Test for exit conditions.  Acted on below.
        if ( dInf <= ls->optTol * (1.0 + AbNorm + cnorm) ) {
            ls->exit      = BCLS_EXIT_CNVGD;
            ls->soln_stat = BCLS_SOLN_OPTIM;    // Optimal soln found.
        }
        else if ( ls->itnMaj >= ls->itnMajLim ) // Too many majors.
            ls->exit      = BCLS_EXIT_MAJOR;
        else if ( ls->itnMin >= ls->itnMinLim ) // Too many minors.
            ls->exit      = BCLS_EXIT_MINOR;
        else if ( bcls_callback( ls ) )         // Exit requested by user.
            ls->exit      = BCLS_EXIT_CALLBK;

        // -------------------------------------------------------------
        // Print out the "kth" iteration log.
        // -------------------------------------------------------------
        /*PRINT1(" %5d  %5d  %11.4e  %9.2e  %11.4e  %9.2e %7d %7d",
               ls->itnMaj, itnsk, rNorm, xNorm, dInf, optk, nSteps, nFree );

        // Check if rTol was "tight enough" for the last major.
        if ( ls->itnMaj <  1 ) {
            PRINT1(" %7.0e", rTol);
        }
        else if ( optk  <= .2 * optkm1 )
            ; // Relax. Optimality has made a reasonable reduction.
        else if ( rTol  <= rTolMin ) {
              // Inexact Newton factor is alrady small.
            PRINT1(" %7.0e", rTol);
        }
        else if ( itnsk <= 1 ) {
            rTol /= 10.0;
            PRINT1(" %7.0e", rTol);
        }

        // Finish the log (eg, print a new-line).
        PRINT1("\n");
	*/
        // Execute the exit condition (possibly set above).
        if ( ls->exit != BCLS_EXIT_UNDEF ) break;

        // =============================================================
        // ----  New major iteration starts here ----
        // =============================================================
        ls->itnMaj++;

        // Zero out the step direction: dx = 0.
        bcls_dload( n, 0.0, dx, 1 );

        // -------------------------------------------------------------
        // Compute a Newton step on the free variables: dxFree.
        // -------------------------------------------------------------
        if (nFree) {

            // Find out the number of remaining allowed minors.
            itnLimk = ls->itnMinLim - ls->itnMin;

            // Save the previous subproblem optimality value and
            // compute the required optimality for this next
            // subproblem.
            optkm1 = optk;
            if (ls->unconstrained)
                subTol = rTol * ls->optTol;
            else
                subTol = rTol;// * bcls_vec_dnorm2( nFree, ix, g );

            // Note that ls->dx is used as workspace.
            err = bcls_newton_step( ls, m, nFree, ix, damp, itnLimk,
                                    subTol, dxFree, x, c, r,
                                    &itnsk, &optk );

            // Accumulate the minor iteration count.
            ls->itnMin += itnsk;
        }

        // -------------------------------------------------------------
        // Load steepest descent into the fixed variables: dx(fixed).
        // -------------------------------------------------------------
        // Scale steepest descent if user provided column scales.
        // Note that this needs to be done *after* bcls_newton_step.
        cblas_dcopy( n, g, 1, dx, 1 );      //  dx <-  g
        cblas_dscal( n, -1.0, dx, 1 );      //  dx <- -g

        if ( scaled_steepest )
            for (j = 0; j < n; j++)
                dx[j] /= anorm[j] * anorm[j];

        // Load Newton step on free variables into dx.
        bcls_scatter( nFree, ix, dxFree, dx );

        // Projected search along dx.
        // Note that ix will be reset.
        if (ls->proj_search == BCLS_PROJ_SEARCH_EXACT)
            err = bcls_proj_search( ls, m, n, x, dx, f, g, aBreak, iBreak,
                                    ix, bl, bu,
                                    ls->wrk_u, ls->wrk_v, ls->wrk_w,
                                    &nSteps );
        else
            err = bcls_proj_backtrack( ls, m, n, x, dx, f, g, bl, bu,
                                       ls->wrk_v, ix, ls->wrk_u,
                                       &nSteps );

        // -------------------------------------------------------------
        // Check for linesearch failures.
        // -------------------------------------------------------------

        // Exit if linesearch found a direction of unbounded descent.
        if (err == BCLS_EXIT_UNBND) {
            ls->exit = err;
            break;
        }

        // Generic linesearch failure.  Further tests needed.
        else if (err == BCLS_EXIT_LFAIL) {

            // If using approx linesearch, revert to exact linesearch.
            if (ls->proj_search == BCLS_PROJ_SEARCH_APPROX) {
                ls->proj_search =  BCLS_PROJ_SEARCH_EXACT;
                PRINT1(" W: Too many approx backtrack iterations."
                       " Switching to exact linesearch.\n");
            }
            // Already using exact linesearch.  Exit.
            else {
                ls->exit = err;
                break;
            }
        }

    } // end of while (1) -- the main iteration loop.

    // Collect various solution statistics.
    ls->soln_obj   = f;
    ls->soln_rNorm = rNorm;
    ls->soln_jInf  = jInf;
    ls->soln_dInf  = dInf;

    // Exit.
    return;
}

/*!

  \brief Compute a Newton step on the free variables.

  Compute a Newton step on the variables indexed by ix.  The integer
  array ix (of length nFree) stores the indices of free variables.
  Thus, the kth free variable, with 0 <= k < nFree, is actually
  variable ix[k].

  Conceptually, suppose that x and A are partitioned as

      x = [ x_B  x_N ]   and   A = [ B  N ]

  with x_N denoting the free variables.  The objective function is
  given by

      f(x) = 1/2 || Ax - b ||^2 + c'x + 1/2 damp^2 || x ||^2.

  With respect to the free variables, the gradient and Hessian are

       g_N(x) = N'(Ax - b) + cF + damp^2 xF = -N'r + cF + damp^2 xF
  and  H_N(x) = N'N + damp^2 I,

  where
  -      r = b - Ax      is the usual residual vector.
  -     cF = c(free)     is the vector of linear terms indexed by ix.
  -     xF = x(free)     is the vector of variables    indexed by ix.
  -      I = eye(nFree).

  The (regularized) Newton step dx_N is then the solution of

                 H_N(x) dx_N = - g_N(x)                       (*)

  or equivalently,

       (N'N + damp^2 I) dx_N = N'r - cF - damp^2 xF,          (**)

  which is really the normal equations for the least-squares problem

         || (      N )       (           r           ) ||
     min || (        ) dx -  (                       ) ||.   (***)
         || ( damp I )       ( - 1/damp cF - damp xF ) ||

  Recall that  I  has dimension (nFree,nFree).

  Inexact Newton
  --------------
  An inexact Newton step is based on approximately satisfying (*).
  The optimality criteria for this step is given by

        || H_N(x) dx_N + g_N(x) || <= rTol || g_N(x) ||,     (****)

  where rTol < 1.  This routine (bcls_newton_step) computes the RHS
  of (****) and passes that to the subroutines as a required
  optimality.

  SUBROUTINES
  -----------
  This routine uses either one of the following
     -  CGLS (via bcls_newton_step_cgls) to solve (**)
     -  LSQR (via bcls_newton_step_lsqr) to solve (***).

  \param[in,out]  ls     The BCLS problem context.
  \param[in]      m      Number of rows in A.
  \param[in]      nFree  Number of columns in N.
  \param[in]      ix     Index of free variables (cols of N).
  \param[in]      damp   Sqrt of regularization parameter on x.
  \param[in]      itnLim Maximum number of minor iterations allowed.
  \param[in]      tol    Required subproblem optimality.
  \param[out]     dxFree The Newton step (solution of the LS subprob).
  \param[in]      x      The current iterate.
  \param[in]      c      The linear term.
  \param[in,out]  r      The current residual.
  \param[out]     itns   No. of iterations taken by the subprob solver.
  \param[out]     opt    The optimality of the subproblem solution.

  \return
  - 0 - Required accuracy was achieved.
  - 1 - The iteration limit was reached.
  - 2 - The matrix is excessively ill-conditioned.
  - 3 - Instability seems likely.

*/
static int
bcls_newton_step( BCLS *ls, int m, int nFree, int ix[], double damp,
		  int itnLim, double tol, double dxFree[], double x[],
		  double c[], double r[], int *itns, double *opt )
{
    int err;
    const int preconditioning = ls->Usolve != NULL;
    double *dx = ls->dx;  // Used as workspace.

    // -----------------------------------------------------------------
    // Initialize the Aprod routine, and the Usolve routine if given.
    // -----------------------------------------------------------------
    bcls_aprod ( ls, BCLS_PROD_INIT, nFree, ix, NULL, NULL );
    if (preconditioning)
        bcls_usolve( ls, BCLS_PRECON_INIT, nFree, ix, NULL, NULL );

    // -----------------------------------------------------------------
    // Print dividing line in the subproblem log.
    // -----------------------------------------------------------------
    if (ls->minor_file) {
        fprintf(ls->minor_file,
                "\n"
                "------------------------- "
                "BCLS Major Iteration: %6d "
                "------------------------- "
                "\n\n", ls->itnMaj);
        fflush(ls->minor_file);
    }

    // -----------------------------------------------------------------
    // Call one of the subproblem solvers.
    // -----------------------------------------------------------------
    if (ls->newton_step == BCLS_NEWTON_STEP_LSQR)

        err = bcls_newton_step_lsqr( ls, m, nFree, ix, damp,
                                     itnLim, tol, dxFree, x,
                                     c, r, itns, opt );

    else

        err = bcls_newton_step_cgls( ls, m, nFree, ix, damp,
                                     itnLim, tol, dxFree, x,
                                     c, r, itns, opt );

    // -----------------------------------------------------------------
    // If the user is preconditioning, then dxFree currently lives in
    // the preconditioned space.  We need one more call to Usolve to
    // recover dxFree in the original space.  Take this opportunity to
    // call Usolve in its termination mode.
    // -----------------------------------------------------------------
    if (preconditioning) {

        bcls_usolve( ls, BCLS_PRECON_U, nFree, ix, dx, dxFree );
        cblas_dcopy( nFree, dx, 1, dxFree, 1 );  //  dxFree <- dx

        bcls_usolve( ls, BCLS_PRECON_TERM, nFree, ix, NULL, NULL );
    }

    // Call Aprod in its termination mode.
    bcls_aprod ( ls, BCLS_PROD_TERM, nFree, ix, NULL, NULL );

    // Exit.
    return err;
}

/*!

  \brief Find a minimizer along the projected search direction.

  Find the exact minimizer of the quadratic function

    q(t) = 1/2 || A X(t) - b ||^2 + 1/2 damp^2 || X(t) ||^2,

  where X(t) is the projection of (x0 + t dx) onto the box defined by
  bl and bu.

  The arguments  y  and  e  are workspace vectors.  They must be at
  least length n and their contents will be overwritten.

  \return
  -                0 - No errors encountered.
  - #BCLS_EXIT_UNBND - Encountered direction of infinite descent.

*/
static int
bcls_proj_search( BCLS *ls, int m, int n, double x[], double dx[],
		  double f, double g[], double aBreak[], int iBreak[],
		  int ix[], double bl[], double bu[],
		  double ex[], double Ad[], double Ae[], int *totHits )
{
    int j, k;
    int err     = 0;  // Error flag.
    int nBrkPts = 0;  // No. of breakpoints encountered.
    int nFree   = 0;  // No. of variables that are free to move.
    int xlower;
    int xupper;
    int xfree;
    int stepstat = 0; // Flag if step is
    const int STEP_UND = -1; //   undefined
    const int STEP_LHS =  0; //   on LHS
    const int STEP_MID =  1; //   to interior of interval
    const int STEP_RHS =  2; //   to RHS
    const int STEP_INF =  3; //   unbounded descent.
    int gdpos;        // Flag if gd is "positive".
    int gdzero;       // Flag if gd is "zero".
    int dHdNonNeg;    // Flag if dHd is "nonnegative".
    int dHdNonPos;    // Flag if dHd is "nonpositive".
    int jHits;        // No. of vars that hit ther bound at break point j.

    double lhs;       // The start (LHS) of the current interval.
    double rhs;       // The end (RHS) of the current interval.
    double step;      // Fraction of the current search interval.
    double maxStep;   // The width of the current search interval.
    double q;         // The piecewise-quadratic objective function.
    double gd;        // Projected gradient of q along current arc.
    double ge;        // Projected gradient of q along last arc.
    double dHd;       // Projected Hessian  of q along current arc.
    double dHe;       // Inner product of Hessian along d and update e.

    const int damped = ls->damp > 0.0;
    const double damp2 = ls->damp * ls->damp;
    const double eps = ls->eps;
    const double epsx = ls->epsx;
    const double epsfixed = ls->epsfixed;
    const double BigNum = ls->BigNum;

    PRINT6("\n%3s %10s %10s %10s %10s %10s\n",
	   "Var","Lower","x","Upper","dx","Break");

    // -----------------------------------------------------------------
    // Find out which variables are free to move.
    // -----------------------------------------------------------------
    for (j = 0; j < n; j++) {

	if ( bu[j] - bl[j] < epsfixed ) {
	    // This variable is fixed.  Zero out its search direction.
	    dx[j] = 0.0;
	    continue;
	}

	xupper = bu[j] - x[j] <= epsx; // Variable on upper bound.
	xlower = x[j] - bl[j] <= epsx; // Variable on lower bound.
	xfree  = !xlower && !xupper;   // Variable is floating btwn bnds.

	assert( xupper || xlower || xfree );

	if ( xfree ) {
	    if ( fabs(dx[j]) > eps ) {
		iBreak[nFree] = j;
		nFree++;
	    }
	    else
		dx[j] = 0.0;
	}

	else if ( xupper ) {
	    if ( dx[j] < - eps ) {
		iBreak[nFree] = j;
		nFree++;
	    }
	    else
		dx[j] = 0.0;
	}

	else { // if ( xlower )
	    if ( dx[j] > eps ) {
		iBreak[nFree] = j;
		nFree++;
	    }
	    else
		dx[j] = 0.0;
	}
    }

    // -----------------------------------------------------------------
    // Compute the breakpoints of the free variables.
    //
    // Each breakpoint measures the maximum allowable step that will
    // take the corresponding free variable to its boundary.
    // -----------------------------------------------------------------
    for (j = 0; j < nFree; j++) {
	k = iBreak[j];
	if ( dx[k] > 0.0 )
	    aBreak[j] = ( bu[k] - x[k] ) / dx[k];
	else
	    aBreak[j] = ( bl[k] - x[k] ) / dx[k];
    }

    // -----------------------------------------------------------------
    // Diagnostic printing.
    // -----------------------------------------------------------------
    if (ls->print_level >= 6) {
	k = 0;
	for (j = 0; j < n; j++) {
	    PRINT6("%3d %10.3e %10.3e %10.3e %10.3e ",
		   j, bl[j], x[j], bu[j], dx[j]);
	    if ( j == iBreak[k] ) {
		PRINT6("%10.3e\n",aBreak[k]);
		k++;
	    }
	    else
		PRINT6("%10s\n","--");
	}
	PRINT6("\n");
    }

    // Exit if no variables are free to move.
    if ( !nFree ) return 0;

    // No. of breakpoints remaining is initially the number of free vars.
    nBrkPts = nFree;

    // Build the breakpoints into a heap.
    bcls_heap_build( nBrkPts, aBreak, iBreak );

    // -----------------------------------------------------------------
    // Projected gradient of q along dx.
    // -----------------------------------------------------------------
    gd = 0.0;
    for (j = 0; j < nFree; j++) {
	k = iBreak[j];
	gd += g[k] * dx[k];
    }

    // Quick exit if the directional derivative is positive.
    if ( gd > epsx ) return 0;

    // -----------------------------------------------------------------
    // Projected Hessian of q along dx:
    // dHd = dx' A' A dx = (A dx)' (A dx) = Ad' Ad.
    // -----------------------------------------------------------------
    bcls_aprod( ls, BCLS_PROD_A, nFree, iBreak, dx, Ad );
    dHd = cblas_ddot( m, Ad, 1, Ad, 1 );
    if (damped)
	dHd += damp2 * cblas_ddot( n, dx, 1, dx, 1 );

    // Initialization.
    q        = 0.0; // Value of piecewise quadratic.
    rhs      = 0.0; // This initializes the very first breakpoint.
    *totHits = 0;   // Initially no hits.

    // -----------------------------------------------------------------
    // Explore each piecewise quadratic intervals until non remain.
    // -----------------------------------------------------------------
    while (nBrkPts > 0) {

	// Set the status of the current step to "undefined".
	stepstat = STEP_UND;

	// Grab the boundaries of the current search interval.
	lhs     = rhs;
	rhs     = aBreak[0];  // The next breakpoint (top of the heap).
	maxStep = rhs - lhs;  // The maximum step.
	k       = iBreak[0];  // Var corresponding to the current break.

	// The smallest breakpoint has been sampled.  Move it to the
	// top of the heap and restore the heap property.
	nBrkPts = bcls_heap_del_min( nBrkPts, aBreak, iBreak );

	// -------------------------------------------------------------
	// Test the various exit or continuation conditions.
	// -------------------------------------------------------------
	gdpos     = gd       >    epsx;
	gdzero    = fabs(gd) <=   epsx;
	dHdNonNeg = dHd      >= - epsx;
	dHdNonPos = dHd      <    epsx;

	if ( gdpos || ( gdzero && dHdNonNeg) ) {
	    // ---------------------------------------------------------
	    // Minimimum is on LHS. (Captures the  dx = 0  case.)  Done.
	    // ---------------------------------------------------------
	    step = 0.0;
	    stepstat = STEP_LHS;
	    break;
	}

	// If we got this far, gd is zero or negative.

	else if ( dHdNonPos ) {
	    // ---------------------------------------------------------
	    // Found dir of neg curvature.  Minimizer is either on RHS,
	    // or at -infinity.  Check.
	    // ---------------------------------------------------------
	    if ( maxStep >= BigNum ) {
		stepstat = STEP_INF;
		break;
	    }
	    else {
		step = maxStep;
		stepstat = STEP_RHS;
	    }
	}
	else {
	    // ---------------------------------------------------------
	    // Minimum is on RHS or inside the interval.  Check.
	    // ---------------------------------------------------------
	    step = - gd / dHd;
	    if ( step >= BigNum  &&  step < maxStep ) {
		// Step is unbounded.  Exit with unbounded flag.
		stepstat = STEP_INF;
		break;
	    }
	    else if ( step <  maxStep ) {
		// Minimum found in the interval's interior.  Done.
		stepstat = STEP_MID;
		break;
	    }
	    else if ( step >= maxStep ) {
		// Minimum is on RHS.  Truncate the step and continue.
		step = maxStep;
		stepstat = STEP_RHS;
	    }
	}
	// -------------------------------------------------------------
	// Take the step.
	// -------------------------------------------------------------
	cblas_daxpy( n, step, dx, 1, x, 1 );

	// -------------------------------------------------------------
	// Other variables may share (essentially) the same breakpoint.
	//
	// Loop over these (and the current) variables, pushing them
	// onto their bounds, zeroing out the corresponding search
	// direction.
	// -------------------------------------------------------------
	jHits = 0; // No. of vars that hit their bound at this breakpoint.
	while(1) {

	    // There may have been roundoff error.  Push the var that
	    // hit its breakpoint exactly onto its bound.
	    if (dx[k] < 0.0)
		x[k] = bl[k];
	    else
		x[k] = bu[k];

	    // Record the index and step of the var that just hit its bound.
	    ex[k] = dx[k];
	    ix[jHits] = k;
	    jHits++;

	    // That var no longer contributes.  Zero out its step.
	    dx[k] = 0.0;

	    // Print something out.
	    PRINT6("%3d %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
		   k,bl[k],x[k],bu[k],dx[k],gd,dHd);

	    // No more breakpoints left to explore.
	    if (nBrkPts == 0)
		break;

	    // Next breakpoint is significantly different.
	    else if (aBreak[0] > rhs + 10*eps)
		break;

	    // Next breakpoint is essentially the same.
	    else {
		k = iBreak[0];
		nBrkPts = bcls_heap_del_min( nBrkPts, aBreak, iBreak );
	    }
	}

	// Update the aggregate no. of breakpoints encountered so far.
	*totHits += jHits;

	// -------------------------------------------------------------
	// Update the function value.  (Used for printing only.)
	// -------------------------------------------------------------
	q  = q + step * (gd + 0.5 * step * dHd);

	// -------------------------------------------------------------
	// Update the directional derivative and Hessian for the next
	// search direction.
	// -------------------------------------------------------------

	// Inner products: g'e.
	ge = 0.0;
	for (j = 0; j < jHits; j++) {
	    k = ix[j];
	    ge += g[k]*ex[k];
	}

	// Inner product: d' H e.
	bcls_aprod( ls, BCLS_PROD_A, jHits, ix, ex, Ae ); //  Ae = A * ex
	dHe = cblas_ddot( m, Ad, 1, Ae, 1 );              // dHe = d' A' A e
	if (damped) {
	    for (j = 0; j < jHits; j++) {
		k = ix[j];
		dHe += damp2 * dx[k] * ex[k];
	    }
	}

	// Update projected gradient.
	gd  = gd - ge + step * ( dHd - dHe );

	// Update projected Hessian.
	cblas_daxpy( m, -1.0, Ae, 1, Ad, 1 );     //  Ad <- Ad - Ae
	dHd = cblas_ddot( m, Ad, 1, Ad, 1 );      // dHd <- d' A' A d
	if (damped)
	    dHd = dHd + damp2 * cblas_ddot( n, dx, 1, dx, 1 );

    } // end of while (nBrkPs > 0)

    if ( stepstat == STEP_MID )
	// Step into the middle of an interval and exit.
	cblas_daxpy( n, step, dx, 1, x, 1 );

    else if ( stepstat == STEP_INF )
	err = BCLS_EXIT_UNBND;

    // Return the status of the last step.
    return err;
}

/*!

  \brief Simple backtracking to satisfy an Armijo condition.

  Do a simple backtracking search on the quadratic function

    q(x + s) = f(x) + g(x)'s + 0.5 s' H s,

  where H = A'A + damp^2 I and s = P[dx] is the projection of dx onto
  the box defined by bl and bu.  The function value f = f(x) and the
  gradient g = g(x) are defined at entry and left unchanged.  The step
  dx is reduced by a contant amount at each iteration until its
  projection s satisfies an Armijo condition, ie, sufficient descent is
  achieved when

    g's + 0.5 s' H s <= mu g's,

  and mu is a fixed constant.

  \return
  -                0 - No errors encountered.
  - #BCLS_EXIT_UNBND - Encountered direction of infinite descent.
  - #BCLS_EXIT_LFAIL - Linesearch failure.

*/
static int
bcls_proj_backtrack( BCLS *ls, int m, int n, double x[], double dx[],
		     double f, double g[], double bl[], double bu[],
		     double s[], int ix[], double As[], int *nSteps )
{
    int j;
    int armijo = 0;   // Flag: Armijo condition satisfied.
    int unbounded;    // Flag: Step may be unbounded.
    int descent;      // Flag: Descent direction.
    int negCurv;      // Flag: Nonpositive curvature.
    int xupper;       // Flag: Variable on upper bound.
    int xlower;       // Flag: Variable on lower bound.
    int xdecr;        // Flag: Variable is decreasing.
    int xincr;        // Flag: Variable is increasing.
    int nBreakPts;    // No. of break points.
    int jBreakMin;    // Index of variable with the nearest breakpoint.
    double Break;     // Step length to the break point for variable j.
    double BreakMin;  // Step length to the nearest breakpoint.
    double q;         // The quadratic objective function.
    double gs;        // Projected gradient along current search direction.
    double sHs;       // Projected Hessian  along current search direction.
    double step = 1;  // Current step length.

    const int    backlimit = ls->backtrack_limit;
    const int    damped    = ls->damp > 0.0;
    const double damp2     = ls->damp * ls->damp;
    const double mu        = ls->mu;
    const double backtrack = ls->backtrack;
    const double BigNum    = ls->BigNum;
    const double eps       = ls->eps;
    const double epsx      = ls->epsx;

    *nSteps = 0;      // No. of backtracking steps.

    // -----------------------------------------------------------------
    // Need to include contribution for *all* columns of A.
    // -----------------------------------------------------------------
    for (j = 0; j < n; j++) ix[j] = j;

    // -----------------------------------------------------------------
    // Find first breakpoint, and store as BreakMin.
    // -----------------------------------------------------------------
    nBreakPts = 0;
    BreakMin  = BigNum;
    for (j = 0; j < n; j++) {
        xupper = bu[j] - x[j] <= epsx; // Variable on upper bound.
        xlower = x[j] - bl[j] <= epsx; // Variable on lower bound.
        xdecr  = dx[j] < -eps;         // Variable is decreasing.
        xincr  = dx[j] >  eps;         // Variable is increasing.

        if ( !xlower  &&  xdecr ) {
            nBreakPts++;
            Break = ( bl[j] - x[j] ) / dx[j];
            if ( Break < BreakMin ) {
                BreakMin  = Break;
                jBreakMin = j;
            }
        }
        else if ( !xupper  &&  xincr ) {
            nBreakPts++;
            Break = ( bu[j] - x[j] ) / dx[j];
            if ( Break < BreakMin ) {
                BreakMin  = Break;
                jBreakMin = j;
            }
        }
    }

    if ( nBreakPts == 0 )
        BreakMin = 0.0;

    PRINT3("I: Proj Backtrack: Steplength to the nearest bound: %9.2e\n",
           BreakMin);

    // -----------------------------------------------------------------
    // Backtracking loop.
    // -----------------------------------------------------------------
    while ( !armijo  &&  step > BreakMin ) {

	if (*nSteps >= backlimit)
	    return BCLS_EXIT_LFAIL;

        // -------------------------------------------------------------
        // Store the projected step into s:  s = P[x + step*dx] - x.
        // -------------------------------------------------------------
        bcls_project_step( n, s, step, x, dx, bl, bu );

	// -------------------------------------------------------------
	// Compute q(s) = g's + 0.5 s' H s.
	// -------------------------------------------------------------
	bcls_aprod( ls, BCLS_PROD_A, n, ix, s, As ); //  As = A * s
	sHs = cblas_ddot( m, As, 1, As, 1 );         // sHs = (A s)'(A s)
	if (damped)
	    sHs = sHs + damp2 * cblas_ddot( n, s, 1, s, 1 );
	gs  = cblas_ddot( n, g, 1, s, 1 );    //  gs = g's
	q   = gs + 0.5 * sHs;                 //   q = g's + 0.5 s' H s

	// -------------------------------------------------------------
	// If sHs <= 0 and gs < 0, there may be a direction of
	// unbounded descent.  Check.
	// -------------------------------------------------------------
	descent  =  gs  <= -epsx;
	negCurv  =  sHs <=  epsx;

	if ( descent  &&  negCurv ) {

	    unbounded = 1;
	    for (j = 0; j < n; j++) {
		if (s[j] > eps)
		    unbounded = bu[j] >=   BigNum;

		else if (s[j] < eps)
		    unbounded = bl[j] <= - BigNum;

		if (!unbounded)
                    break;
	    }

	    if (unbounded)
                return BCLS_EXIT_UNBND;
	}

	// -------------------------------------------------------------
	// Check if the sufficient decrease criteria was satisfied.
	// Otherwise, reduce the step and continue.
	// -------------------------------------------------------------
        PRINT3("I: Proj Backtrack: step = %9.2e  q = %10.2e  gs = %10.2e\n",
               step, q, gs);

	armijo = ( q <= mu * gs );

	if (!armijo) {
	    step = backtrack * step;
	    *nSteps += 1;
	}
    }

    // -----------------------------------------------------------------
    // Force a variable onto its nearest bound: If the armijo
    // condition is satisfied by a truncated step that falls short of
    // the nearest breakpoint, then we need to be a bit more
    // aggressive and continue the step until the nearest boundary.
    //
    // 17 Dec 06: Removed this step.  Why was I doing it?!
    // -----------------------------------------------------------------
#ifdef UNDEF
    if ( step < 1.0  &&  step < BreakMin ) {
        step = BreakMin;
        PRINT3("I: Proj Backtrack: Force index %4i onto bound\n",
               jBreakMin );
    }
#endif

    // Compute the final projected step.
    bcls_project_step( n, s, step, x, dx, bl, bu );

    // Take the step: x = x + s.
    cblas_daxpy( n, 1.0, s, 1, x, 1 );

    // Push x onto bounds that are very close.
    for (j = 0; j < n; j++) {
        if      (x[j] < bl[j] + 10*eps)
                 x[j] = bl[j];
        else if (x[j] > bu[j] - 10*eps)
                 x[j] = bu[j];
    }

    return 0;
}

max_stars_repo_name: echristakopoulou/glslim

max_issues_repo_name: echristakopoulou/glslim

max_forks_repo_name: echristakopoulou/glslim


text:
/*
Copyright (c) 2015, Patrick Weltevrede
All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_sf.h>
#include "psrsalsa.h"
#define AmoebaAlgorithm 0
int readVonMisesModel(char *filename, vonMises_collection_definition *components, verbose_definition verbose)
{
  int i;
  FILE *fin;
  if(verbose.verbose) printf("Opening %s for reading\n", filename);
  fin = fopen(filename, "r");
  if(fin == NULL) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR readVonMisesModel: Cannot open %s", filename);
    return 0;
  }
  for(components->nrcomponents = 0; components->nrcomponents < maxNrVonMisesComponents; (components->nrcomponents)++) {
    i = fscanf(fin, "%lf %lf %lf", &(components->centre[components->nrcomponents]), &(components->concentration[components->nrcomponents]), &(components->height[components->nrcomponents]));
    if(i != 3) {
      break;
    }else if(verbose.verbose) {
      printf("  component %d: phase=%e concentration=%e amplitude=%e\n", components->nrcomponents, components->centre[components->nrcomponents], components->concentration[components->nrcomponents], components->height[components->nrcomponents]);
    }
  }
  if(verbose.verbose) printf("Closing %s\n", filename);
  fclose(fin);
  if(components->nrcomponents == 0) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR readVonMisesModel: No components found in %s", filename);
    return 0;
  }
  return 1;
}
double calcVonMisesFunction2(double centre, double concentration, double height, double phase, double shift)
{
  double y;
  y = exp((cos(2.0*M_PI*(phase-centre-shift))-1.0)*concentration) * height;
  return y;
}
double integrateVonMisesFunction2(double concentration, double height)
{
  double area;
  area = 2*M_PI*height*gsl_sf_bessel_I0_scaled(concentration);
  return area;
}
double widthVonMisesFunction2(double concentration, double ampfrac)
{
  double value = log(ampfrac)/concentration + 1.0;
  if(value < -1.0 || value >= 1.0) {
    return sqrt(-1);
  }
  return 2.0*acos(value);
}
double calcVonMisesFunction(vonMises_collection_definition *components, double phase, double shift)
{
  int n;
  double y;
  y = 0;
  if(components->nrcomponents > 0) {
    for(n = 0; n < components->nrcomponents; n++) {
      y += calcVonMisesFunction2(components->centre[n], components->concentration[n], components->height[n], phase, shift);
    }
  }
  return y;
}
void calcVonMisesProfile(vonMises_collection_definition *components, int nrbins, float *profile, double shift, int normalize)
{
  int i;
  double x, Imax = -1;
  for(i = 0; i < nrbins; i++) {
    profile[i] = 0;
  }
  for(i = 0; i < nrbins; i++) {
    x = i/(double)nrbins;
    profile[i] = calcVonMisesFunction(components, x, shift);
  }
  if(normalize) {
    for(i = 0; i < nrbins; i++) {
      if(profile[i] > Imax)
 Imax = profile[i];
    }
    for(i = 0; i < nrbins; i++) {
      profile[i] /= Imax;
    }
  }
}
float correlateVonMisesFunction(vonMises_collection_definition *components, int nrbins, float *profile, verbose_definition verbose)
{
  float correl_max, *profile2;
  int i;
  int ishift;
  if(verbose.verbose) {
    for(i = 0; i < verbose.indent; i++)
      printf(" ");
    printf("Correlating template with profile\n");
    verbose.indent += 2;
  }
  profile2 = (float *)malloc(nrbins*sizeof(float));
  if(profile2 == NULL) {
    fflush(stdout);
    printerror(verbose.debug, "ERROR correlateVonMisesFunction: Memory allocation error.");
    return 0;
  }
  calcVonMisesProfile(components, nrbins, profile2, 0, 0);
  find_peak_correlation(profile, profile2, nrbins, 0, 1, 1, 1, &ishift, &correl_max, verbose);
  if(verbose.verbose) {
    for(i = 0; i < verbose.indent; i++)
      printf(" ");
    printf("Found a shift of %d bins (%f phase) and max/min correlation is %f.\n", ishift, ishift/(float)nrbins, correl_max);
    for(i = 0; i < verbose.indent; i++)
      printf(" ");
    printf("so you could do a shift by %d bins (%f phase) to rotate profile to align with model.\n", -ishift, -ishift/(float)nrbins);
  }
  free(profile2);
  return ishift/(float)nrbins;
}

max_stars_repo_name: weltevrede/psrsalsa

max_issues_repo_name: weltevrede/psrsalsa

max_forks_repo_name: weltevrede/psrsalsa


text:
// stdafx.h : include file for standard system include files,
// or project specific include files that are used frequently, but
// are changed infrequently
//

#pragma once

#include "targetver.h"
#define _NO_SCRIPT_GUIDS 0
//#define WIN32_LEAN_AND_MEAN             // Exclude rarely-used stuff from Windows headers
// Windows Header Files
#include <windows.h>
#include <d3d11.h>

// C RunTime Header Files
#include <stdlib.h>
#include <malloc.h>
#include <memory.h>
#include <tchar.h>
#include <guiddef.h>

// Spectre
#include <NativeRenderer/Engine.h>
#include <NativeRenderer/Resources/IndexBuffer.h>
#include <NativeRenderer/Resources/ShaderProgram.h>
#include <NativeRenderer/Resources/VertexBuffer.h>
#include <NativeRendererD3D11/RendererD3D11.h>
#include <NativeRendererD3D11/RenderOutputD3D11.h>
#include <Framework/PerformanceTraceLogging.h>
//
// angle
#include <common/platform.h>
#include <angle_gl.h>
#include <EGL/egl.h>
#include <EGL/eglext.h>
#include <libANGLE/Context.h>
#include <libANGLE/renderer/d3d/d3d11/Context11.h>
//
// Chakra
#include <jsrt.h>
//
//// C++/WinRT
//#include <winrt/Windows.ApplicationModel.h>
//#include <winrt/Windows.Graphics.Display.h>
//#include <winrt/Windows.UI.Core.h>

// Arcana
#include <arcana/macros.h>
#include <arcana/string.h>
#include <arcana/experimental/array.h>
#include <arcana/threading/cancellation.h>
#include <arcana/threading/coroutine.h>
#include <arcana/threading/dispatcher.h>
#include <arcana/threading/task.h>
#include <arcana/threading/task_conversions.h>

//// GSL
//#include <gsl/gsl>
//
// C++ standard library
#include <algorithm>
#include <cstddef>
#include <exception>
#include <functional>
#include <memory>
#include <optional>
#include <sstream>
#include <string>
#include <string_view>
#include <unordered_map>
#include <utility>
#include <vector>
#include <variant>

// C standard library
#include <assert.h>

// N-API
#define NODE_ADDON_API_DISABLE_DEPRECATED
#include <napi.h>

max_forks_repo_name: ryantrem/BabylonNative