software/scripts/gradient.sh

Go to the documentation of this file.

#!/bin/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.
# The method getChi2 should set the current value of chi2 to its first positional 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.
# To this end, the methods evaluate and gprint (in this order) can be used.

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\]

    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\]

            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\]

            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\]

            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  ->  CHI[2]
                # X2 =  0 * DS  ->  CHI[3]
                # X3 = +1 * DS  ->  CHI[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\]

                if (( $CHI2[4] < $CHI2[3] )); then

                    CHI2[3]=$CHI2[4]

                else

                    move $((-1.0 * $XS * $DS))

                    getChi2 CHI2\[3\]
                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\]

        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\]

            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\]

            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\]

                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  ->  CHI[2]
                    # X2 =  0 * DS  ->  CHI[3]
                    # X3 = +1 * DS  ->  CHI[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\]

                    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]
}