text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#include <mpi.h>
#include <StGermain/StGermain.h>
#include <StgDomain/StgDomain.h>
#include <StgFEM/StgFEM.h>
#include "types.h"
#include <petsc.h>
#include <petscvec.h>
#include <petscmat.h>
#include <petscksp.h>
#include <petscpc.h>
#include <petscsnes.h>
#include <petscis.h>
#include <petscviewer.h>
#include <petscsys.h>
#include <petscversion.h>

#if ( (PETSC_VERSION_MAJOR >= 3) && (PETSC_VERSION_MINOR >=3) )
  #if (PETSC_VERSION_MINOR >=6)
     #include "petsc/private/kspimpl.h"
  #else
     #include "petsc-private/kspimpl.h"   /*I "petscksp.h" I*/
  #endif
#else
  #include "private/kspimpl.h"   /*I "petscksp.h" I*/
#endif

//#include "ksptypes.h"
#include "ksp-register.h"
#include "StokesBlockKSPInterface.h"

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

/* Macro for checking number integrity - i.e. checks if number is infinite or "not a number" */
#define SBKSP_isGoodNumber( number ) ( (! isnan( number ) ) && ( ! isinf( number ) ) )
#define SBKSP_GetPetscMatrix( matrix ) ( (Mat)(matrix) )
#define SBKSP_GetPetscVector( vector ) ( (Vec)(vector) )
#define SBKSP_GetPetscKSP( solver ) ( (KSP)(solver)  )
PetscErrorCode _BlockSolve( void* solver, void* _stokesSLE );

const Type StokesBlockKSPInterface_Type = "StokesBlockKSPInterface";

void* _StokesBlockKSPInterface_DefaultNew( Name name ) {
    SizeT                                              _sizeOfSelf = sizeof(StokesBlockKSPInterface);
    Type                                                      type = StokesBlockKSPInterface_Type;
    Stg_Class_DeleteFunction*                              _delete = _SLE_Solver_Delete;
    Stg_Class_PrintFunction*                                _print = _SLE_Solver_Print;
    Stg_Class_CopyFunction*                                  _copy = _SLE_Solver_Copy;
    Stg_Component_BuildFunction*                            _build = _StokesBlockKSPInterface_Build;
    Stg_Component_InitialiseFunction*                  _initialise = _StokesBlockKSPInterface_Initialise;
    Stg_Component_ExecuteFunction*                        _execute = _SLE_Solver_Execute;
    Stg_Component_DestroyFunction*                        _destroy = _SLE_Solver_Destroy;
    SLE_Solver_GetResidualFunc*                       _getResidual = NULL;
    Stg_Component_DefaultConstructorFunction*  _defaultConstructor = _StokesBlockKSPInterface_DefaultNew;
    Stg_Component_ConstructFunction*                    _construct = _StokesBlockKSPInterface_AssignFromXML;
    SLE_Solver_SolverSetupFunction*                   _solverSetup = _StokesBlockKSPInterface_SolverSetup;
    SLE_Solver_SolveFunction*                               _solve = _StokesBlockKSPInterface_Solve;

    AllocationType  nameAllocationType = NON_GLOBAL /* default value NON_GLOBAL */;
    return (void*) _StokesBlockKSPInterface_New( STOKESBLOCKKSPINTERFACE_PASSARGS );
}

/* Creation implementation / Virtual constructor */
/* Set up function pointers */
StokesBlockKSPInterface* _StokesBlockKSPInterface_New( STOKESBLOCKKSPINTERFACE_DEFARGS )
{
    StokesBlockKSPInterface* self;
    /* Allocate memory */
    assert( _sizeOfSelf >= sizeof(StokesBlockKSPInterface) );

    self = (StokesBlockKSPInterface*) _SLE_Solver_New( SLE_SOLVER_PASSARGS );

    /* Virtual info */
    return self;
}

void _StokesBlockKSPInterface_Init(
		StokesBlockKSPInterface*      self,
		StiffnessMatrix*   preconditioner,
		Stokes_SLE *       st_sle,
		PETScMGSolver *    mg,
        Name   filename,
        char * string,
		StiffnessMatrix*  k2StiffMat,
		StiffnessMatrix*  mStiffMat,
		ForceVector*	  f2ForceVec,
		ForceVector*	  jForceVec,
		double            penaltyNumber,
		double            hFactor,
		StiffnessMatrix*  vmStiffMat,
		ForceVector*	  vmForceVec  )
{
	self->preconditioner = preconditioner;
	self->st_sle = st_sle;
	self->mg     = mg;
    self->optionsFile = filename;
    self->optionsString = string;
	self->k2StiffMat = k2StiffMat;
	self->f2ForceVec = f2ForceVec;
	self->penaltyNumber = penaltyNumber;
	self->hFactor       = hFactor;
	self->mStiffMat = mStiffMat;
	self->jForceVec = jForceVec;
	self->vmStiffMat = vmStiffMat;
	self->vmForceVec = vmForceVec;
	/* add the vecs and matrices to the Base SLE class's dynamic lists, so they can be
	initialised and built properly */

	/*
	if (k2StiffMat )
	    SystemLinearEquations_AddStiffnessMatrix( st_sle, k2StiffMat );

	if (f2ForceVec )
	    SystemLinearEquations_AddForceVector( st_sle, f2ForceVec );

	if (mStiffMat )
	    SystemLinearEquations_AddStiffnessMatrix( st_sle, mStiffMat );

	if (jForceVec )
	    SystemLinearEquations_AddForceVector( st_sle, jForceVec );

	if (vmStiffMat )
	    SystemLinearEquations_AddStiffnessMatrix( st_sle, vmStiffMat );

	if (vmForceVec )
	    SystemLinearEquations_AddForceVector( st_sle, vmForceVec );
    */
}

void _StokesBlockKSPInterface_Build( void* solver, void* sle ) {/* it is the sle here being passed in*/
	StokesBlockKSPInterface*	self  = (StokesBlockKSPInterface*)solver;

	Stream_IndentBranch( StgFEM_Debug );

	/* Build Preconditioner */
	if ( self->preconditioner ) {
		Stg_Component_Build( self->preconditioner, sle, False );
		SystemLinearEquations_AddStiffnessMatrix( self->st_sle, self->preconditioner );

	}
	if( self->mg ){
	    Stg_Component_Build( self->mg, sle, False );
	}

	Stream_UnIndentBranch( StgFEM_Debug );
}

void _StokesBlockKSPInterface_AssignFromXML( void* solver, Stg_ComponentFactory* cf, void* data ) {
	StokesBlockKSPInterface* self         = (StokesBlockKSPInterface*) solver;
	//double                  tolerance;
	//Iteration_Index         maxUzawaIterations, minUzawaIterations;
	StiffnessMatrix*  preconditioner;
	StiffnessMatrix*  k2StiffMat;
	ForceVector*	  f2ForceVec;
	StiffnessMatrix*  mStiffMat;
	ForceVector*	  jForceVec;
	double            penaltyNumber;
	double            hFactor;
	StiffnessMatrix*  vmStiffMat;
	ForceVector*	  vmForceVec;
	//Bool                    useAbsoluteTolerance;
	//Bool                    monitor;
	Stokes_SLE *            st_sle;
	PETScMGSolver *         mg;
	//Name                filename;
	//char* 		        string;

	_SLE_Solver_AssignFromXML( self, cf, data );

	preconditioner = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"Preconditioner", StiffnessMatrix, False, data  );
	st_sle  = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"stokesEqn", Stokes_SLE, True, data  );
	mg      = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"mgSolver", PETScMGSolver, False, data);
	k2StiffMat = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"2ndStressTensorMatrix", StiffnessMatrix, False, data  );
	f2ForceVec = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"2ndForceVector", ForceVector, False, data  );
	penaltyNumber   = Stg_ComponentFactory_GetDouble( cf, self->name, (Dictionary_Entry_Key)"penaltyNumber", 0.0  );
	hFactor         = Stg_ComponentFactory_GetDouble( cf, self->name, (Dictionary_Entry_Key)"hFactor", 0.0  );
	mStiffMat = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"MassMatrix", StiffnessMatrix, False, data  );
	jForceVec = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"JunkForceVector", ForceVector, False, data  );
	vmStiffMat = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"VelocityMassMatrix", StiffnessMatrix, False, data  );
	vmForceVec = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"VMassForceVector", ForceVector, False, data  );

	_StokesBlockKSPInterface_Init( self, preconditioner, st_sle, mg, NULL, NULL, k2StiffMat, mStiffMat,
                                   f2ForceVec, jForceVec, penaltyNumber, hFactor, vmStiffMat, vmForceVec);

}

void _StokesBlockKSPInterface_Initialise( void* solver, void* data ) {
	StokesBlockKSPInterface* self = (StokesBlockKSPInterface*) solver;
	Stokes_SLE*             sle  = (Stokes_SLE*)  self->st_sle;

	/* Initialise Parent */
	_SLE_Solver_Initialise( self, sle );

	KSPRegisterAllKSP("Solvers/KSPSolvers/src");
}

/* SolverSetup */

void _StokesBlockKSPInterface_SolverSetup( void* solver, void* stokesSLE ) {
	StokesBlockKSPInterface* self = (StokesBlockKSPInterface*) solver;
	//Stokes_SLE*             sle  = (Stokes_SLE*)             stokesSLE;

 	Journal_DPrintf( self->debug, "In %s:\n", __func__ );
	Stream_IndentBranch( StgFEM_Debug );

	Stream_UnIndentBranch( StgFEM_Debug );
}
void SBKSP_SetSolver( void* solver, void* stokesSLE ) {
	SLE_Solver* self = (SLE_Solver*) solver;
	Stokes_SLE*              sle  = (Stokes_SLE*) stokesSLE;

    sle->solver=self;

}
void SBKSP_SetPenalty( void* solver, double penalty ) {
  StokesBlockKSPInterface* self = (StokesBlockKSPInterface*) solver;
  self->penaltyNumber=penalty;
}

int SBKSP_GetPressureIts(void *solver){
  StokesBlockKSPInterface* self = (StokesBlockKSPInterface*) solver;
  return self->stats.pressure_its;
}

/***********************************************************************************************************/
/***********************************************************************************************************/
/***********************************************************************************************************/

