software/scripts/gradient.sh

Go to the documentation of this file.

#!/usr/bin/env zsh
script=${0##*/}
 
zmodload zsh/mathfunc
 
if [[ ! $ZSH_EVAL_CONTEXT =~ :file$ ]]; then
    echo "Auxiliary script for running the conjugate gradient fit in zsh."
    echo "This script should be sourced by the steering script."
    exit
fi
 
# This script contains auxiliary functions for running the conjugate gradient fit.
# It should be sourced by the steering script.
# The steering script should provide for the method getChi2 and set the fit parameters before calling method gradient.
# On return, the method getChi2 should set the current value of chi2 to its first positional parameter.
# The second positional parameter corresponds to an option that can be used within method getChi2.
# The value of the option corresponds to the following cases.
#  0 => step wise improvement of the chi2;
#  1 => evaluation of the chi2 before the determination of the gradient of the chi2; and
#  2 => evaluation of the derivative of the chi2 to each fit parameter.
# The fit parameters are defined by the associative array PARAMETERS.
# In this, each key corresponds to a valid zsh command that can be used for changing a specific fit parameter.
# The command should contain the wild card '%' which is repeatedly replaced by a value, referred to as step size.
# The associated value is used to evaluate the derivative of chi2 per unit step size of each fit parameter.
# These step sizes should be tuned beforehand to yield similar absolute values of the derivatives.
 
typeset -a    CHI2                                # internal chi2 values; index [1] corresponds to best value
typeset -a    GRADIENT                            # derivative of chi2 per unit step size of each fit parameter
let          "EPSILON               =  5.0E-4"    # default maximal distance to minimum
let          "NUMBER_OF_ITERATIONS  =  1000"      # default maximum number of iterations
typeset -A    PARAMETERS                          # fit parameters
 
if [[ -z "$DEBUG" ]]; then
    DEBUG=0
fi
 
 
#
# Auxiliary function to evaluate chi2 and gradient thereof at current position.
#
# The external method getChi2 is used to evaluate the chi2.
# The value corresponding to each key in the associative array PARAMETERS 
# is used as unit step size for the evaluation of the derivative of the chi2.
# The results are stored in the internal variables CHI2[1] and GRADIENT, respectively.
#
function evaluate()
{
    GRADIENT=()
 
    getChi2 CHI2\[1\] 1
 
    let "WIDTH = 1"
 
    for KEY in ${(k)PARAMETERS}; do
        if (( ${#KEY} > $WIDTH )); then
            let "WIDTH = ${#KEY}"
        fi
    done
 
    let "V  = 0.0"
    let "i  = 1"
 
    for KEY in ${(ko)PARAMETERS}; do
 
        BUFFER=(${(s/;/)${(Q)KEY}})
 
        VALUE=$PARAMETERS[$KEY]
 
        if (( $VALUE != 0.0 )); then
 
            for K1 in ${(o)BUFFER}; do
                eval ${(Q)${(qq)K1//\%/$((+0.5 * $VALUE))}}
            done
 
            getChi2 CHI2\[2\] 2
 
            for K1 in ${(O)BUFFER}; do
                eval ${(Q)${(qq)K1//\%/$((-0.5 * $VALUE))}}
            done
 
            for K1 in ${(o)BUFFER}; do
                eval ${(Q)${(qq)K1//\%/$((-0.5 * $VALUE))}}
            done
 
            getChi2 CHI2\[3\] 2
 
            for K1 in ${(O)BUFFER}; do
                eval ${(Q)${(qq)K1//\%/$((+0.5 * $VALUE))}}
            done
 
            GRADIENT[$i]=$(($CHI2[2] - $CHI2[3]))
        else
 
            GRADIENT[$i]=0.0
        fi
 
        if (( $DEBUG >= 3 )); then
            printf "%-${WIDTH}s %12.5f %12.5f\n" $KEY $PARAMETERS[$KEY] $GRADIENT[$i]
        fi
 
        let "V += $GRADIENT[$i]*$GRADIENT[$i]"
        let "i += 1"
    done
 
    if (( $DEBUG >= 3 )); then
        printf "%-${WIDTH}s %12s %12.5f\n" "|gradient|" "" $((sqrt($V)))
    fi
}
 
 
#
# Move current position along gradient by given scaling factor of the unit step sizes.
#
# \param 1     scaling factor
#
function move()
{
    if (( $1 > 0.0 )); then
 
        let "i  = 1"
 
        for KEY in ${(ko)PARAMETERS}; do
 
            BUFFER=(${(s/;/)${(Q)KEY}})
 
            for K1 in ${(o)BUFFER}; do
                eval ${(Q)${(qq)K1//\%/$(($GRADIENT[$i] * $1 * $PARAMETERS[$KEY]))}}
            done
 
            let "i += 1"
        done
 
    else
 
        let "i  = ${#GRADIENT}"
 
        for KEY in ${(kO)PARAMETERS}; do
 
            BUFFER=(${(s/;/)${(Q)KEY}})
 
            for K1 in ${(O)BUFFER}; do
                eval ${(Q)${(qq)K1//\%/$(($GRADIENT[$i] * $1 * $PARAMETERS[$KEY]))}}
            done
 
            let "i -= 1"
        done
    fi
}
 
 
#
# Main fit function.
#
# The algorithm is based on the conjugate gradient method.
# The position is moved to the optimal value and the final chi2 is returned.
#
# If the option is one, the step sizes are fine-tuned according to the gradient of the chi2.
# The value corresponding to each key in the associative array PARAMETERS is then 
# multiplied by the corresponding element in GRADIENT and only the sign thereof is maintained.
#
# \param 1     maximum number of extra steps
# \param 2     variable containing final chi2 value on return
#
function gradient()
{
    typeset -a G
    typeset -a H
 
    if (( ${#PARAMETERS} == 0 )); then
        return
    fi
 
    evaluate
 
    # initialise
 
    for (( i = 1; i <= ${#GRADIENT}; ++i)); do
        G[$i]=$((-1.0 * $GRADIENT[$i]))
        H[$i]=$G[$i]
        GRADIENT[$i]=$H[$i]
    done
 
    for (( N = 0; $N != $NUMBER_OF_ITERATIONS; N += 1 )); do
 
        if (( $DEBUG >= 3 )); then
            printf "chi2[1]  %4d %12.5f\n" $N $CHI2[1]
        fi
 
        # minimise chi2 in direction of gradient
 
        CHI2[2]=$CHI2[1]
        CHI2[3]=$CHI2[2]
 
        for (( DS = 1.0, M = 0 ; $DS > 1.0e-3; )); do
 
            move $DS
 
            getChi2 CHI2\[4\] 0
 
            if (( $DEBUG >= 3 )); then
                printf "chi2[4]  %4d %12.5f %12.5e\n" $M $CHI2[4] $DS
            fi
 
            if (( $CHI2[4] < $CHI2[3] )); then
 
                # continue
 
                CHI2[2]=$CHI2[3]
                CHI2[3]=$CHI2[4]
 
                let "M   = 0"
 
                continue
            fi
 
            # if chi2 increases, try additional steps first to overcome possible hurdle, then try single step with reduced size
 
            if (( $DS == 1.0 )); then
 
                if (( $M == 0 )); then
                    CHI2[5]=$CHI2[4]
                fi
 
                if (( $M != $1 )); then
 
                    let "M  += 1"
 
                    continue;
 
                else
 
                    for (( ; $M != 0; --M )); do                            
                        move $((-1.0 * $DS))
                    done
 
                    CHI2[4]=$CHI2[5]
                fi
            fi
 
            move $((-1.0 * $DS))
 
            if (( $CHI2[3] < $CHI2[2] )); then
 
                # final step based on parabolic interpolation through following points
                #
                # X1 = -1 * DS  ->  CHI2[2]
                # X2 =  0 * DS  ->  CHI2[3]
                # X3 = +1 * DS  ->  CHI2[4]
 
                let "F21 = $CHI2[3] - $CHI2[2]"   # F(X2) - F(X1)
                let "F23 = $CHI2[3] - $CHI2[4]"   # F(X2) - F(X3)
 
                let "XS  = 0.5 * ($F21 - $F23) / ($F23 + $F21)"
 
                move $((+1.0 * $XS * $DS))
 
                getChi2 CHI2\[4\] 0
 
                if (( $CHI2[4] < $CHI2[3] )); then
 
                    CHI2[3]=$CHI2[4]
 
                else
 
                    move $((-1.0 * $XS * $DS))
 
                    getChi2 CHI2\[3\] 0
                fi
 
                if (( $DEBUG >= 3 )); then
                    printf "chi2[3]  %4d %12.5f %12.5e\n" $N $CHI2[3] $DS
                fi
 
                break
 
            else
 
                # reduce step size
                
                let "DS *= 0.5"
            fi
        done
 
        if (( fabs($CHI2[3] - $CHI2[1]) < $EPSILON * 0.5 * (fabs($CHI2[1]) + fabs($CHI2[3])) )); then
 
            CHI2[1]=$CHI2[3]
 
            break;
        fi
 
        evaluate
 
        let "GG   = 0.0"
        let "DGG  = 0.0"
 
        for (( i = 1; $i <= ${#GRADIENT}; ++i )); do
            let "GG  +=  $G[$i]*$G[$i]"
            let "DGG += ($GRADIENT[$i] + $G[$i]) * $GRADIENT[$i]"
        done
 
        if (( $GG == 0.0 )); then
            break
        fi
 
        let "GAM  = $DGG / $GG"
        
        for (( i = 1; $i <= ${#GRADIENT}; ++i )); do
            G[$i]=$((-1.0 * $GRADIENT[$i]))
            H[$i]=$(($G[$i] + $GAM * $H[$i]))
            GRADIENT[$i]=$H[$i]
        done
    done
 
    if (( $DEBUG >= 3 )); then
        printf "chi2[1]  %4d %12.5f\n" $N $CHI2[1]
    fi
 
    eval $2=$CHI2[1]
}
 
 
#
# Additional fit function to be called after gradient.
#
# The algorithm is based on sorting the derivatives of the chi2 and 
# moving the position one-by-one to the optimal value.
# The final chi2 is returned.
#
# \param 1     variable containing final chi2 value on return
#
function simplex()
{
    typeset -A G
 
 
    let "WIDTH = 1"
 
    for KEY in ${(k)PARAMETERS}; do
        if (( ${#KEY} > $WIDTH )); then
            let "WIDTH = ${#KEY}"
        fi
    done
 
 
    for (( N = 0; $N != $NUMBER_OF_ITERATIONS; N += 1 )); do
 
        G=()
 
        GRADIENT=()
 
        getChi2 CHI2\[1\] 1
 
        let "i  = 1"
 
        for KEY in ${(ko)PARAMETERS}; do
 
            VALUE=$PARAMETERS[$KEY]
 
            BUFFER=(${(s/;/)${(Q)KEY}})
 
            for K1 in ${(o)BUFFER}; do
                eval ${(Q)${(qq)K1//\%/$((+1.0 * $VALUE))}}
            done
 
            getChi2 CHI2\[2\] 2
 
            for K1 in ${(O)BUFFER}; do
                eval ${(Q)${(qq)K1//\%/$((-1.0 * $VALUE))}}
            done
 
            for K1 in ${(o)BUFFER}; do
                eval ${(Q)${(qq)K1//\%/$((-1.0 * $VALUE))}}
            done
 
            getChi2 CHI2\[3\] 2
 
            for K1 in ${(O)BUFFER}; do
                eval ${(Q)${(qq)K1//\%/$((+1.0 * $VALUE))}}
            done
 
            if (( $DEBUG >= 3 )); then
                printf "%-${WIDTH}s %12.5f %12.5f -> %+9.5f <- %+9.5f\n" $KEY $VALUE $CHI2[1] $(($CHI2[2] - $CHI2[1])) $(($CHI2[3] - $CHI2[1]))
            fi
 
            # map index to derivative of chi2 and maintain sign of step
 
            if   (($CHI2[2] < $CHI2[1] && $CHI2[2] < $CHI2[3])); then
 
                G[$i]=$(($CHI2[2] - $CHI2[1]))
 
                GRADIENT[$i]=+1.0
 
            elif (($CHI2[3] < $CHI2[1] && $CHI2[3] < $CHI2[2])); then
 
                G[$i]=$(($CHI2[3] - $CHI2[1]))
 
                GRADIENT[$i]=-1.0
#
#           elif (($CHI2[3] - $CHI2[2] > 0.5*($CHI2[2] + $CHI2[3]) - $CHI2[1])); then
#
#               G[$i]=$(($CHI2[2] - $CHI2[1]))
#
#               GRADIENT[$i]=+0.33
#
#           elif (($CHI2[2] - $CHI2[3] > 0.5*($CHI2[2] + $CHI2[3]) - $CHI2[1])); then
#
#               G[$i]=$(($CHI2[3] - $CHI2[1]))
#
#               GRADIENT[$i]=-0.33
            fi
 
            let "i += 1"
        done
 
 
        if (( ${#G} == 0 )); then
            break
        fi
 
 
        # sort derivatives of chi2 and optimise position one-by-one
 
        CHI2[2]=$CHI2[1]
 
        for i VALUE in `printf "%d %12.3e\n" ${(kv)G} | sort -g -k2`; do
 
            KEY=${${(ko)PARAMETERS}[$i]}
 
            BUFFER=(${(s/;/)${(Q)KEY}})
 
            # initial step size
            
            let "DS  = $GRADIENT[$i] * $PARAMETERS[$KEY]"
            
            if (($DEBUG >= 3)); then
                printf "%-${WIDTH}s %12.5f %12.5f\n" $KEY $DS $G[$i]
            fi
 
            CHI2[3]=$CHI2[2]
 
            for (( n = 1; ; ++n )); do
 
                for K1 in ${(o)BUFFER}; do
                    eval ${(Q)${(qq)K1//\%/$((+1.0 * $DS))}}
                done
 
                getChi2 CHI2\[4\] 0
 
                if (( $DEBUG >= 3 )); then
                    printf "chi2[4]  %4d %12.5f\n" $n $CHI2[4]
                fi
 
                if (( $CHI2[4] < $CHI2[3] )); then
 
                    CHI2[3]=$CHI2[4]
 
                    continue
                fi
 
                for K1 in ${(O)BUFFER}; do
                    eval ${(Q)${(qq)K1//\%/$((-1.0 * $DS))}}
                done
 
                if (( $CHI2[3] < $CHI2[2] )); then
 
                    # final step based on parabolic interpolation through following points
                    #
                    # X1 = -1 * DS  ->  CHI2[2]
                    # X2 =  0 * DS  ->  CHI2[3]
                    # X3 = +1 * DS  ->  CHI2[4]
 
                    let "F21 = $CHI2[3] - $CHI2[2]"   # F(X2) - F(X1)
                    let "F23 = $CHI2[3] - $CHI2[4]"   # F(X2) - F(X3)
 
                    let "XS  = 0.5 * ($F21 - $F23) / ($F23 + $F21)"
 
                    for K1 in ${(o)BUFFER}; do
                        eval ${(Q)${(qq)K1//\%/$((+1.0 * $XS * $DS))}}
                    done
                    
                    getChi2 CHI2\[4\] 0
 
                    if (( $CHI2[4] < $CHI2[3] )); then
 
                        CHI2[3]=$CHI2[4]
 
                    else
 
                        for K1 in ${(O)BUFFER}; do
                            eval ${(Q)${(qq)K1//\%/$((-1.0 * $XS * $DS))}}
                        done
                    fi
                fi
 
                if (( $DEBUG >= 3 )); then
                    printf "chi2[3]  %4d %12.5f\n" $n $CHI2[3]
                fi
 
                break
            done
 
            CHI2[2]=$CHI2[3]
        done
 
        if (( fabs($CHI2[2] - $CHI2[1]) < $EPSILON * 0.5 * (fabs($CHI2[1]) + fabs($CHI2[2])) )); then
 
            CHI2[1]=$CHI2[2]
 
            break;
        fi
    done
 
    if (( $DEBUG >= 3 )); then
        printf "chi2[1]  %4d %12.5f\n" $N $CHI2[1]
    fi
 
    eval $1=$CHI2[1]
}