void SBKSP_GetStokesOperators(
		Stokes_SLE *stokesSLE,
		Mat *K,Mat *G,Mat *D,Mat *C,Mat *approxS,
		Vec *f,Vec *h,Vec *u,Vec *p )
{

	*K = *G = *D = *C = PETSC_NULL;
	if (stokesSLE->kStiffMat){      *K = SBKSP_GetPetscMatrix( stokesSLE->kStiffMat->matrix );     }
	if (stokesSLE->gStiffMat){      *G = SBKSP_GetPetscMatrix( stokesSLE->gStiffMat->matrix );     }
	if (stokesSLE->dStiffMat){      *D = SBKSP_GetPetscMatrix( stokesSLE->dStiffMat->matrix );     }
	if (stokesSLE->cStiffMat){      *C = SBKSP_GetPetscMatrix( stokesSLE->cStiffMat->matrix );     }

	/* preconditioner */
	*approxS = PETSC_NULL;
	if( ((StokesBlockKSPInterface*)stokesSLE->solver)->preconditioner ) {
		StiffnessMatrix *preconditioner;

		preconditioner = ((StokesBlockKSPInterface*)stokesSLE->solver)->preconditioner;
		*approxS = SBKSP_GetPetscMatrix( preconditioner->matrix );
	}

	*f = *h = PETSC_NULL;
	if (stokesSLE->fForceVec){      *f = SBKSP_GetPetscVector( stokesSLE->fForceVec->vector );     }
	if (stokesSLE->hForceVec){      *h = SBKSP_GetPetscVector( stokesSLE->hForceVec->vector );     }

	*u = *p = PETSC_NULL;
	if (stokesSLE->uSolnVec){       *u = SBKSP_GetPetscVector( stokesSLE->uSolnVec->vector );      }
	if (stokesSLE->pSolnVec){       *p = SBKSP_GetPetscVector( stokesSLE->pSolnVec->vector );      }

}
/***********************************************************************************************************/
/***********************************************************************************************************/
/***********************************************************************************************************/
/* Sets up Solver to be a custom ksp (KSP_BSSCR) solve by default: */
void _StokesBlockKSPInterface_Solve( void* solver, void* _stokesSLE ) {
  StokesBlockKSPInterface* self    = (StokesBlockKSPInterface*)solver;
  PetscLogDouble flopsA,flopsB;
  PetscTruth found, get_flops;

  found = PETSC_FALSE;
  get_flops = PETSC_FALSE;
  PetscOptionsGetTruth( PETSC_NULL, "-get_flops", &get_flops, &found);
  if(get_flops){
    PetscGetFlops(&flopsA); }

  _BlockSolve(solver, _stokesSLE);

  if(get_flops){
    PetscGetFlops(&flopsB);
    self->stats.total_flops=(double)(flopsB-flopsA); }
}
PetscErrorCode _BlockSolve( void* solver, void* _stokesSLE ) {
  Stokes_SLE*  stokesSLE  = (Stokes_SLE*)_stokesSLE;
  StokesBlockKSPInterface* Solver    = (StokesBlockKSPInterface*)solver;

  /* Create shortcuts to stuff needed on sle */
  Mat       K;
  Mat       G;
  Mat       Gt;
  Mat       D;
  Mat       C;
  Mat       approxS;
  Vec       u;
  Vec       p;
  Vec       f;
  Vec       h;
  Mat stokes_P;
  Mat stokes_A;
  Vec stokes_x;
  Vec stokes_b;
  Mat a[2][2];
  Vec x[2];
  Vec b[2];
  KSP stokes_ksp;
  PC  stokes_pc;
  PetscTruth sym,flg;
  PetscErrorCode ierr;

  PetscInt   N,n;

  SBKSP_GetStokesOperators( stokesSLE, &K,&G,&D,&C, &approxS, &f,&h, &u,&p );

  /* create Gt */
  if( !D ) {
    ierr = MatTranspose( G, MAT_INITIAL_MATRIX, &Gt);CHKERRQ(ierr);
    sym = PETSC_TRUE;
    Solver->DIsSym = sym;
  }
  else {
    Gt = D;
    sym = PETSC_FALSE;
    Solver->DIsSym = sym;
  }
  flg=PETSC_FALSE;
  PetscOptionsHasName(PETSC_NULL,"-use_petsc_ksp",&flg);
  if (flg) {
    if( !C ) {
      /* Everything in this bracket, dependent on !C, is to build
         a matrix with diagonals of 0 for C the previous comment ways

      need a 'zero' matrix to keep fieldsplit happy in petsc? */
      MatType mtype;
      Vec V;
      //MatGetSize( G, &M, &N );
      VecGetSize(p, &N);
      VecGetLocalSize( p, &n );
      MatCreate( PetscObjectComm((PetscObject) K), &C );
      MatSetSizes( C, PETSC_DECIDE ,PETSC_DECIDE, N, N );
#if (((PETSC_VERSION_MAJOR==3) && (PETSC_VERSION_MINOR>=3)) || (PETSC_VERSION_MAJOR>3) )
      MatSetUp(C);
#endif
      MatGetType( G, &mtype );
      MatSetType( C, mtype );
      MatGetVecs( G, &V, PETSC_NULL );
      VecSet(V, 0.0);
      //VecSet(h, 1.0);
      ierr = VecAssemblyBegin( V );CHKERRQ(ierr);
      ierr = VecAssemblyEnd  ( V );CHKERRQ(ierr);
      ierr = MatDiagonalSet(C,V,INSERT_VALUES);CHKERRQ(ierr);
      ierr = MatAssemblyBegin( C, MAT_FINAL_ASSEMBLY );CHKERRQ(ierr);
      ierr = MatAssemblyEnd  ( C, MAT_FINAL_ASSEMBLY );CHKERRQ(ierr);
    }
  }
  a[0][0]=K;  a[0][1]=G;
  a[1][0]=Gt; a[1][1]=C;
  ierr = MatCreateNest(PetscObjectComm((PetscObject) K), 2, NULL, 2, NULL, (Mat *)a, &stokes_A);CHKERRQ(ierr);
  ierr = MatAssemblyBegin( stokes_A, MAT_FINAL_ASSEMBLY );CHKERRQ(ierr);
  ierr = MatAssemblyEnd( stokes_A, MAT_FINAL_ASSEMBLY );CHKERRQ(ierr);

  x[0]=u;
  x[1]=p;
  ierr = VecCreateNest(PetscObjectComm((PetscObject) u), 2, NULL, x, &stokes_x);CHKERRQ(ierr);
  ierr = VecAssemblyBegin( stokes_x );CHKERRQ(ierr);
  ierr = VecAssemblyEnd( stokes_x);CHKERRQ(ierr);

  b[0]=f;
  b[1]=h;
  ierr = VecCreateNest(PetscObjectComm((PetscObject) f), 2, NULL, b, &stokes_b);CHKERRQ(ierr);
  ierr = VecAssemblyBegin( stokes_b );CHKERRQ(ierr);
  ierr = VecAssemblyEnd( stokes_b);CHKERRQ(ierr);

  /* if( approxS ) { */
  /*   a[0][0]=K;    a[0][1]=G; */
  /*   a[1][0]=NULL; a[1][1]=approxS; */
  /*   ierr = MatCreateNest(PetscObjectComm((PetscObject) K), 2, NULL, 2, NULL, (Mat *)a, &stokes_P);CHKERRQ(ierr); */
  /*   ierr = MatAssemblyBegin( stokes_P, MAT_FINAL_ASSEMBLY );CHKERRQ(ierr); */
  /*   ierr = MatAssemblyEnd( stokes_P, MAT_FINAL_ASSEMBLY );CHKERRQ(ierr); */
  /* } */
  /* else { */
    stokes_P = stokes_A;
  /* } */

  /* probably should make a Destroy function for these two */
  /* Update options from file and/or string here so we can change things on the fly */
  //PetscOptionsInsertFile(PETSC_COMM_WORLD, Solver->optionsFile, PETSC_FALSE);
  //PetscOptionsInsertString(Solver->optionsString);

  ierr = KSPCreate( PETSC_COMM_WORLD, &stokes_ksp );CHKERRQ(ierr);
  Stg_KSPSetOperators( stokes_ksp, stokes_A, stokes_P, SAME_NONZERO_PATTERN );
  ierr = KSPSetType( stokes_ksp, "bsscr" );/* i.e. making this the default solver : calls KSPCreate_XXX */CHKERRQ(ierr);

  ierr = KSPGetPC( stokes_ksp, &stokes_pc );CHKERRQ(ierr);
  ierr = PCSetType( stokes_pc, PCNONE );CHKERRQ(ierr);
  ierr = KSPSetInitialGuessNonzero( stokes_ksp, PETSC_TRUE );CHKERRQ(ierr);
  ierr = KSPSetFromOptions( stokes_ksp );CHKERRQ(ierr);

  /*
    Doing this so the KSP Solver has access to the StgFEM Multigrid struct (PETScMGSolver).
    As well as any custom stuff on the Stokes_SLE struct
  */
  if( stokes_ksp->data ){/* then ksp->data has been created in a KSpSetUp_XXX function */
    /* testing for our KSP types that need the data that is on Solver... */
    /* for the moment then, this function not completely agnostic about our KSPs */
    //if(!strcmp("bsscr",stokes_ksp->type_name)){/* if is bsscr then set up the data on the ksp */
    flg=PETSC_FALSE;
    PetscOptionsHasName(PETSC_NULL,"-use_petsc_ksp",&flg);
    if (!flg) {
      ((KSP_COMMON*)(stokes_ksp->data))->st_sle         = Solver->st_sle;
      ((KSP_COMMON*)(stokes_ksp->data))->mg             = Solver->mg;
      ((KSP_COMMON*)(stokes_ksp->data))->DIsSym         = Solver->DIsSym;
      ((KSP_COMMON*)(stokes_ksp->data))->preconditioner = Solver->preconditioner;
      ((KSP_COMMON*)(stokes_ksp->data))->solver         = Solver;
    }
  }

  ierr = KSPSolve( stokes_ksp, stokes_b, stokes_x );CHKERRQ(ierr);

  Stg_KSPDestroy(&stokes_ksp );
  //if( ((StokesBlockKSPInterface*)stokesSLE->solver)->preconditioner )
  if(stokes_P != stokes_A) { Stg_MatDestroy(&stokes_P ); }

  Stg_MatDestroy(&stokes_A );

  Stg_VecDestroy(&stokes_x);
  Stg_VecDestroy(&stokes_b);

  if(!D){ Stg_MatDestroy(&Gt); }
  if(C && (stokesSLE->cStiffMat->matrix != C) ){ Stg_MatDestroy(&C); }

  PetscFunctionReturn(0);
}

max_stars_repo_name: StuartRClark/mantle

max_issues_repo_name: StuartRClark/mantle

max_forks_repo_name: StuartRClark/mantle


text:
/*
  Copyright [2021] [IBM Corporation]
  Licensed under the Apache License, Version 2.0 (the "License");
  you may not use this file except in compliance with the License.
  You may obtain a copy of the License at
  http://www.apache.org/licenses/LICENSE-2.0
  Unless required by applicable law or agreed to in writing, software
  distributed under the License is distributed on an "AS IS" BASIS,
  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  See the License for the specific language governing permissions and
  limitations under the License.
*/

#ifndef _MCAS_NUMA_NODE_MASK__H_
#define _MCAS_NUMA_NODE_MASK__H_

#include <numa.h> /* bitmask */
#include <common/string_view.h>
#include <gsl/pointers> /* not_null */
#include <memory> /* unique_ptr */

struct numa_node_mask_dtor
{
	void operator()(bitmask *b) { numa_free_nodemask(b); }
};

struct numa_node_mask
{
	numa_node_mask(common::string_view mask_)
		: _mask(std::unique_ptr<bitmask, numa_node_mask_dtor>(mask_not_null(std::string(mask_))))
	{
	}
	numa_node_mask(const bitmask *b)
		: _mask(std::unique_ptr<bitmask, numa_node_mask_dtor>(numa_allocate_nodemask()))
	{
		copy_bitmask_to_bitmask(const_cast<bitmask *>(b), get());
	}
	const bitmask *get() const { return _mask.get().get(); }
	bitmask *get() { return _mask.get().get(); }
	/* get first 64 bits in mask */
	uint64_t get64()
	{
		uint64_t r = 0;
		for ( unsigned i = 0; i != 64; ++i )
		{
			if ( numa_bitmask_isbitset(get(), i) )
			{
				r |= ( 1U << i );
			}
		}
		return r;
	}
private:
	static bitmask *mask_not_null(common::string_view mask_)
	{
		auto m = numa_parse_nodestring(std::string(mask_).c_str());
		return m ? m : numa_allocate_nodemask();
	}
	gsl::not_null<std::unique_ptr<bitmask, numa_node_mask_dtor>> _mask;
};

#endif

ext: h

max_stars_repo_name: IBM/artemis

max_issues_repo_name: IBM/artemis

max_forks_repo_name: IBM/artemis


text:
/****************************************************************
 *
 * Copyright (C) Max Planck Institute
 * for Biological Cybernetics, Tuebingen, Germany
 *
 * Author: Gabriele Lohmann, 2015
 *
 * 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 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.
 *
 * 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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 *****************************************************************/

/*
** trial average using spline interpolation
**
** G.Lohmann, May 2014
*/
#include <viaio/Vlib.h>
#include <viaio/VImage.h>
#include <viaio/mu.h>
#include <viaio/option.h>

#include <gsl/gsl_cblas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_spline.h>

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>

#define ABS(x) ((x) > 0 ? (x) : -(x))
#define SQR(x) ((x)*(x))

#define LEN     10000   /* buffer length        */
#define NTRIALS 10000   /* max number of trials */

typedef struct SpointStruct {
  VShort x;
  VShort y;
  VShort z;
} SPoint;

typedef struct TrialStruct {
  int   id;
  float onset;
  float duration;
  float height;
} Trial;

Trial trial[NTRIALS];
int ntrials=0;
int nevents=0;

int test_ascii(int val)
{
  if (val >= 'a' && val <= 'z') return 1;
  if (val >= 'A' && val <= 'Z') return 1;
  if (val >= '0' && val <= '9') return 1;
  if (val ==  ' ') return 1;
  if (val == '\0') return 1;
  if (val == '\n') return 1;
  if (val == '\r') return 1;
  if (val == '\t') return 1;
  if (val == '\v') return 1;
  return 0;
}

/* parse design file */
void ReadDesign(VString designfile)
{
  FILE *fp=NULL;
  int  i,j,k,id;
  char buf[LEN];
  float onset=0,duration=0,height=0;

  fp = fopen(designfile,"r");
  if (!fp) VError(" error opening design file %s",designfile);

  i = ntrials = nevents = 0;
  while (!feof(fp)) {
    for (j=0; j<LEN; j++) buf[j] = '\0';
    if (fgets(buf,LEN,fp) == NULL) break;
    if (strlen(buf) < 2) continue;
    if (buf[0] == '%' || buf[0] == '#') continue;
    if (! test_ascii((int)buf[0])) VError(" input file must be a text file");

    /* remove non-alphanumeric characters */
    for (j=0; j<strlen(buf); j++) {
      k = (int)buf[j];
      if (!isgraph(k) && buf[j] != '\n' && buf[j] != '\r' && buf[j] != '\0') {
	buf[j] = ' ';
      }
      if (buf[j] == '\v') buf[j] = ' '; /* remove tabs */
      if (buf[j] == '\t') buf[j] = ' ';
    }

    if (sscanf(buf,"%d %f %f %f",&id,&onset,&duration,&height) != 4)
      VError(" line %d: illegal input format",i+1);

    if (duration < 0.5 && duration >= -0.0001) duration = 0.5;
    trial[i].id       = id;
    trial[i].onset    = onset;
    trial[i].duration = duration;
    trial[i].height   = height;
    i++;
    if (i > NTRIALS) VError(" too many trials %d",i);

    if (id > nevents) nevents = id;
  }
  fclose(fp);

  ntrials = i;
}

int main (int argc,char *argv[])
{
  static VString  designfile = "";
  static VShort   cond_id = 1;
  static VShort   trial_id = 0;
  static VDouble  temporal_resolution = 1.0;
  static VDouble  start = 0;
  static VDouble  length = 20;
  static VOptionDescRec  options[] = {
    {"design", VStringRepn, 1, & designfile, VRequiredOpt, NULL,"Design file (ascii)" },
    {"cond", VShortRepn, 1, & cond_id, VOptionalOpt, NULL,"Id of experimental condition"},
    {"trial", VShortRepn, 1, & trial_id, VOptionalOpt, NULL,"Id of trial (starts at 0)"},
    {"resolution", VDoubleRepn, 1, & temporal_resolution, VOptionalOpt, NULL," output temporal resolution in secs"},
    {"start", VDoubleRepn, 1, & start, VOptionalOpt, NULL, "start relative to beginning of trial in secs"},
    {"length", VDoubleRepn, 1, & length, VOptionalOpt, NULL, "trial length in seconds"},
  };
  FILE *in_file=NULL,*out_file=NULL;
  VAttrList list=NULL;
  VAttrListPosn posn;
  int i,j,k,slice,row,col;
  int nslices=0,nrows=0,ncols=0,ntimesteps=0;
  double t=0;
  char *prg = GetLipsiaName("vcuttrials");
  fprintf (stderr, "%s\n", prg);

  /*  parse command line */
  VParseFilterCmd (VNumber (options),options,argc,argv,&in_file,&out_file);

  /* get image dimensions, read functional data */
  if (! (list = VReadFile (in_file, NULL))) exit (1);
  fclose(in_file);

  nslices = 0;
  VImage tmp=NULL;
  for (VFirstAttr (list, & posn); VAttrExists (& posn); VNextAttr (& posn)) {
    if (VGetAttrRepn (& posn) != VImageRepn) continue;
    VGetAttrValue (& posn, NULL,VImageRepn, & tmp);
    /* if (VPixelRepn(tmp) != VShortRepn) continue; */
    nslices++;
  }
  if (nslices < 1) VError(" no slices");

  VAttrList geoinfo = VGetGeoInfo(list);
  VImage *src = (VImage *) VCalloc(nslices,sizeof(VImage));
  i = ntimesteps = nrows = ncols = 0;
  for (VFirstAttr (list, & posn); VAttrExists (& posn); VNextAttr (& posn)) {
    if (VGetAttrRepn (& posn) != VImageRepn) continue;
    VGetAttrValue (& posn, NULL,VImageRepn, & src[i]);
    /* if (VPixelRepn(src[i]) != VShortRepn) continue; */
    if (VImageNBands(src[i]) > ntimesteps) ntimesteps = VImageNBands(src[i]);
    if (VImageNRows(src[i])  > nrows)  nrows = VImageNRows(src[i]);
    if (VImageNColumns(src[i]) > ncols) ncols = VImageNColumns(src[i]);
    i++;
  }
  fprintf(stderr," nslices: %d, nrows: %d, ncols: %d,  ntimesteps: %d\n",nslices,nrows,ncols,ntimesteps);

  /* read repetition time */
  double tr = 0;
  if (VGetAttr (VImageAttrList (src[0]), "repetition_time", NULL,
		VDoubleRepn, (VPointer) & tr) != VAttrFound) {
    VError(" attribute 'repetition_time' missing");
  }
  tr /= 1000.0;
  double experiment_duration = tr * (double)ntimesteps;

  /* read design file */
  ReadDesign(designfile);

  /* select trial id */
  int tid=0,j0=-1;
  for (j=0; j<ntrials; j++) {
    if (trial[j].id == cond_id && tid == trial_id) {
      j0 = j;
      break;
    }
    if (trial[j].id == cond_id) tid++;
  }
  if (j0 < 0) VError("trial not found");
  fprintf(stderr," trial onset: %f\n",trial[j0].onset);

  /* adjust start time */
  trial[j0].onset += start;
  if (trial[j0].onset < 0) {
    VWarning(" negative onset, reset to zero");
    trial[j0].onset = 0;
  }
  if (trial[j0].onset + length >= experiment_duration) {
    VWarning(" experiment_duration exceeded",experiment_duration);
    length = experiment_duration - trial[j0].onset-0.5;
    VWarning(" parameter '-length' set to %f\n",length);
  }

  int nt = (int)(length/temporal_resolution + 0.5);
  fprintf(stderr," number of output timesteps: %d\n",nt);

  /* ini output data structs  */
  VImage *dest = (VImage *) VCalloc(nslices, sizeof(VImage));
  VAttrList out_list = VCreateAttrList();
  if (geoinfo != NULL) VSetGeoInfo(geoinfo,out_list);
  for (slice=0; slice<nslices; slice++) {
    dest[slice] = VCreateImage(nt,nrows,ncols,VPixelRepn(src[slice]));
    VFillImage(dest[slice],VAllBands,0);
    VCopyImageAttrs (src[slice], dest[slice]);
    VAppendAttr(out_list,"image",NULL,VImageRepn,dest[slice]);
  }

  /* spline interpolation  */
  double tstep = temporal_resolution;
  double *xx = (double *) VCalloc(ntimesteps,sizeof(double));
  double *yy = (double *) VCalloc(ntimesteps,sizeof(double));

  gsl_interp_accel *acc = gsl_interp_accel_alloc ();
  gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, ntimesteps);

  for (slice=0; slice<nslices; slice++) {
    for (row=0; row<nrows; row++) {
      for (col=0; col<ncols; col++) {

	for (k=0; k<ntimesteps; k++) {
	  xx[k] = tr*((double)k);
	  yy[k] = (double)VGetPixel(src[slice],k,row,col);
	}
	gsl_spline_init (spline, xx, yy, ntimesteps);

	k=0;
	for (t=0; t<length; t += tstep) {
	  if (k >= nt) break;
	  double xi = trial[j].onset + t;
	  double yi = 0;
	  if (xi < experiment_duration && xi >= 0) {
	    if (gsl_spline_eval_e (spline,xi,acc,&yi) == GSL_EDOM) continue;
	  }
	  VSetPixel(dest[slice],k,row,col,yi);
	  k++;
	}
      }
    }
  }

  /* write to disk */
  if (! VWriteFile (out_file, out_list)) exit (1);
  fprintf (stderr, "%s: done.\n", argv[0]);
  exit(0);
}

ext: c

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
/* err/stream.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
 *
 * 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 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.
 *
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <stddef.h>
#include <stdlib.h>
#include <stdio.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_message.h>

FILE * gsl_stream = NULL ;
gsl_stream_handler_t * gsl_stream_handler = NULL;

void
gsl_stream_printf (const char *label, const char *file, int line,
		   const char *reason)
{
  if (gsl_stream == NULL)
    {
      gsl_stream = stderr;
    }
  if (gsl_stream_handler)
    {
      (*gsl_stream_handler) (label, file, line, reason);
      return;
    }
  fprintf (gsl_stream, "gsl: %s:%d: %s: %s\n", file, line, label, reason);

}

gsl_stream_handler_t *
gsl_set_stream_handler (gsl_stream_handler_t * new_handler)
{
  gsl_stream_handler_t * previous_handler = gsl_stream_handler;
  gsl_stream_handler = new_handler;
  return previous_handler;
}

FILE *
gsl_set_stream (FILE * new_stream)
{
  FILE * previous_stream;
  if (gsl_stream == NULL) {
    gsl_stream = stderr;
  }
  previous_stream = gsl_stream;
  gsl_stream = new_stream;
  return previous_stream;
}


text:
// The MIT License (MIT)
//
// Copyright (c) 2018 Mateusz Pusz
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.

// Formatting library for C++ - the core API for char/UTF-8
//
// Copyright (c) 2012 - present, Victor Zverovich
// All rights reserved.
//
// For the license information refer to format.h.

#include <gsl/gsl-lite.hpp>
#include <units/bits/fmt_hacks.h>
#include <concepts>
#include <limits>
#include <string_view>

// most of the below code is based on/copied from libfmt

namespace units::detail {

struct auto_id {};

enum class fmt_align { none, left, right, center };
enum class fmt_sign { none, minus, plus, space };
enum class arg_id_kind { none, index, name };

template<typename Char>
struct fill_t {
private:
  static constexpr size_t max_size = 4 / sizeof(Char);
  // At most one codepoint (so one char32_t or four utf-8 char8_t)
  Char data_[max_size] = {Char{' '}};
  unsigned char size_ = 1;

public:
  constexpr void operator=(std::basic_string_view<Char> s)
  {
    auto size = s.size();
    if (size > max_size) return throw STD_FMT::format_error("invalid fill");
    for (size_t i = 0; i < size; ++i) data_[i] = s[i];
    size_ = static_cast<unsigned char>(size);
  }

  [[nodiscard]] constexpr size_t size() const { return size_; }
  [[nodiscard]] constexpr const Char* data() const { return data_; }

  [[nodiscard]] constexpr Char& operator[](size_t index) { return data_[index]; }
  [[nodiscard]] constexpr const Char& operator[](size_t index) const { return data_[index]; }
};

template<typename T>
inline constexpr bool is_integer = std::is_integral<T>::value && !std::is_same<T, bool>::value &&
                                   !std::is_same<T, char>::value && !std::is_same<T, wchar_t>::value;

template<typename Char>
[[nodiscard]] constexpr bool is_ascii_letter(Char c)
{
  return (c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z');
}

// Converts a character to ASCII. Returns a number > 127 on conversion failure.
template<std::integral Char>
[[nodiscard]] constexpr Char to_ascii(Char value)
{
  return value;
}

template<typename Char>
  requires std::is_enum_v<Char>
[[nodiscard]] constexpr auto to_ascii(Char value) -> std::underlying_type_t<Char> { return value; }

struct width_checker {
  template<typename T>
  [[nodiscard]] constexpr unsigned long long operator()(T value) const
  {
    if constexpr (is_integer<T>) {
      if constexpr (std::numeric_limits<T>::is_signed) {
        if (value < 0) throw STD_FMT::format_error("negative width");
      }
      return static_cast<unsigned long long>(value);
    } else {
      throw STD_FMT::format_error("width is not integer");
    }
  }
};

struct precision_checker {
  template<typename T>
  [[nodiscard]] constexpr unsigned long long operator()(T value) const
  {
    if constexpr (is_integer<T>) {
      if constexpr (std::numeric_limits<T>::is_signed) {
        if (value < 0) throw STD_FMT::format_error("negative precision");
      }
      return static_cast<unsigned long long>(value);
    } else {
      throw STD_FMT::format_error("precision is not integer");
    }
  }
};

// Format specifiers for built-in and string types.
template<typename Char>
struct basic_format_specs {
  int width = 0;
  int precision = -1;
  char type = '\0';
  fmt_align align : 4 = fmt_align::none;
  fmt_sign sign : 3 = fmt_sign::none;
  bool alt : 1 = false;  // Alternate form ('#').
  bool localized : 1 = false;
  fill_t<Char> fill;
};

// Format specifiers with width and precision resolved at formatting rather
// than parsing time to allow re-using the same parsed specifiers with
// different sets of arguments (precompilation of format strings).
template<typename Char>
struct dynamic_format_specs : basic_format_specs<Char> {
  int dynamic_width_index = -1;
  int dynamic_precision_index = -1;
};

[[nodiscard]] constexpr int verify_dynamic_arg_index_in_range(size_t idx)
{
  if (idx > static_cast<size_t>(std::numeric_limits<int>::max())) {
    throw STD_FMT::format_error("Dynamic width or precision index too large.");
  }
  return static_cast<int>(idx);
}

template<typename CharT>
[[nodiscard]] constexpr int on_dynamic_arg(size_t arg_id, STD_FMT::basic_format_parse_context<CharT>& context)
{
  context.check_arg_id(FMT_TO_ARG_ID(arg_id));
  return verify_dynamic_arg_index_in_range(arg_id);
}

template<typename CharT>
[[nodiscard]] constexpr int on_dynamic_arg(auto_id, STD_FMT::basic_format_parse_context<CharT>& context)
{
  return verify_dynamic_arg_index_in_range(FMT_FROM_ARG_ID(context.next_arg_id()));
}

template<class Handler, typename FormatContext>
[[nodiscard]] constexpr int get_dynamic_spec(int index, FormatContext& ctx)
{
  const unsigned long long value =
    STD_FMT::visit_format_arg(Handler{}, ctx.arg(FMT_TO_ARG_ID(static_cast<size_t>(index))));
  if (value > static_cast<unsigned long long>(std::numeric_limits<int>::max())) {
    throw STD_FMT::format_error("number is too big");
  }
  return static_cast<int>(value);
}

// Parses the range [begin, end) as an unsigned integer. This function assumes
// that the range is non-empty and the first character is a digit.
template<std::input_iterator It, std::sentinel_for<It> S>
[[nodiscard]] constexpr It parse_nonnegative_int(It begin, S end, size_t& value)
{
  gsl_Expects(begin != end && '0' <= *begin && *begin <= '9');
  constexpr auto max_int = static_cast<unsigned>(std::numeric_limits<int>::max());
  constexpr auto big_int = max_int / 10u;
  value = 0;

  do {
    if (value > big_int) {
      value = max_int + 1;
      break;
    }
    value = value * 10 + static_cast<unsigned int>(*begin - '0');
    ++begin;
  } while (begin != end && '0' <= *begin && *begin <= '9');

  if (value > max_int) throw STD_FMT::format_error("Number is too big");

  return begin;
}

template<std::input_iterator It, std::sentinel_for<It> S>
[[nodiscard]] constexpr It parse_nonnegative_int(It begin, S end, int& value)
{
  size_t val_unsigned = 0;
  begin = parse_nonnegative_int(begin, end, val_unsigned);
  // Never invalid because parse_nonnegative_integer throws an error for values that don't fit in signed integers
  value = static_cast<int>(val_unsigned);
  return begin;
}

template<std::input_iterator It, std::sentinel_for<It> S, typename IDHandler>
[[nodiscard]] constexpr It do_parse_arg_id(It begin, S end, IDHandler&& handler)
{
  gsl_Expects(begin != end);
  auto c = *begin;
  if (c >= '0' && c <= '9') {
    size_t index = 0;
    if (c != '0')
      begin = parse_nonnegative_int(begin, end, index);
    else
      ++begin;
    if (begin == end || (*begin != '}' && *begin != ':'))
      throw STD_FMT::format_error("invalid format string");
    else
      handler(index);
    return begin;
  }
  throw STD_FMT::format_error("invalid format string");
}

template<std::input_iterator It, std::sentinel_for<It> S, typename IDHandler>
[[nodiscard]] constexpr It parse_arg_id(It begin, S end, IDHandler&& handler)
{
  auto c = *begin;
  if (c != '}' && c != ':') return do_parse_arg_id(begin, end, handler);
  handler();
  return begin;
}

template<std::input_iterator It, std::sentinel_for<It> S, typename Handler>
[[nodiscard]] constexpr It parse_sign(It begin, S end, Handler&& handler)
{
  gsl_Expects(begin != end);
  switch (to_ascii(*begin)) {
    case '+':
      handler.on_sign(fmt_sign::plus);
      ++begin;
      break;
    case '-':
      handler.on_sign(fmt_sign::minus);
      ++begin;
      break;
    case ' ':
      handler.on_sign(fmt_sign::space);
      ++begin;
      break;
    default:
      break;
  }
  return begin;
}

template<std::input_iterator It, std::sentinel_for<It> S, typename Handler>
[[nodiscard]] constexpr It parse_width(It begin, S end, Handler&& handler)
{
  struct width_adapter {
    Handler& handler;
    constexpr void operator()() { handler.on_dynamic_width(auto_id{}); }
    constexpr void operator()(size_t id) { handler.on_dynamic_width(id); }
  };

  gsl_Expects(begin != end);
  if ('0' <= *begin && *begin <= '9') {
    int width = 0;
    begin = parse_nonnegative_int(begin, end, width);
    if (width != -1)
      handler.on_width(width);
    else
      throw STD_FMT::format_error("number is too big");
  } else if (*begin == '{') {
    ++begin;
    if (begin != end) begin = parse_arg_id(begin, end, width_adapter{handler});
    if (begin == end || *begin != '}') throw STD_FMT::format_error("invalid format string");
    ++begin;
  }
  return begin;
}

template<std::input_iterator It, std::sentinel_for<It> S, typename Handler>
[[nodiscard]] constexpr It parse_precision(It begin, S end, Handler&& handler)
{
  struct precision_adapter {
    Handler& handler;
    constexpr void operator()() { handler.on_dynamic_precision(auto_id{}); }
    constexpr void operator()(size_t id) { handler.on_dynamic_precision(id); }
  };

  ++begin;
  auto c = begin != end ? *begin : std::iter_value_t<It>();
  if ('0' <= c && c <= '9') {
    auto precision = 0;
    begin = parse_nonnegative_int(begin, end, precision);
    if (precision != -1)
      handler.on_precision(precision);
    else
      throw STD_FMT::format_error("number is too big");
  } else if (c == '{') {
    ++begin;
    if (begin != end) begin = parse_arg_id(begin, end, precision_adapter{handler});
    if (begin == end || *begin++ != '}') throw STD_FMT::format_error("invalid format string");
  } else {
    throw STD_FMT::format_error("missing precision specifier");
  }
  return begin;
}

template<std::input_iterator It>
constexpr int code_point_length(It begin)
{
  if constexpr (sizeof(std::iter_value_t<It>) != 1) return 1;
  constexpr char lengths[] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                              0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 3, 3, 4, 0};
  int len = lengths[static_cast<unsigned char>(*begin) >> 3];

  // Compute the pointer to the next character early so that the next
  // iteration can start working on the next character. Neither Clang
  // nor GCC figure out this reordering on their own.
  return len + !len;
}

// Parses fill and alignment.
template<std::input_iterator It, std::sentinel_for<It> S, typename Handler>
[[nodiscard]] constexpr It parse_align(It begin, S end, Handler&& handler)
{
  gsl_Expects(begin != end);
  auto align = fmt_align::none;
  auto p = begin + code_point_length(begin);
  if (p >= end) p = begin;
  for (;;) {
    switch (to_ascii(*p)) {
      case '<':
        align = fmt_align::left;
        break;
      case '>':
        align = fmt_align::right;
        break;
      case '^':
        align = fmt_align::center;
        break;
      default:
        break;
    }
    if (align != fmt_align::none) {
      if (p != begin) {
        auto c = *begin;
        if (c == '{') throw STD_FMT::format_error("invalid fill character '{'");
        handler.on_fill(std::basic_string_view<std::iter_value_t<It>>(begin, static_cast<size_t>(p - begin)));
        begin = p + 1;
      } else
        ++begin;
      handler.on_align(align);
      break;
    } else if (p == begin) {
      break;
    }
    p = begin;
  }
  return begin;
}

// Parses standard format specifiers and sends notifications about parsed
// components to handler.
template<std::input_iterator It, std::sentinel_for<It> S, typename SpecHandler>
[[nodiscard]] constexpr It parse_format_specs(It begin, S end, SpecHandler&& handler)
{
  if (begin + 1 < end && begin[1] == '}' && is_ascii_letter(*begin) && *begin != 'L') {
    handler.on_type(*begin++);
    return begin;
  }

  if (begin == end) return begin;

  begin = ::units::detail::parse_align(begin, end, handler);
  if (begin == end) return begin;

  // Parse sign.
  begin = ::units::detail::parse_sign(begin, end, handler);
  if (begin == end) return begin;

  if (*begin == '#') {
    handler.on_hash();
    if (++begin == end) return begin;
  }

  // Parse zero flag.
  if (*begin == '0') {
    handler.on_zero();
    if (++begin == end) return begin;
  }

  begin = ::units::detail::parse_width(begin, end, handler);
  if (begin == end) return begin;

  // Parse precision.
  if (*begin == '.') {
    begin = ::units::detail::parse_precision(begin, end, handler);
    if (begin == end) return begin;
  }

  if (*begin == 'L') {
    handler.on_localized();
    ++begin;
  }

  // Parse type.
  if (begin != end && *begin != '}') handler.on_type(*begin++);
  return begin;
}

// A format specifier handler that sets fields in basic_format_specs.
template<typename Char>
class specs_setter {
protected:
  basic_format_specs<Char>& specs_;
public:
  constexpr explicit specs_setter(basic_format_specs<Char>& specs) : specs_(specs) {}
  constexpr void on_align(fmt_align align) { specs_.align = align; }
  constexpr void on_fill(std::basic_string_view<Char> fill) { specs_.fill = fill; }
  constexpr void on_sign(fmt_sign s) { specs_.sign = s; }
  constexpr void on_hash() { specs_.alt = true; }
  constexpr void on_localized() { specs_.localized = true; }
  constexpr void on_zero() { specs_.fill[0] = Char('0'); }
  constexpr void on_width(int width) { specs_.width = width; }
  constexpr void on_precision(int precision) { specs_.precision = precision; }
  constexpr void on_type(Char type) { specs_.type = static_cast<char>(type); }
};

// Format spec handler that saves references to arguments representing dynamic
// width and precision to be resolved at formatting time.
template<typename ParseContext>
class dynamic_specs_handler : public specs_setter<typename ParseContext::char_type> {
public:
  using char_type = TYPENAME ParseContext::char_type;

  constexpr dynamic_specs_handler(dynamic_format_specs<char_type>& specs, ParseContext& ctx) :
      specs_setter<char_type>(specs), specs_(specs), context_(ctx)
  {}

  template<typename T>
  constexpr void on_dynamic_width(T t)
  {
    specs_.dynamic_width_index = on_dynamic_arg(t, context_);
  }

  template<typename T>
  constexpr void on_dynamic_precision(T t)
  {
    specs_.dynamic_precision_index = on_dynamic_arg(t, context_);
  }
private:
  dynamic_format_specs<char_type>& specs_;
  ParseContext& context_;
};

}  // namespace units::detail

max_stars_repo_name: hofbi/units

max_issues_repo_name: hofbi/units

max_forks_repo_name: hofbi/units


text:
#include <stdio.h>
#include <string.h>
#include <errno.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

struct color {
	double value[3];
};

struct color cartesian_to_polar(struct color source_color) {
	struct color target_color;
	target_color.value[0] = source_color.value[0];
	target_color.value[1] = hypot(source_color.value[1], source_color.value[2]);
	target_color.value[2] = atan2(source_color.value[2], source_color.value[1])*180.0/M_PI;

	if (target_color.value[2] < 0) {
		target_color.value[2] = target_color.value[2] + 360;
	}

	return target_color;
}

struct color polar_to_cartesian(struct color source_color) {
	struct color target_color;
	double angle;
	angle = source_color.value[2]/180*(M_PI);
	target_color.value[0] = source_color.value[0];
	target_color.value[1] = cos(angle)*source_color.value[1];
	target_color.value[2] = sin(angle)*source_color.value[1];
	return target_color;
}

/* XYZ = CIE 1931 XYZ */

struct color XYZ_values_to_OSA74(double X, double Y, double Z) {
	struct color target_color;
	double x = X/(X+Y+Z);
	double y = Y/(X+Y+Z);
	double Y0 = (4.4934*x*x + 4.3034*y*y - 4.276*x*y - 1.3744*x - 2.5643*y + 1.8103)*Y;
	double Lambda = 5.9*(cbrt(Y0) - 2.0/3 + 0.042*cbrt(Y0-30));
	double C = Lambda/(5.9*(cbrt(Y0) - 2.0/3));
	double R = 0.7990*X + 0.4194*Y - 0.1648*Z;
	double G = -0.4493*X + 1.3265*Y + 0.0927*Z;
	double B = -0.1149*X + 0.3394*Y + 0.7170*Z;
	/* target_color.value[0] = (Lambda - 14.3993)/sqrt(2); */
	target_color.value[0] = (Lambda - 14.4)/sqrt(2);
	target_color.value[1] = C*(1.7*cbrt(R) + 8*cbrt(G) - 9.7*cbrt(B));
	target_color.value[2] = C*(-13.7*cbrt(R) + 17.7*cbrt(G) - 4*cbrt(B));
	return target_color;
}

struct color XYZ_values_to_OSA90(double X, double Y, double Z) {
	struct color target_color;
	double x = X/(X+Y+Z);
	double y = Y/(X+Y+Z);
	double Y0 = (4.4934*x*x + 4.3034*y*y - 4.276*x*y - 1.3744*x - 2.5643*y + 1.8103)*Y;
	double Lambda = 5.9*(cbrt(Y0) - 2.0/3 + 0.06*cbrt(Y0-30));
	double C = Lambda/(5.9*(cbrt(Y0) - 2.0/3));
	double R = 0.9285*X + 0.3251*Y - 0.1915*Z;
	double G = -0.4493*X + 1.3265*Y + 0.0927*Z;
	double B = -0.2032*X + 0.6*Y + 0.5523*Z;
	/* target_color.value[0] = (Lambda - 14.3993)/sqrt(2); */
	target_color.value[0] = (Lambda - 14.4)/sqrt(2);
	target_color.value[1] = C*(-1.3*cbrt(R) + 17*cbrt(G) - 15.7*cbrt(B));
	target_color.value[2] = C*(-12.7*cbrt(R) + 19*cbrt(G) - 6.3*cbrt(B));
	return target_color;
}

struct color XYZ_to_OSA74(struct color source_color) {
	double X, Y, Z;
	X = source_color.value[0];
	Y = source_color.value[1];
	Z = source_color.value[2];
	return XYZ_values_to_OSA74(X, Y, Z);
}

struct color XYZ_to_OSA90(struct color source_color) {
	double X, Y, Z;
	X = source_color.value[0];
	Y = source_color.value[1];
	Z = source_color.value[2];
	return XYZ_values_to_OSA90(X, Y, Z);
}

int print_matrix(gsl_matrix *matrix) {
	int row, column;
	for(row = 0; row < matrix->size1; row++) {
		for(column = 0; column < matrix->size2; column++) {
			printf("%8.8f", gsl_matrix_get(matrix, row, column));
		}
		putchar('\n');
	}
	return 0;
}

struct color OSA74_to_XYZ(struct color source_color) {
	double L, j, g;
	L = source_color.value[0];
	j = source_color.value[1];
	g = source_color.value[2];
	struct color target_color;
	int max_iterations = 600;
	double X, Y, Z, previous_X, previous_Y, previous_Z, Delta;
	double L0, L1, L2, L3, j0, j1, j2, j3, g0, g1, g2, g3;
	struct color color0, color1, color2, color3;

	gsl_matrix *matrix0 = gsl_matrix_alloc(3, 3);
	gsl_matrix *inverse_matrix0 = gsl_matrix_alloc(3, 3);
	gsl_matrix *matrix1 = gsl_matrix_alloc(3, 1);
	gsl_matrix *matrix2 = gsl_matrix_alloc(3, 1);
	gsl_permutation *permutation = gsl_permutation_alloc(3);
	int sign = 0;

	/* start values */
	/*X = 28.4; Y = 30.0; Z = 32.2; Delta = 0.5;*/
	X = 28.4; Y = 30.0; Z = 32.2; Delta = 0.005;
	gsl_matrix_set(matrix2, 0, 0, X);
	gsl_matrix_set(matrix2, 1, 0, Y);
	gsl_matrix_set(matrix2, 2, 0, Z);

	int i = 0; do {
		/* calculate colors for influence matrix */
		color0 = XYZ_values_to_OSA74(X, Y, Z);
		color1 = XYZ_values_to_OSA74(X+Delta, Y, Z);
		color2 = XYZ_values_to_OSA74(X, Y+Delta, Z);
		color3 = XYZ_values_to_OSA74(X, Y, Z+Delta);
		L0 = color0.value[0]; j0 = color0.value[1]; g0 = color0.value[2];
		L1 = color1.value[0]; j1 = color1.value[1]; g1 = color1.value[2];
		L2 = color2.value[0]; j2 = color2.value[1]; g2 = color2.value[2];
		L3 = color3.value[0]; j3 = color3.value[1]; g3 = color3.value[2];

		gsl_matrix_set(matrix0, 0, 0, (L1-L0)/Delta);
		gsl_matrix_set(matrix0, 0, 1, (L2-L0)/Delta);
		gsl_matrix_set(matrix0, 0, 2, (L3-L0)/Delta);
		gsl_matrix_set(matrix0, 1, 0, (j1-j0)/Delta);
		gsl_matrix_set(matrix0, 1, 1, (j2-j0)/Delta);
		gsl_matrix_set(matrix0, 1, 2, (j3-j0)/Delta);
		gsl_matrix_set(matrix0, 2, 0, (g1-g0)/Delta);
		gsl_matrix_set(matrix0, 2, 1, (g2-g0)/Delta);
		gsl_matrix_set(matrix0, 2, 2, (g3-g0)/Delta);

		gsl_matrix_set(matrix1, 0, 0, L-L0);
		gsl_matrix_set(matrix1, 1, 0, j-j0);
		gsl_matrix_set(matrix1, 2, 0, g-g0);

		/* update matrix2:  matrix2 = inverse_matrix0  matrix1 + matrix2 */
		/* incremental term (inverse_matrix0 × matrix1) is divided by (5+cbrt(i)) */
		gsl_linalg_LU_decomp(matrix0, permutation, &sign);
		gsl_linalg_LU_invert(matrix0, permutation, inverse_matrix0);
		gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1/(5+cbrt(i)), inverse_matrix0, matrix1, 1.0, matrix2);

		previous_X = X; previous_Y = Y; previous_Z = Z;
		X = gsl_matrix_get(matrix2, 0, 0);
		Y = gsl_matrix_get(matrix2, 1, 0);
		Z = gsl_matrix_get(matrix2, 2, 0);

		/* printf("\n%d\n", i);
		print_matrix(matrix2); */

		i++;
	} while ((i < max_iterations) && ((fabs(X - previous_X) > 0.0001) || (fabs(Y - previous_Y) > 0.0001) || (fabs(Z - previous_Z) > 0.0001)));

	/*printf("Iterations: %i\n", i);*/

	/*free memory*/
	gsl_matrix_free(matrix0);
	gsl_matrix_free(inverse_matrix0);
	gsl_matrix_free(matrix1);
	gsl_matrix_free(matrix2);
	/*gsl_matrix_free(matrix_product);*/
	gsl_permutation_free(permutation);

	if (i == max_iterations) {
		fprintf(stderr, "Error: Maximum iterations reached.\n");
	}

	target_color.value[0] = X;
	target_color.value[1] = Y;
	target_color.value[2] = Z;
	return target_color;
}

struct color OSA90_to_XYZ(struct color source_color) {
	double L, j, g;
	L = source_color.value[0];
	j = source_color.value[1];
	g = source_color.value[2];
	struct color target_color;
	int max_iterations = 600;
	double X, Y, Z, previous_X, previous_Y, previous_Z, Delta;
	double L0, L1, L2, L3, j0, j1, j2, j3, g0, g1, g2, g3;
	struct color color0, color1, color2, color3;

	gsl_matrix *matrix0 = gsl_matrix_alloc(3, 3);
	gsl_matrix *inverse_matrix0 = gsl_matrix_alloc(3, 3);
	gsl_matrix *matrix1 = gsl_matrix_alloc(3, 1);
	gsl_matrix *matrix2 = gsl_matrix_alloc(3, 1);
	gsl_permutation *permutation = gsl_permutation_alloc(3);
	int sign = 0;

	/* start values */
	/*X = 28.4; Y = 30.0; Z = 32.2; Delta = 0.5;*/
	X = 28.4; Y = 30.0; Z = 32.2; Delta = 0.005;
	gsl_matrix_set(matrix2, 0, 0, X);
	gsl_matrix_set(matrix2, 1, 0, Y);
	gsl_matrix_set(matrix2, 2, 0, Z);

	int i = 0; do {
		/* calculate colors for influence matrix */
		color0 = XYZ_values_to_OSA90(X, Y, Z);
		color1 = XYZ_values_to_OSA90(X+Delta, Y, Z);
		color2 = XYZ_values_to_OSA90(X, Y+Delta, Z);
		color3 = XYZ_values_to_OSA90(X, Y, Z+Delta);
		L0 = color0.value[0]; j0 = color0.value[1]; g0 = color0.value[2];
		L1 = color1.value[0]; j1 = color1.value[1]; g1 = color1.value[2];
		L2 = color2.value[0]; j2 = color2.value[1]; g2 = color2.value[2];
		L3 = color3.value[0]; j3 = color3.value[1]; g3 = color3.value[2];

		gsl_matrix_set(matrix0, 0, 0, (L1-L0)/Delta);
		gsl_matrix_set(matrix0, 0, 1, (L2-L0)/Delta);
		gsl_matrix_set(matrix0, 0, 2, (L3-L0)/Delta);
		gsl_matrix_set(matrix0, 1, 0, (j1-j0)/Delta);
		gsl_matrix_set(matrix0, 1, 1, (j2-j0)/Delta);
		gsl_matrix_set(matrix0, 1, 2, (j3-j0)/Delta);
		gsl_matrix_set(matrix0, 2, 0, (g1-g0)/Delta);
		gsl_matrix_set(matrix0, 2, 1, (g2-g0)/Delta);
		gsl_matrix_set(matrix0, 2, 2, (g3-g0)/Delta);

		gsl_matrix_set(matrix1, 0, 0, L-L0);
		gsl_matrix_set(matrix1, 1, 0, j-j0);
		gsl_matrix_set(matrix1, 2, 0, g-g0);

		/* update matrix2:  matrix2 = inverse_matrix0  matrix1 + matrix2 */
		/* incremental term (inverse_matrix0 × matrix1) is divided by (5+cbrt(i)) */
		gsl_linalg_LU_decomp(matrix0, permutation, &sign);
		gsl_linalg_LU_invert(matrix0, permutation, inverse_matrix0);
		gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1/(5+cbrt(i)), inverse_matrix0, matrix1, 1.0, matrix2);

		previous_X = X; previous_Y = Y; previous_Z = Z;
		X = gsl_matrix_get(matrix2, 0, 0);
		Y = gsl_matrix_get(matrix2, 1, 0);
		Z = gsl_matrix_get(matrix2, 2, 0);

		/* printf("\n%d\n", i);
		print_matrix(matrix2); */

		i++;
	} while ((i < max_iterations) && ((fabs(X - previous_X) > 0.0001) || (fabs(Y - previous_Y) > 0.0001) || (fabs(Z - previous_Z) > 0.0001)));

	/*printf("Iterations: %i\n", i);*/

	/*free memory*/
	gsl_matrix_free(matrix0);
	gsl_matrix_free(inverse_matrix0);
	gsl_matrix_free(matrix1);
	gsl_matrix_free(matrix2);
	/*gsl_matrix_free(matrix_product);*/
	gsl_permutation_free(permutation);

	if (i == max_iterations) {
		fprintf(stderr, "Error: Maximum iterations reached.\n");
	}

	target_color.value[0] = X;
	target_color.value[1] = Y;
	target_color.value[2] = Z;
	return target_color;
}

double npow(double base, double exponent) {
	if (base < 0) {
		return -pow(fabs(base), exponent);
	} else {
		return pow(base, exponent);
	}
}

const double c1 = 3424/pow(2, 12);
const double c2 = 2413/pow(2, 7);
const double c3 = 2392/pow(2, 7);
const double n = 2610/pow(2, 14);
const double p = 1.7*2523/pow(2, 5);

double perceptual_quantizer(double cone_response) {
	return npow((c1 + c2*npow(cone_response/10000, n)) / (1 + c3*npow(cone_response/10000, n)), p);
}

double invert_perceptual_quantizer(double cone_response) {
	/*printf("input: %f\n", cone_response);*/
	/*printf("return: %f\n", 10000*npow((c1 - npow(cone_response, 1/p))/(c3*npow(cone_response, 1/p) - c2), 1/n));*/

	return 10000*npow((c1 - npow(cone_response, 1/p))/(c3*npow(cone_response, 1/p) - c2), 1/n);
}

const double d = -0.56;
const double d0 = 1.6295499532821566E-11;

struct color XYZ_to_Jab(struct color source_color) {
	double X = source_color.value[0];
	double Y = source_color.value[1];
	double Z = source_color.value[2];
	double adjusted_X, adjusted_Y, L, M, S, adjusted_L, adjusted_M, adjusted_S, I, J, a, b, C, h;
	struct color target_color;

	adjusted_X = 1.15*X - 0.15*Z;
	adjusted_Y = 0.66*Y + 0.34*X;

	L = 0.41478972*adjusted_X + 0.579999*adjusted_Y + 0.0146480*Z;
	M = -0.2015100*adjusted_X + 1.120649*adjusted_Y + 0.0531008*Z;
	S = -0.0166008*adjusted_X + 0.264800*adjusted_Y + 0.6684799*Z;

	adjusted_L = perceptual_quantizer(L);
	adjusted_M = perceptual_quantizer(M);
	adjusted_S = perceptual_quantizer(S);

	I = 0.5*adjusted_L + 0.5*adjusted_M;
	J = ((1+d)*I)/(1+d*I) - d0;

	a = 3.524000*adjusted_L - 4.066708*adjusted_M + 0.542708*adjusted_S;
	b = 0.199076*adjusted_L + 1.096799*adjusted_M - 1.295875*adjusted_S;

	target_color.value[0] = J;
	target_color.value[1] = a;
	target_color.value[2] = b;
	return target_color;
}

struct color Jab_to_XYZ(struct color source_color) {
	double J = source_color.value[0];
	double a = source_color.value[1];
	double b = source_color.value[2];
	double I, adjusted_L, adjusted_M, adjusted_S, L, M, S, adjusted_X, adjusted_Y, X, Y, Z;
	struct color target_color;

	I = ((J + d0)/(1 + d - d*(J + d0)));

	adjusted_L = I + 0.138605*a + 0.0580473*b;
	adjusted_M = I - 0.138605*a - 0.0580473*b;
	adjusted_S = I - 0.0960192*a - 0.811892*b;

	L = invert_perceptual_quantizer(adjusted_L);
	M = invert_perceptual_quantizer(adjusted_M);
	S = invert_perceptual_quantizer(adjusted_S);

	adjusted_X = 1.92423*L - 1.00479*M + 0.0376514*S;
	adjusted_Y = 0.350317*L + 0.726481*M - 0.0653844*S;
	Z = -0.0909828*L - 0.312728*M + 1.52277*S;

	X = (adjusted_X + 0.15*Z)/1.15;
	Y = (adjusted_Y - 0.34*X)/0.66;

	target_color.value[0] = X;
	target_color.value[1] = Y;
	target_color.value[2] = Z;
	return target_color;
}

const double J_scale = 0.167174631034780;

struct color scale_Jab(struct color source_color) {
	struct color target_color;

	target_color.value[0] = source_color.value[0]/J_scale*100;
	target_color.value[1] = source_color.value[1]/J_scale*100;
	target_color.value[2] = source_color.value[2];
	return target_color;
}

struct color unscale_Jab(struct color source_color) {
	struct color target_color;

	target_color.value[0] = source_color.value[0]*J_scale/100;
	target_color.value[1] = source_color.value[1]*J_scale/100;
	target_color.value[2] = source_color.value[2];
	return target_color;
}

double apply_gamma_correction(double value) {
	if (value > 0.0031308) {
		return 1.055*pow(value, (1.0/2.4)) - 0.055;
	} else {
		return 12.92*value;
	}
}

double remove_gamma_correction(double value) {
	if (value > 0.04045) {
		return pow(((value+0.055)/1.055), 2.4);
	} else {
		return value/12.92;
	}
}

struct color XYZ_to_sRGB(struct color source_color) {
	double X, Y, Z;
	X = source_color.value[0];
	Y = source_color.value[1];
	Z = source_color.value[2];
	struct color target_color;
	double linear_R, linear_G, linear_B, R, G, B;

	/* XYZ source values range = [0, 100] */
	X = X/100.0; Y = Y/100.0; Z = Z/100.0;

	/* 2° observer, D65 illuminant */
	linear_R = X* 3.2406 + Y*-1.5372 + Z*-0.4986;
	linear_G = X*-0.9689 + Y* 1.8758 + Z* 0.0415;
	linear_B = X* 0.0557 + Y*-0.2040 + Z* 1.0570;

	R = apply_gamma_correction(linear_R);
	G = apply_gamma_correction(linear_G);
	B = apply_gamma_correction(linear_B);

	if ((R<0.0) || (R>1.0) || (G<0.0) || (G>1.0) || (B<0.0) || (B>1.0)) {
		errno = ERANGE;
	}

	/* RGB target values range = [0, 255] */
	target_color.value[0] = R*255;
	target_color.value[1] = G*255;
	target_color.value[2] = B*255;

	return target_color;
}

struct color sRGB_to_XYZ(struct color source_color) {
	double R, G, B;
	R = source_color.value[0];
	G = source_color.value[1];
	B = source_color.value[2];
	struct color target_color;
	double linear_R, linear_G, linear_B, X, Y, Z;

	/* RGB source values range = [0, 255] */
	R = R/255.0; G = G/255.0; B = B/255.0;

	linear_R = remove_gamma_correction(R);
	linear_G = remove_gamma_correction(G);
	linear_B = remove_gamma_correction(B);

	/* 2° observer, D65 illuminant */
	X = linear_R*0.4124 + linear_G*0.3576 + linear_B*0.1805;
	Y = linear_R*0.2126 + linear_G*0.7152 + linear_B*0.0722;
	Z = linear_R*0.0193 + linear_G*0.1192 + linear_B*0.9505;

	/* XYZ target values range = [0, 100] */
	target_color.value[0] = 100*X;
	target_color.value[1] = 100*Y;
	target_color.value[2] = 100*Z;
	return target_color;
}

const double Y_threshold = pow((6/29.0), 3);
const double D65_X = 95.0156;
const double D65_Y = 100;
const double D65_Z = 108.8199;

/*const double D65_X = 97.31273;*/
/*const double D65_Y = 100.00;*/
/*const double D65_Z = 138.59590;*/

struct color XYZ_to_Luv(struct color source_color) {
	/* L, u, v stand for L*, u*, v* */
	double X, Y, Z, prime_u, prime_v, relative_Y, white_prime_u, white_prime_v, L, u, v;
	struct color target_color;

	X = source_color.value[0];
	Y = source_color.value[1];
	Z = source_color.value[2];

	prime_u = (4*X) / (X + 15*Y + 3*Z);
	prime_v = (9*Y) / (X + 15*Y + 3*Z);

	white_prime_u = (4*D65_X) / (D65_X + 15*D65_Y + 3*D65_Z);
	white_prime_v = (9*D65_Y) / (D65_X + 15*D65_Y + 3*D65_Z);

	relative_Y = Y/D65_Y;

	if (relative_Y > pow((6/29.0), 3)) {
		L = cbrt(relative_Y)*116 - 16;
	} else {
		L = relative_Y*pow((29/3.0), 3);
	}

	u = 13*L*(prime_u - white_prime_u);
	v = 13*L*(prime_v - white_prime_v);

	target_color.value[0] = L;
	target_color.value[1] = u;
	target_color.value[2] = v;
	return target_color;
}

struct color Luv_to_XYZ(struct color source_color) {
	/* L, u, v stand for L*, u*, v* */
	double L, u, v, white_prime_u, white_prime_v, prime_u, prime_v, X, Y, Z;
	struct color target_color;

	L = source_color.value[0];
	u = source_color.value[1];
	v = source_color.value[2];

	white_prime_u = (4*D65_X) / (D65_X + 15*D65_Y + 3*D65_Z);
	white_prime_v = (9*D65_Y) / (D65_X + 15*D65_Y + 3*D65_Z);

	prime_u = u/(13*L) + white_prime_u;
	prime_v = v/(13*L) + white_prime_v;

	if (L > 8.0) {
		Y = D65_Y*pow((L + 16)/116.0, 3);
	} else {
		Y = D65_Y*L*pow(3/29.0, 3);
	}

	X = Y * (9*prime_u)/(4.0*prime_v);
	Z = Y * (12 - 3*prime_u - 20*prime_v)/(4.0*prime_v);

	target_color.value[0] = X;
	target_color.value[1] = Y;
	target_color.value[2] = Z;
	return target_color;
}

/* OSA74 range: -19.1032--10.0938 (29.1970) */

struct color OSA74_to_pOSA74(struct color source_color) {
	double L, j, g, new_L, chroma, hue;
	struct color target_color;

	L = source_color.value[0];
	j = source_color.value[1];
	g = source_color.value[2];

	new_L = (L*sqrt(2) + 19.1032)*100/29.1970;
	chroma = hypot(j, g)*100/29.1970;
	hue = atan2(g, j)*180.0/M_PI;

	if (hue < 0) {
		hue = hue + 360;
	}

	target_color.value[0] = new_L;
	target_color.value[1] = chroma;
	target_color.value[2] = hue;
	return target_color;
}

struct color pOSA74_to_OSA74(struct color source_color) {
	double new_L, chroma, hue, L, angle, j, g;
	struct color target_color;

	new_L = source_color.value[0];
	chroma = source_color.value[1];
	hue = source_color.value[2];

	L = (new_L*29.1970/100 - 19.1032)/sqrt(2);

	angle = hue/180*(M_PI);
	j = cos(angle)*chroma*29.1970/100;
	g = sin(angle)*chroma*29.1970/100;

	target_color.value[0] = L;
	target_color.value[1] = j;
	target_color.value[2] = g;
	return target_color;
}

int RGB_to_hex(struct color RGB_color) {
	int R = (int)(RGB_color.value[0] + 0.5);
	int G = (int)(RGB_color.value[1] + 0.5);
	int B = (int)(RGB_color.value[2] + 0.5);
	return (R<<16) | (G<<8) | B;
}

struct color hex_to_RGB(int hex_value) {
	struct color RGB_color;
	RGB_color.value[0] = ((hex_value>>16) & 0xFF)/255.0;
	RGB_color.value[1] = ((hex_value>>8) & 0xFF)/255.0;
	RGB_color.value[2] = (hex_value & 0xFF)/255.0;
	return RGB_color;
}

int main(int argc, char *argv[]) {
	int i;
	struct color source_color;
	struct color intermediate_color;
	struct color target_color;
	char next_char;
	char *source_format = argv[1];
	char *target_format = argv[2];

	/* check arguments */

	/* check argument count */
	if (argc != 6) {
		fprintf(stderr, "Usage: %s source_format target_format value1 value2 value3\n", argv[0]);
		return EXIT_FAILURE;
	}

	if (strcmp(source_format, target_format) == 0) {
		fprintf(stderr, "Source and target formats are identical.\n");
		return EXIT_FAILURE;
	}

	/* check if color values are numbers */
	for (i = 3; i < argc; i++) {
		if (sscanf(argv[i], "%lf%c", &source_color.value[i-3], &next_char) != 1) {
			fprintf(stderr, "Invalid color values, values must be numbers.\n");
			return EXIT_FAILURE;
		}
	}

	/* OSA/polar conversions do not need intermediate conversion to XYZ */
	if ((strcmp(source_format, "pOSA74") == 0) && (strcmp(target_format, "OSA74") == 0)) {
		target_color = pOSA74_to_OSA74(source_color);
		printf("%.2f %.2f %.2f\n", target_color.value[0], target_color.value[1], target_color.value[2]);
		return EXIT_SUCCESS;
	} else if ((strcmp(source_format, "OSA74") == 0) && (strcmp(target_format, "pOSA74") == 0)) {
		target_color = OSA74_to_pOSA74(source_color);
		printf("%.2f %.2f %.2f\n", target_color.value[0], target_color.value[1], target_color.value[2]);
		return EXIT_SUCCESS;
	}

	/* convert from source format to XYZ */
	if (strcmp(source_format, "Jab") == 0) {
		intermediate_color = Jab_to_XYZ(polar_to_cartesian(unscale_Jab(source_color)));
		/*intermediate_color = Jab_to_XYZ(polar_to_cartesian(source_color));*/
		/*intermediate_color = Jab_to_XYZ(source_color);*/
	} else if (strcmp(source_format, "LCh") == 0) {
		intermediate_color = Luv_to_XYZ(polar_to_cartesian(source_color));
	} else if (strcmp(source_format, "Luv") == 0) {
		intermediate_color = Luv_to_XYZ(source_color);
	} else if (strcmp(source_format, "OSA74") == 0) {
		intermediate_color = OSA74_to_XYZ(source_color);
	} else if (strcmp(source_format, "OSA90") == 0) {
		intermediate_color = OSA90_to_XYZ(source_color);
	} else if (strcmp(source_format, "pOSA74") == 0) {
		/*intermediate_color = OSA90_to_XYZ(polar_to_OSA(source_color));*/
		intermediate_color = OSA74_to_XYZ(pOSA74_to_OSA74(source_color));
	} else if (strcmp(source_format, "sRGB") == 0) {
		intermediate_color = sRGB_to_XYZ(source_color);
	} else if (strcmp(source_format, "XYZ") == 0) {
		intermediate_color = source_color;
	} else {
		fprintf(stderr, "Invalid source format.");
		return EXIT_FAILURE;
	}

	/* convert from XYZ to target format */
	if (strcmp(target_format, "Jab") == 0) {
		target_color = scale_Jab(cartesian_to_polar(XYZ_to_Jab(intermediate_color)));
		/*target_color = cartesian_to_polar(XYZ_to_Jab(intermediate_color));*/
		/*target_color = XYZ_to_Jab(intermediate_color);*/
	} else if (strcmp(target_format, "LCh") == 0) {
		target_color = cartesian_to_polar(XYZ_to_Luv(intermediate_color));
	} else if (strcmp(target_format, "Luv") == 0) {
		target_color = XYZ_to_Luv(intermediate_color);
	} else if (strcmp(target_format, "OSA74") == 0) {
		target_color = XYZ_to_OSA74(intermediate_color);
	} else if (strcmp(target_format, "OSA90") == 0) {
		target_color = XYZ_to_OSA90(intermediate_color);
	} else if (strcmp(target_format, "pOSA74") == 0) {
		target_color = OSA74_to_pOSA74(XYZ_to_OSA74(intermediate_color));
	} else if (strcmp(target_format, "polarOSA74") == 0) {
		target_color = cartesian_to_polar(XYZ_to_OSA74(intermediate_color));
	} else if (strcmp(target_format, "sRGB") == 0) {
		target_color = XYZ_to_sRGB(intermediate_color);
	} else if (strcmp(target_format, "XYZ") == 0) {
		target_color = intermediate_color;
	} else {
		fprintf(stderr, "Invalid target format.");
		return EXIT_FAILURE;
	}
	if (errno == ERANGE) {
		fprintf(stderr, "The requested color lies outside of the %s gamut.\n", target_format);
		/*return EXIT_FAILURE;*/
		/*printf("%d %d %d\n", 149, 149, 148);*/
		return EXIT_SUCCESS;
	}

	/* round RGB values to integers */
	if (strcmp(target_format, "sRGB") == 0) {
		/*printf("%.0f %.0f %.0f\n", target_color.value[0], target_color.value[1], target_color.value[2]);*/
		printf("#%06x\n", RGB_to_hex(target_color));
	} else {
		printf("%.4f %.4f %.4f\n", target_color.value[0], target_color.value[1], target_color.value[2]);
	}

	return EXIT_SUCCESS;
}

max_stars_repo_name: joshuakraemer/color-converter

max_issues_repo_name: joshuakraemer/color-converter

max_forks_repo_name: joshuakraemer/color-converter


text:
/******************************************************************************
 * gsl_sprng.h: rewrite of gsl-sprng.h to use sprng streams                   *
 * Author: Yan Y. Liu <yanliu@illinois.edu>                                   *
 * Date: 2014/08/17                                                           *
 * Copyright and license of this source file are specified in LICENSE.TXT     *
 * under the root directory of this software package.                         *
 ******************************************************************************/

/*
 * gsl-sprng.h: rewrite of gsl-sprng.h to use sprng streams
 *
 * This header file is based on  Darren Wilkinson's version for simple interface,
 * which does not use streams. We also refer to IceCube's C++ extension based on
 * sprng 4.0.
 *
 * Darren Wilkinson  http://www.staff.ncl.ac.uk/d.j.wilkinson/
 * IceCube project: http://icecube.wisc.edu/
 *
 * To use, just add the line:
 * #include "gsl-sprng.h"
 * immediately after the line:
 * #include <gsl/gsl_rng.h>
 * near the start of your code.
 * Make sure you alloc the rng on each processor. If you wish to
 * set a seed, you should set it to be the same on each processor.
 */

#ifndef GSL_SPRNG_H
#define GSL_SPRNG_H

#define USE_MPI
#include "sprng.h"

/***************************/
/* gsl definition of sprng */
/***************************/

typedef struct
{
    int streamnum;
    int nstreams;
    int seed;
    int *stream;
} gsl_sprng_state_t;

/* seed is gsl_rng_default_seed
 * either 0 by default or GSL_RNG_SEED environment variable
 */
void sprng_set(void * vstate, unsigned long int seed);
/* the get() function to get integer random number */
unsigned long sprng_get(void * vstate);
/* the get() function to get double random number */
double sprng_get_double(void * vstate);

/***************************/
/*      gsl_sprng API      */
/***************************/
/* replace all gsl_rng_* func with corresponding one below */

/* initialize gsl_rng */
gsl_rng *gsl_sprng_alloc (int streamnum, int nstreams, int seed, int * stream);
/* finalize gsl_rng */
void gsl_sprng_free (gsl_rng * r);
/* print out gsl_rng stream info for debug purpose */
void gsl_sprng_print (gsl_rng * r);
/* load gsl_rng from a file */
int gsl_sprng_fread (char * fpath, gsl_rng * r);
/* store gsl_rng to a file */
int gsl_sprng_fwrite (char * fpath, const gsl_rng * r);
#endif


text:
// https://github.com/attractivechaos/matmul
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <stdio.h>

/**********************************
 * Pseudo-random number generator *
 **********************************/

static uint64_t mat_rng[2] = { 11ULL, 1181783497276652981ULL };

static inline uint64_t xorshift128plus(uint64_t s[2])
{
	uint64_t x, y;
	x = s[0], y = s[1];
	s[0] = y;
	x ^= x << 23;
	s[1] = x ^ y ^ (x >> 17) ^ (y >> 26);
	y += s[1];
	return y;
}

double mat_drand(void)
{
	return (xorshift128plus(mat_rng)>>11) * (1.0/9007199254740992.0);
}

/*******************************************
 * Helper routines for matrix manipulation *
 *******************************************/

float **mat_init(int n_rows, int n_cols)
{
	float **m;
	int i;
	m = (float**)malloc(n_rows * sizeof(float*));
	m[0] = (float*)calloc(n_rows * n_cols, sizeof(float));
	for (i = 1; i < n_rows; ++i)
		m[i] = m[i-1] + n_cols;
	return m;
}

void mat_destroy(float **m)
{
	free(m[0]); free(m);
}

float **mat_gen_random(int n_rows, int n_cols)
{
	float **m;
	int i, j;
	m = mat_init(n_rows, n_cols);
	for (i = 0; i < n_rows; ++i)
		for (j = 0; j < n_cols; ++j)
			m[i][j] = mat_drand();
	return m;
}

float **mat_transpose(int n_rows, int n_cols, float *const* a)
{
	int i, j;
	float **m;
	m = mat_init(n_cols, n_rows);
	for (i = 0; i < n_rows; ++i)
		for (j = 0; j < n_cols; ++j)
			m[j][i] = a[i][j];
	return m;
}

float sdot_1(int n, const float *x, const float *y)
{
	int i;
	float s = 0.0f;
	for (i = 0; i < n; ++i) s += x[i] * y[i];
	return s;
}

float sdot_8(int n, const float *x, const float *y)
{
	int i, n8 = n>>3<<3;
	float s, t[8];
	t[0] = t[1] = t[2] = t[3] = t[4] = t[5] = t[6] = t[7] = 0.0f;
	for (i = 0; i < n8; i += 8) {
		t[0] += x[i+0] * y[i+0];
		t[1] += x[i+1] * y[i+1];
		t[2] += x[i+2] * y[i+2];
		t[3] += x[i+3] * y[i+3];
		t[4] += x[i+4] * y[i+4];
		t[5] += x[i+5] * y[i+5];
		t[6] += x[i+6] * y[i+6];
		t[7] += x[i+7] * y[i+7];
	}
	for (s = 0.0f; i < n; ++i) s += x[i] * y[i];
	s += t[0] + t[1] + t[2] + t[3] + t[4] + t[5] + t[6] + t[7];
	return s;
}

#ifdef __SSE__
#include <xmmintrin.h>

float sdot_sse(int n, const float *x, const float *y)
{
	int i, n8 = n>>3<<3;
	__m128 vs1, vs2;
	float s, t[4];
	vs1 = _mm_setzero_ps();
	vs2 = _mm_setzero_ps();
	for (i = 0; i < n8; i += 8) {
		__m128 vx1, vx2, vy1, vy2;
		vx1 = _mm_loadu_ps(&x[i]);
		vx2 = _mm_loadu_ps(&x[i+4]);
		vy1 = _mm_loadu_ps(&y[i]);
		vy2 = _mm_loadu_ps(&y[i+4]);
		vs1 = _mm_add_ps(vs1, _mm_mul_ps(vx1, vy1));
		vs2 = _mm_add_ps(vs2, _mm_mul_ps(vx2, vy2));
	}
	for (s = 0.0f; i < n; ++i) s += x[i] * y[i];
	_mm_storeu_ps(t, vs1);
	s += t[0] + t[1] + t[2] + t[3];
	_mm_storeu_ps(t, vs2);
	s += t[0] + t[1] + t[2] + t[3];
	return s;
}
#endif

/*************************
 * Matrix multiplication *
 *************************/

float **mat_mul0(int n_a_rows, int n_a_cols, float *const *a, int n_b_cols, float *const *b)
{
	int i, j, k;
	float **m;
	m = mat_init(n_a_rows, n_b_cols);
	for (i = 0; i < n_a_rows; ++i) {
		for (j = 0; j < n_b_cols; ++j) {
			float t = 0.0;
			for (k = 0; k < n_a_cols; ++k)
				t += a[i][k] * b[k][j];
			m[i][j] = t;
		}
	}
	return m;
}

float **mat_mul1(int n_a_rows, int n_a_cols, float *const *a, int n_b_cols, float *const *b)
{
	int i, j, k, n_b_rows = n_a_cols;
	float **m, **bT;
	m = mat_init(n_a_rows, n_b_cols);
	bT = mat_transpose(n_b_rows, n_b_cols, b);
	for (i = 0; i < n_a_rows; ++i) {
		const float *ai = a[i];
		float *mi = m[i];
		for (j = 0; j < n_b_cols; ++j) {
			float t = 0.0f, *bTj = bT[j];
			for (k = 0; k < n_a_cols; ++k)
				t += ai[k] * bTj[k];
			mi[j] = t;
		}
	}
	mat_destroy(bT);
	return m;
}

#ifdef __SSE__
float **mat_mul2(int n_a_rows, int n_a_cols, float *const *a, int n_b_cols, float *const *b)
{
	int i, j, n_b_rows = n_a_cols;
	float **m, **bT;
	m = mat_init(n_a_rows, n_b_cols);
	bT = mat_transpose(n_b_rows, n_b_cols, b);
	for (i = 0; i < n_a_rows; ++i)
		for (j = 0; j < n_b_cols; ++j)
			m[i][j] = sdot_sse(n_a_cols, a[i], bT[j]);
	mat_destroy(bT);
	return m;
}
float **mat_mul7(int n_a_rows, int n_a_cols, float *const *a, int n_b_cols, float *const *b)
{
	int i, j, ii, jj, x = 16, n_b_rows = n_a_cols;
	float **m, **bT;
	m = mat_init(n_a_rows, n_b_cols);
	bT = mat_transpose(n_b_rows, n_b_cols, b);
	for (i = 0; i < n_a_rows; i += x) {
		for (j = 0; j < n_b_cols; j += x) {
			int je = n_b_cols < j + x? n_b_cols : j + x;
			int ie = n_a_rows < i + x? n_a_rows : i + x;
			for (ii = i; ii < ie; ++ii)
				for (jj = j; jj < je; ++jj)
					m[ii][jj] += sdot_sse(n_a_cols, a[ii], bT[jj]);
		}
	}
	mat_destroy(bT);
	return m;
}
#endif

float **mat_mul3(int n_a_rows, int n_a_cols, float *const *a, int n_b_cols, float *const *b)
{
	int i, j, n_b_rows = n_a_cols;
	float **m, **bT;
	m = mat_init(n_a_rows, n_b_cols);
	bT = mat_transpose(n_b_rows, n_b_cols, b);
	for (i = 0; i < n_a_rows; ++i)
		for (j = 0; j < n_b_cols; ++j)
			m[i][j] = sdot_8(n_a_cols, a[i], bT[j]);
	mat_destroy(bT);
	return m;
}

float **mat_mul4(int n_a_rows, int n_a_cols, float *const *a, int n_b_cols, float *const *b)
{
	int i, j, n_b_rows = n_a_cols;
	float **m, **bT;
	m = mat_init(n_a_rows, n_b_cols);
	bT = mat_transpose(n_b_rows, n_b_cols, b);
	for (i = 0; i < n_a_rows; ++i)
		for (j = 0; j < n_b_cols; ++j)
			m[i][j] = sdot_1(n_a_cols, a[i], bT[j]);
	mat_destroy(bT);
	return m;
}

#ifdef HAVE_CBLAS
#include <cblas.h>

float **mat_mul5(int n_a_rows, int n_a_cols, float *const *a, int n_b_cols, float *const *b)
{
	int i, j, n_b_rows = n_a_cols;
	float **m, **bT;
	m = mat_init(n_a_rows, n_b_cols);
	bT = mat_transpose(n_b_rows, n_b_cols, b);
	for (i = 0; i < n_a_rows; ++i)
		for (j = 0; j < n_b_cols; ++j)
			m[i][j] = cblas_sdot(n_a_cols, a[i], 1, bT[j], 1);
	mat_destroy(bT);
	return m;
}

float **mat_mul6(int n_a_rows, int n_a_cols, float *const *a, int n_b_cols, float *const *b)
{
	float **m, n_b_rows = n_a_cols;
	m = mat_init(n_a_rows, n_b_cols);
	cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n_a_rows, n_b_cols, n_a_cols, 1.0f, a[0], n_a_rows, b[0], n_b_rows, 0.0f, m[0], n_a_rows);
	return m;
}
#endif

/*****************
 * Main function *
 *****************/

#include <unistd.h>
#include <time.h>

int main(int argc, char *argv[])
{
	int c, n = 1000, algo = 2;
	clock_t t;
	float **a, **b, **m = 0;

	while ((c = getopt(argc, argv, "n:a:h")) >= 0) {
		if (c == 'n') n = atoi(optarg);
		else if (c == 'a') algo = atoi(optarg);
		else if (c == 'h') {
			fprintf(stderr, "Usage: mat-eval [options]\n");
			fprintf(stderr, "Options:\n");
			fprintf(stderr, "  -n INT    size of the square matrix [%d]\n", n);
			fprintf(stderr, "  -a INT    matrix multiplication implementation [%d]\n", algo);
			fprintf(stderr, "            0: naive - no optimization\n");
			fprintf(stderr, "            1: transposing the second matrix\n");
#ifdef __SSE__
			fprintf(stderr, "            2: explicitly vectorized sdot() with SSE\n");
			fprintf(stderr, "            7: explicitly SSE sdot() plus loop tiling\n");
#endif
			fprintf(stderr, "            3: implicitly vectorized sdot()\n");
			fprintf(stderr, "            4: no vectorization hints\n");
#ifdef HAVE_CBLAS
			fprintf(stderr, "            5: with sdot() from an external CBLAS library\n");
			fprintf(stderr, "            6: with sgemm() from an external CBLAS library\n");
#endif
			fprintf(stderr, "  -h        this help message\n");
			return 1;
		}
	}

	a = mat_gen_random(n, n);
	b = mat_gen_random(n, n);

	t = clock();
	if (algo == 0) {
		m = mat_mul0(n, n, a, n, b);
	} else if (algo == 1) {
		m = mat_mul1(n, n, a, n, b);
#ifdef __SSE__
	} else if (algo == 2) {
		m = mat_mul2(n, n, a, n, b);
	} else if (algo == 7) {
		m = mat_mul7(n, n, a, n, b);
#endif
	} else if (algo == 3) {
		m = mat_mul3(n, n, a, n, b);
	} else if (algo == 4) {
		m = mat_mul4(n, n, a, n, b);
#ifdef HAVE_CBLAS
	} else if (algo == 5) {
		m = mat_mul5(n, n, a, n, b);
	} else if (algo == 6) {
		m = mat_mul6(n, n, a, n, b);
#endif
	} else {
		fprintf(stderr, "ERROR: unknown algorithm %d\n", algo);
		return 1;
	}
	fprintf(stderr, "CPU time: %g\n", (double)(clock() - t) / CLOCKS_PER_SEC);
	fprintf(stderr, "Central cell: %g\n", m[n/2][n/2]);

	if (m) mat_destroy(m);
	mat_destroy(b); mat_destroy(a);
	return 0;
}

max_stars_repo_name: chrisroman/llvm-pass-skeleton

max_issues_repo_name: chrisroman/llvm-pass-skeleton

max_forks_repo_name: chrisroman/llvm-pass-skeleton