package ICC::Support::bern;
use strict;
use Carp;
our $VERSION = 0.41;
# revised 2018-09-04
#
# Copyright © 2004-2020 by William B. Birkett
# inherit from Shared
use parent qw(ICC::Shared);
# enable static variables
use feature 'state';
# create new bern object
# hash keys are: 'input', 'output', 'fit', 'fix_hl', 'fix_sh'
# default values for 'input' and 'output' are []
# parameters: ([ref_to_attribute_hash])
# returns: (ref_to_object)
sub new {
        # get object class
        my $class = shift();
        # create empty bern object
        my $self = [
                {},     # object header
                [],     # parameter array
                [],     # min/max t-values
                []      # min/max i/o-values
        ];
        # if one parameter, a hash reference
        if (@_ == 1 && ref($_[0]) eq 'HASH') {
                 
                # set object contents from hash
                _new_from_hash($self, $_[0]);
                 
        } elsif (@_) {
                 
                # error
                croak('invalid parameter(s)');
                 
        }
        # return blessed object
        return(bless($self, $class));
}
# create inverse object
# swaps the 'input' and 'output' arrays
# returns: (ref_to_object)
sub inv {
        # clone new object
        my $self = Storable::dclone(shift());
        # swap 'input' and 'output' arrays
        ($self->[1][0], $self->[1][1]) = ($self->[1][1], $self->[1][0]);
        ($self->[2][0], $self->[2][1]) = ($self->[2][1], $self->[2][0]);
        ($self->[3][0], $self->[3][1]) = ($self->[3][1], $self->[3][0]);
        # return inverted object
        return($self);
}
# write bern object to ICC profile
# note: writes an equivalent 'curv' object
# parameters: (ref_to_parent_object, file_handle, ref_to_tag_table_entry)
sub write_fh {
        # verify 4 parameters
        (@_ == 4) or croak('wrong number of parameters');
        # write bern data to profile
        goto &_writeICCbern;
}
# get tag size (for writing to profile)
# note: writes an equivalent curv object
# returns: (tag_size)
sub size {
        # get parameters
        my ($self) = @_;
        # get count value (defaults to 4370)
        my $n = $self->[0]{'curv_points'} // 4370;
        # return size
        return(12 + $n * 2);
}
# compute curve function
# parameters: (input_value)
# returns: (output_value)
sub transform {
        # get parameters
        my ($self, $in) = @_;
        # verify input
        (Scalar::Util::looks_like_number($in)) or croak('invalid transform input');
        # compute intermediate value(s)
        my @t = _rev($self->[1][0], $self->[2][0], $self->[3][0], $in);
        # warn if multiple values
        (@t == 1) or carp('multiple transform values');
        # return first output value
        return(_fwd($self->[1][1], $t[0]));
}
# compute inverse curve function
# parameters: (input_value)
# returns: (output_value)
sub inverse {
        # get parameters
        my ($self, $in) = @_;
        # verify input
        (Scalar::Util::looks_like_number($in)) or croak('invalid transform input');
        # compute intermediate value(s)
        my @t = _rev($self->[1][1], $self->[2][1], $self->[3][1], $in);
        # warn if multiple values
        (@t == 1) or carp('multiple transform values');
        # return first output value
        return(_fwd($self->[1][0], $t[0]));
}
# compute curve derivative
# parameters: (input_value)
# returns: (derivative_value)
sub derivative {
        # get parameters
        my ($self, $in) = @_;
        # return forward derivative
        return(_derivative($self, 0, $in));
}
# compute parametric partial derivatives
# parameters: (input_value)
# returns: (partial_derivative_array)
sub parametric {
        # get parameters
        my ($self, $in) = @_;
        # return forward parametric partial derivatives
        return(_parametric($self, 0, $in));
}
# get/set header hash reference
# parameters: ([ref_to_new_hash])
# returns: (ref_to_hash)
sub header {
        # get object reference
        my $self = shift();
        # if there are parameters
        if (@_) {
                 
                # if one parameter, a hash reference
                if (@_ == 1 && ref($_[0]) eq 'HASH') {
                         
                        # set header to copy of hash
                        $self->[0] = Storable::dclone(shift());
                         
                } else {
                         
                        # error
                        croak('parameter must be a hash reference');
                         
                }
                 
        }
        # return reference
        return($self->[0]);
}
# get/set input array reference
# an array of two values is a range
# otherwise, it contains input values
# parameters: ([ref_to_array])
# returns: (ref_to_array)
sub input {
        # get object reference
        my $self = shift();
        # if one parameter supplied
        if (@_ == 1) {
                 
                # verify 'input' vector (between 1 and 21 numeric values)
                (ICC::Shared::is_num_vector($_[0]) && 1 <= @{$_[0]} && 21 >= @{$_[0]}) or croak('invalid input values');
                 
                # store output values
                $self->[1][0] = [@{$_[0]}];
                 
                # check for min/max values
                _minmax($self, 0);
                 
        } elsif (@_) {
                 
                # error
                carp('too many parameters');
                 
        }
        # return array reference
        return($self->[1][0]);
}
# get/set output array reference
# parameters: ([ref_to_array])
# returns: (ref_to_array)
sub output {
        # get object reference
        my $self = shift();
        # if one parameter supplied
        if (@_ == 1) {
                 
                # verify 'output' vector (between 1 and 21 numeric values)
                (ICC::Shared::is_num_vector($_[0]) && 1 <= @{$_[0]} && 21 >= @{$_[0]}) or croak('invalid input values');
                 
                # store output values
                $self->[1][1] = [@{$_[0]}];
                 
                # check for min/max values
                _minmax($self, 1);
                 
        } elsif (@_) {
                 
                # error
                carp('too many parameters');
                 
        }
        # return array reference
        return($self->[1][1]);
}
# normalize transform
# adjusts Bernstein coefficients linearly
# adjusts 'input' and 'output' by default
# hash keys: 'input', 'output'
# parameter: ([hash_ref])
# returns: (ref_to_object)
sub normalize {
        # get parameters
        my ($self, $hash) = @_;
        # local variables
        my ($m, $b, $val, $src, $x0, $x1, $y0, $y1, $f0, $f1, @minmax);
        # set hash key array
        my @key = qw(input output);
        # for ('input', 'output')
        for my $i (0, 1) {
                 
                # undefine slope
                $m = undef;
                 
                # if no parameter (default)
                if (! defined($hash)) {
                         
                        # if array has 2 or more coefficients
                        if (@{$self->[1][$i]} > 1 && ($x0 = $self->[1][$i][0]) != ($x1 = $self->[1][$i][-1])) {
                                 
                                # compute adjustment values
                                $m = abs(1/($x0 - $x1));
                                $b = -$m * ($x0 < $x1 ? $x0 : $x1);
                                 
                                # set endpoints
                                ($y0, $y1) = $x0 < $x1 ? (0, 1) : (1, 0);
                                 
                        }
                         
                } else {
                         
                        # if hash contains key
                        if (defined($val = $hash->{$key[$i]})) {
                                 
                                # if array has 2 or more coefficients
                                if (@{$self->[1][$i]} > 1 && $self->[1][$i][0] != $self->[1][$i][-1]) {
                                         
                                        # if hash value is an array reference
                                        if (ref($val) eq 'ARRAY') {
                                                 
                                                # if two parameters
                                                if (@{$val} == 2) {
                                                         
                                                        # get parameters
                                                        ($y0, $y1) = @{$val};
                                                         
                                                        # get endpoints
                                                        $x0 = $self->[1][$i][0];
                                                        $x1 = $self->[1][$i][-1];
                                                         
                                                # if three parameters
                                                } elsif (@{$val} == 3) {
                                                         
                                                        # get parameters
                                                        ($src, $y0, $y1) = @{$val};
                                                         
                                                        # if source is 'endpoints'
                                                        if (! ref($src) && $src eq 'endpoints') {
                                                                 
                                                                # get endpoints
                                                                $x0 = $self->[1][$i][0];
                                                                $x1 = $self->[1][$i][-1];
                                                                 
                                                        # if source is 'minmax'
                                                        } elsif (! ref($src) && $src eq 'minmax') {
                                                                 
                                                                # get all possible min/max values
                                                                @minmax = ($self->[1][$i][0], $self->[1][$i][-1], @{$self->[3][$i]});
                                                                 
                                                                # get min/max values
                                                                $x0 = List::Util::min(@minmax);
                                                                $x1 = List::Util::max(@minmax);
                                                                 
                                                        }
                                         
                                                # if four parameters
                                                } elsif (@{$val} == 4) {
                                                         
                                                        # get parameters
                                                        ($x0, $x1, $y0, $y1) = @{$val};
                                                         
                                                }
                                                 
                                                # if x/y values are numeric
                                                if ((4 == grep {Scalar::Util::looks_like_number($_)} ($x0, $x1, $y0, $y1)) && ($x0 != $x1)) {
                                                         
                                                        # compute slope
                                                        $m = ($y0 - $y1)/($x0 - $x1);
                                                         
                                                        # compute offset
                                                        $b = $y0 - $m * $x0;
                                                         
                                                } else {
                                                         
                                                        # warning
                                                        carp("invalid '$key[$i]' parameter(s)");
                                                         
                                                }
                                         
                                        }
                                 
                                } else {
                                         
                                        # warning
                                        carp("can't normalize '$key[$i]' coefficients");
                                         
                                }
                                 
                        }
                 
                }
                 
                # if slope is defined
                if (defined($m)) {
                         
                        # set endpoint flags
                        $f0 = $self->[1][$i][0] == $x0;
                        $f1 = $self->[1][$i][-1] == $x1;
                         
                        # adjust Bernstein coefficients (y = mx + b)
                        @{$self->[1][$i]} = map {$m * $_ + $b} @{$self->[1][$i]};
                         
                        # set endpoints (to ensure exact values)
                        $self->[1][$i][0] = $y0 if ($f0);
                        $self->[1][$i][-1] = $y1 if ($f1);
                         
                }
                 
        }
        # re-calculate min/max values
        _minmax($self);
        # return
        return($self);
}
# check if monotonic
# returns min/max values
# returns t-values by default
# returns 'input' values if format is 'input'
# returns 'output' values if format is 'output'
# curve is monotonic if no min/max values
# parameters: ([format])
# returns: (array_of_values)
sub monotonic {
        # get object reference
        my ($self, $fmt) = @_;
        # local variables
        my (@s);
        # compute min/max values
        _minmax($self);
        # if format undefined
        if (! defined($fmt)) {
                 
                # return sorted t-values ('input' and 'output')
                return(sort {$a <=> $b} (@{$self->[2][0]}, @{$self->[2][1]}));
                 
        # if format 'input'
        } elsif ($fmt eq 'input') {
                 
                # return input values
                return(map {_fwd($self->[1][0], $_)} sort {$a <=> $b} (@{$self->[2][0]}, @{$self->[2][1]}));
                 
        # if format 'output'
        } elsif ($fmt eq 'output') {
                 
                # return output values
                return(map {_fwd($self->[1][1], $_)} sort {$a <=> $b} (@{$self->[2][0]}, @{$self->[2][1]}));
                 
        } else {
                 
                # error
                croak('unsupported format for min/max values');
                 
        }
         
}
# update internal object elements
# this method used primarily when optimizing
# call after changing 'input' or 'output' values
sub update {
        # recompute min/max values
        goto &_minmax;
}
# make table for profile objects
# assumes curve domain/range is (0 - 1)
# parameters: (number_of_table_entries, [direction])
# returns: (ref_to_table_array)
sub table {
        # get parameters
        my ($self, $n, $dir) = @_;
        # local variables
        my ($up, $table);
        # validate number of table entries
        ($n == int($n) && $n >= 2) or carp('invalid number of table entries');
        # purify direction flag
        $dir = ($dir) ? 1 : 0;
        # array upper index
        $up = $n - 1;
        # for each table entry
        for my $i (0 .. $up) {
                 
                # compute table value
                $table->[$i] = _transform($self, $dir, $i/$up);
                 
        }
        # return table reference
        return($table);
}
# make equivalent 'curv' object
# assumes curve domain/range is (0 - 1)
# parameters: (number_of_table_entries, [direction])
# returns: (ref_to_curv_object)
sub curv {
        # return 'curv' object reference
        return(ICC::Profile::curv->new(table(@_)));
}
# print object contents to string
# format is an array structure
# parameter: ([format])
# returns: (string)
sub sdump {
        # get parameters
        my ($self, $p) = @_;
        # local variables
        my ($s, $fmt, $fmt2);
        # resolve parameter to an array reference
        $p = defined($p) ? ref($p) eq 'ARRAY' ? $p : [$p] : [];
        # get format string
        $fmt = defined($p->[0]) && ! ref($p->[0]) ? $p->[0] : 'undef';
        # set string to object ID
        $s = sprintf("'%s' object, (0x%x)\n", ref($self), $self);
        # if input_parameter_array defined
        if (defined($self->[1][0])) {
                 
                # set format
                $fmt2 = join(', ', ('%.3f') x @{$self->[1][0]});
                 
                # append parameter string
                $s .= sprintf("  input parameters  [$fmt2]\n", @{$self->[1][0]});
                 
        }
        # if output_parameter_array defined
        if (defined($self->[1][1])) {
                 
                # set format
                $fmt2 = join(', ', ('%.3f') x @{$self->[1][1]});
                 
                # append parameter string
                $s .= sprintf("  output parameters [$fmt2]\n", @{$self->[1][1]});
                 
        }
        # return
        return($s);
}
# find min/max values within the interval (0 - 1)
# used by '_rev' to determine solution regions
# called by the 'update' method
# i/o_index: 0 - input, 1 - output, undef - both
# parameters: (ref_to_object, [i/o_index])
sub _minmax {
        # get parameters
        my ($self, $io) = @_;
        # local variables
        my (@s, $t, $d);
        my ($par, $fwd, $fwd2, $drv);
        # for input and/or output curves
        for my $i (defined($io) ? ($io ? 1 : 0) : (0, 1)) {
                 
                # initialize s-values array
                @s = ();
                 
                # get ref to Bernstein parameters
                $par = $self->[1][$i];
                 
                # if number of parameters > 2
                if (@{$par} > 2) {
                         
                        # compute forward differences
                        $fwd = [map {$par->[$_] - $par->[$_ - 1]} (1 .. $#{$par})];
                         
                        # if forward differences have different signs (+/-)
                        if (grep {($fwd->[0] < 0) ^ ($fwd->[$_] < 0)} (1 .. $#{$fwd})) {
                                 
                                # compute second forward differences
                                $fwd2 = [map {$fwd->[$_] - $fwd->[$_ - 1]} (1 .. $#{$fwd})];
                                 
                                # compute derivative values over range (0 - 1)
                                $drv = [map {$#{$par} * _fwd($fwd, $_/(4 * $#{$par}))} (0 .. 4 * $#{$par})];
                                 
                                # for each pair of derivatives
                                for my $j (1 .. $#{$drv}) {
                                         
                                        # if derivatives have different signs (+/-)
                                        if (($drv->[$j - 1] < 0) ^ ($drv->[$j] < 0)) {
                                                 
                                                # estimate root from the derivative values
                                                $t = ($j * $drv->[$j - 1] - ($j - 1) * $drv->[$j])/(($drv->[$j - 1] - $drv->[$j]) * 4 * $#{$par});
                                                 
                                                # loop until derivative approaches 0
                                                while (abs($d = $#{$par} * _fwd($fwd, $t)) > 1E-9) {
                                                         
                                                        # adjust root using Newton's method
                                                        $t -= $d/($#{$fwd} * $#{$par} * _fwd($fwd2, $t));
                                                         
                                                }
                                                 
                                                # push root on t-value array
                                                push(@s, $t) if ($t > 0 && $t < 1);
                                                 
                                        }
                                         
                                }
                                 
                        }
                         
                }
                 
                # warn if root values were found
                print "curve $i is not monotonic\n" if (@s && ! defined($self->[2][$i]));
                 
                # save root t-values
                $self->[2][$i] = [@s];
                 
                # save root i/o-values
                $self->[3][$i] = [map {_fwd($par, $_)} @s];
                 
        }
         
}
# combined transform
# nominal domain (0 - 1)
# parameters: (ref_to_object, direction, input_value)
# returns: (output_value)
sub _transform {
        # get parameters
        my ($self, $dir, $in) = @_;
        # local variables
        my ($p0, $p1, $p2, $p3, @t);
        # verify input
        (Scalar::Util::looks_like_number($in)) or croak('invalid transform input');
        # get parameter arrays
        $p0 = $self->[1][$dir];
        $p1 = $self->[2][$dir];
        $p2 = $self->[3][$dir];
        $p3 = $self->[1][1 - $dir];
        # compute intermediate value(s)
        @t = _rev($p0, $p1, $p2, $in);
        # warn if multiple values
        (@t == 1) or carp('multiple transform values');
        # return first output value
        return(_fwd($p3, $t[0]));
}
# combined derivative
# nominal domain (0 - 1)
# parameters: (ref_to_object, direction, input_value)
# returns: (derivative_value)
sub _derivative {
        # get parameters
        my ($self, $dir, $in) = @_;
        # local variables
        my ($p0, $p1, $p2, $p3, @t);
        my ($di, $do);
        # verify input
        (Scalar::Util::looks_like_number($in)) or croak('invalid transform input');
        # get parameter arrays
        $p0 = $self->[1][$dir];
        $p1 = $self->[2][$dir];
        $p2 = $self->[3][$dir];
        $p3 = $self->[1][1 - $dir];
        # compute intermediate value(s)
        @t = _rev($p0, $p1, $p2, $in);
        # warn if multiple values
        (@t == 1) || print "multiple transform values\n";
        # compute input derivative (using first value)
        ($di = _drv($p0, $t[0])) or croak('infinite derivative value');
        # compute output derivative (using first value)
        $do = _drv($p3, $t[0]);
        # return combined derivative
        return($do/$di);
}
# compute parametric partial derivatives
# direction: 0 - normal, 1 - inverse
# parameters: (ref_to_object, direction, input_value)
# returns: (partial_derivative_array)
sub _parametric {
        # get parameters
        my ($self, $dir, $in) = @_;
        # local variables
        my ($p0, $p1, $p2, $p3);
        my (@t, $t, $di, $do, $bfi, $bfo, @pj);
        # verify input
        (Scalar::Util::looks_like_number($in)) or croak('invalid transform input');
        # get parameter arrays
        $p0 = $self->[1][$dir];
        $p1 = $self->[2][$dir];
        $p2 = $self->[3][$dir];
        $p3 = $self->[1][1 - $dir];
        # compute intermediate value(s)
        @t = _rev($p0, $p1, $p2, $in);
        # warn if multiple values
        (@t == 1) || print "multiple transform values\n";
        # use first value
        $t = $t[0];
        # compute input derivative
        $di = _drv($p0, $t);
        # compute input Bernstein polynomials
        $bfi = @{$p0} ? _bernstein($#{$p0}, $t) : undef();
        # compute output derivative
        $do = _drv($p3, $t);
        # compute output Bernstein polynomials
        $bfo = @{$p3} ? _bernstein($#{$p3}, $t) : undef();
        # initialize array
        @pj = ();
        # if there are input parameters
        if (defined($bfi)) {
                 
                # for each parameter
                for my $i (0 .. $#{$bfi}) {
                         
                        # verify non-zero input derivative
                        ($di) or croak('infinite Jacobian element');
                         
                        # compute partial derivative
                        $pj[$i] = -$bfi->[$i] * $do/$di;
                         
                }
                 
        }
        # add output parameters (if any)
        push(@pj, @{$bfo}) if (defined($bfo));
        # return array of partial derivatives
        return(@pj);
}
# direct forward transform
# nominal domain (0 - 1)
# parameters: (parameter_array_reference, input_value)
# returns: (output_value)
sub _fwd {
        # get parameters
        my ($par, $in) = @_;
        # if no parameters
        if (! @{$par}) {
                 
                # return input
                return($in);
                 
        # one parameter
        } elsif (@{$par} == 1) {
                 
                # return constant
                return($par->[0]);
                 
        # two parameters
        } elsif (@{$par} == 2) {
                 
                # return linear interpolated value
                return((1 - $in) * $par->[0] + $in * $par->[1]);
                 
        # three or more parameters
        } else {
                 
                # check if ICC::Support::Lapack module is loaded
                state $lapack = defined($INC{'ICC/Support/Lapack.pm'});
                 
                # if input < 0
                if ($in < 0) {
                         
                        # return extrapolated value
                        return($par->[0] + $in * _drv($par, 0));
                         
                # if input > 1
                } elsif ($in > 1) {
                         
                        # return extrapolated value
                        return($par->[-1] + ($in - 1) * _drv($par, 1));
                         
                } else {
                         
                        # if Lapack module loaded
                        if ($lapack) {
                                 
                                # return Bernstein interpolated value
                                return(ICC::Support::Lapack::bernstein($in, $par));
                                 
                        } else {
                                 
                                # return Bernstein interpolated value
                                return(ICC::Shared::dotProduct(_bernstein($#{$par}, $in), $par));
                                 
                        }
                         
                }
                 
        }
         
}
# direct reverse transform
# nominal range (0 - 1)
# parameters: (parameter_array_reference, x-min/max_array_reference, y-min/max_array_reference, input_value)
# returns: (output_value -or- array_of_output_values)
sub _rev {
        # get parameters
        my ($par, $xminmax, $yminmax, $in) = @_;
        # local variables
        my (@xs, @ys, $slope, $ext, @sol);
        my ($p, $loop, $xlo, $xhi, $xval, $yval);
        # if no parameters
        if (! @{$par}) {
                 
                # return input
                return($in);
                 
        # one parameter
        } elsif (@{$par} == 1) {
                 
                # warn if input not equal to constant
                ($in == $par->[0]) or carp('curve input not equal to constant parameter');
                 
                # return input
                return($in);
                 
        # two parameters
        } elsif (@{$par} == 2) {
                 
                # return linear interpolated value
                return(($in - $par->[0])/($par->[1] - $par->[0]));
                 
        # three or more parameters
        } else {
                 
                # setup segment arrays
                @xs = (0, @{$xminmax}, 1);
                @ys = ($par->[0], @{$yminmax}, $par->[-1]);
                 
                # if slope(0) not 0
                if ($slope = _drv($par, 0)) {
                         
                        # compute extrapolation from 0
                        $ext = ($in - $par->[0])/$slope;
                         
                        # save if less than 0
                        push(@sol, $ext) if ($ext < 0);
                         
                }
                 
                # test first y-value
                push(@sol, $xs[0]) if ($ys[0] == $in);
                 
                # for each array segment
                for my $i (1 .. $#xs) {
                         
                        # if input value falls within segment
                        if (($in > $ys[$i - 1] && $in < $ys[$i]) || ($in < $ys[$i - 1] && $in > $ys[$i])) {
                                 
                                # get segment polarity
                                $p = $ys[$i] > $ys[$i - 1];
                                 
                                # initialize limits
                                $xlo = $xs[$i - 1];
                                $xhi = $xs[$i];
                                 
                                # initialize loop counter
                                $loop = 0;
                                 
                                # bisection method loop
                                while ($xhi - $xlo > 1E-3 && $loop++ < 15) {
                                         
                                        # compute midpoint
                                        $xval = ($xlo + $xhi)/2;
                                         
                                        # compute y-value
                                        $yval = _fwd($par, $xval);
                                         
                                        # if y-value > input -or- y-value < input
                                        if ($p ? ($yval > $in) : ($yval < $in)) {
                                                 
                                                # set high limit
                                                $xhi = $xval;
                                                 
                                        } else {
                                                 
                                                # set low limit
                                                $xlo = $xval;
                                                 
                                        }
                                         
                                }
                                 
                                # initialize loop counter
                                $loop = 0;
                                 
                                # Newton's method loop
                                while (abs($in - $yval) > 1E-12 && $loop++ < 5 && ($slope = _drv($par, $xval))) {
                                         
                                        # adjust x-value
                                        $xval += ($in - $yval)/$slope;
                                         
                                        # compute y-value
                                        $yval = _fwd($par, $xval);
                                         
                                }
                                 
                                # print warning if failed to converge
                                printf("'_rev' function of 'bern' object failed to converge with input of %.2f\n", $in) if (abs($in - $yval) > 1E-12);
                                 
                                # save result
                                push(@sol, $xval);
                                 
                        }
                         
                        # add segment y-value, if a solution
                        push(@sol, $xs[$i]) if ($ys[$i] == $in);
                         
                }
                 
                # if slope(1) not 0
                if ($slope = _drv($par, 1)) {
                         
                        # compute extrapolation from 1
                        $ext = ($in - $par->[-1])/$slope + 1;
                         
                        # save if greater than 1
                        push(@sol, $ext) if ($ext > 1);
                         
                }
                 
        }
        # return result (array or first solution)
        return(wantarray ? @sol : $sol[0]);
}
# forward derivative
# nominal domain is (0 - 1)
# parameters: (parameter_array_reference, input_value)
# returns: (output_value)
sub _drv {
        # get parameters
        my ($par, $in) = @_;
        # if no parameters
        if (! @{$par}) {
                 
                # return identity derivative
                return(1);
                 
        # one parameter
        } elsif (@{$par} == 1) {
                 
                # return constant derivative
                return(0);
                 
        # two parameters
        } elsif (@{$par} == 2) {
                 
                # return linear derivative
                return($par->[1] - $par->[0]);
                 
        # three or more parameters
        } else {
                 
                # check if ICC::Support::Lapack module is loaded
                state $lapack = defined($INC{'ICC/Support/Lapack.pm'});
                 
                # limit input value (0 - 1) to force linear extrapolation
                $in = $in < 0 ? 0 : ($in > 1 ? 1 : $in);
                 
                # compute parameter differences
                my $diff = [map {$par->[$_ + 1] - $par->[$_]} (0 .. $#{$par} - 1)];
                 
                # if Lapack module loaded
                if ($lapack) {
                         
                        # return Bernstein interpolated derivative
                        return($#{$par} * ICC::Support::Lapack::bernstein($in, $diff));
                         
                } else {
                         
                        # return Bernstein interpolated derivative
                        return($#{$par} * ICC::Shared::dotProduct(_bernstein($#{$diff}, $in), $diff));
                         
                }
                 
        }
         
}
# compute Bernstein basis values
# parameters: (bernstein_degree, t)
# returns: (ref_to_bernstein_array)
sub _bernstein {
        # get parameters
        my ($degree, $t0) = @_;
        # local variables
        my ($bern, $t1, $m0, $m1, $n);
        # binomial coefficients
        state $pascal = [
                [1],
                [1, 1],
                [1, 2, 1],
                [1, 3, 3, 1],
                [1, 4, 6, 4, 1],
                [1, 5, 10, 10, 5, 1],
                [1, 6, 15, 20, 15, 6, 1],
                [1, 7, 21, 35, 35, 21, 7, 1],
                [1, 8, 28, 56, 70, 56, 28, 8, 1],
                [1, 9, 36, 84, 126, 126, 84, 36, 9, 1],
                [1, 10, 45, 120, 210, 252, 210, 120, 45, 10, 1],
                [1, 11, 55, 165, 330, 462, 462, 330, 165, 55, 11, 1],
                [1, 12, 66, 220, 495, 792, 924, 792, 495, 220, 66, 12, 1],
                [1, 13, 78, 286, 715, 1287, 1716, 1716, 1287, 715, 286, 78, 13, 1],
                [1, 14, 91, 364, 1001, 2002, 3003, 3432, 3003, 2002, 1001, 364, 91, 14, 1],
                [1, 15, 105, 455, 1365, 3003, 5005, 6435, 6435, 5005, 3003, 1365, 455, 105, 15, 1],
                [1, 16, 120, 560, 1820, 4368, 8008, 11440, 12870, 11440, 8008, 4368, 1820, 560, 120, 16, 1],
                [1, 17, 136, 680, 2380, 6188, 12376, 19448, 24310, 24310, 19448, 12376, 6188, 2380, 680, 136, 17, 1],
                [1, 18, 153, 816, 3060, 8568, 18564, 31824, 43758, 48620, 43758, 31824, 18564, 8568, 3060, 816, 153, 18, 1],
                [1, 19, 171, 969, 3876, 11628, 27132, 50388, 75582, 92378, 92378, 75582, 50388, 27132, 11628, 3876, 969, 171, 19, 1],
                [1, 20, 190, 1140, 4845, 15504, 38760, 77520, 125970, 167960, 184756, 167960, 125970, 77520, 38760, 15504, 4845, 1140, 190, 20, 1]
        ];
        # copy binomial coefficients for this degree
        $bern = [@{$pascal->[$degree]}];
        # compute input complement
        $t1 = 1 - $t0;
        # initialize multipliers
        $m0 = $m1 = 1;
        # get polynomial degree
        $n = $#{$bern};
        # for each polynomial
        for (1 .. $n) {
                 
                # multiply
                $bern->[$_] *= ($m0 *= $t0);
                $bern->[$n - $_] *= ($m1 *= $t1);
                 
        }
        # return array reference
        return($bern);
}
# fit bern object to data
# uses LAPACK dgels function to perform a least-squares fit
# the parameters to be fitted have the value 'undef' and must all be within the
# output side of the parameter array, e.g. [[0, 100], [undef, undef, undef, undef]]
# fit data is an array containing input and output vectors
# parameters: (ref_to_object, fit_data_reference, ref_to_parameter_hash)
# returns: (dgels_info_value)
sub _fit {
        # get parameters
        my ($self, $fit, $hash) = @_;
        # local variables
        my ($m, $n, $d, $t, $p, $span, @sol);
        my (@so, $outz, $bern, $nrhs, $info);
        my ($outz2, $bern2, $fix_hl, $fix_sh, @c2);
        # check if ICC::Support::Lapack module is loaded
        state $lapack = defined($INC{'ICC/Support/Lapack.pm'});
        # verify Lapack module loaded (need to add non-Lapack routine someday)
        ($lapack) or croak("'fit' option requires ICC::Support::Lapack module");
        # get degree
        $d = $#{$self->[1][1]};
        # for each output parameter
        for my $i (0 .. $d) {
                 
                # if parameter is undefined
                if (! defined($self->[1][1][$i])) {
                         
                        # push index
                        push(@so, $i);
                         
                        # set value to 0
                        $self->[1][1][$i] = 0;
                         
                }
                 
        }
        # if no 'undef' parameters
        if (@so == 0) {
                 
                # warning
                print "no 'undef' parameters to fit\n";
                 
                # return
                return(0)
                 
        }
        # check for min/max values
        _minmax($self);
        # for each sample
        for my $i (0 .. $#{$fit->[0]}) {
                 
                # compute output difference
                $outz->[$i] = [$fit->[1][$i] - transform($self, $fit->[0][$i])];
                 
                # get intermediate value
                $t = _rev($self->[1][0], $self->[2][0], $self->[3][0], $fit->[0][$i]);
                 
                # compute Bernstein values for undefined parameters
                $bern->[$i] = [@{_bernstein($d, $t)}[@so]];
                 
        }
        # make copies
        $outz2 = Storable::dclone($outz);
        $bern2 = Storable::dclone($bern);
        # get array sizes
        $m = @{$bern};
        $n = @{$bern->[0]};
        $nrhs = @{$outz->[0]};
        # fit the data, return if failed
        ($info = ICC::Support::Lapack::dgels('N', $m, $n, $nrhs, $bern, $m, $outz, $m)) && return($info);
        # get curve polarity (0 - increasing, 1 - decreasing)
        $p = ($span = ($so[-1] == $d ? $outz->[$n - 1][0] : $self->[1][1][-1]) - ($so[0] == 0 ? $outz->[0][0] : $self->[1][1][0])) < 0;
        # get flags (indicating HL or SH contrast needs to be fixed)
        $fix_hl = $hash->{'fix_hl'} && $so[0] == 0 && $so[1] == 1 && (($outz->[0][0] > $outz->[1][0]) ^ $p);
        $fix_sh = $hash->{'fix_sh'} && $so[-1] == $d && $so[-2] == $d - 1 && (($outz->[$n - 1][0] < $outz->[$n - 2][0]) ^ $p);
        # if degree >= 4 and HL or SH needs fixing
        if ($d >= 4 && ($fix_hl || $fix_sh)) {
                 
                # for each sample
                for my $i (0 .. $#{$bern2}) {
                         
                        # if highlight fix
                        if ($fix_hl) {
                                 
                                # get first 2 Bernstein values
                                @c2 = splice(@{$bern2->[$i]}, 0, 2);
                                 
                                # combine to single value
                                unshift(@{$bern2->[$i]}, ($c2[0] + $c2[1]));
                                 
                        }
                         
                        # if shadow fix
                        if ($fix_sh) {
                                 
                                # get last 2 Bernstein values
                                @c2 = splice(@{$bern2->[$i]}, -2);
                                 
                                # combine to single value
                                push(@{$bern2->[$i]}, ($c2[0] + $c2[1]));
                                 
                        }
                         
                }
                 
                # get array sizes
                $m = @{$bern2};
                $n = @{$bern2->[0]};
                $nrhs = @{$outz2->[0]};
                 
                # fit the data, return if failed
                ($info = ICC::Support::Lapack::dgels('N', $m, $n, $nrhs, $bern2, $m, $outz2, $m)) && return($info);
                 
                # extract solution
                @sol = map {$outz2->[$_][0]} (0 .. $n - 1);
                 
                # replace combined coefficient(s)
                splice(@sol, 1, 0, $sol[0] + $span * 0.001) if ($fix_hl);
                splice(@sol, -1, 0, $sol[-1] - $span * 0.001) if ($fix_sh);
                 
        } else {
                 
                # extract solution
                @sol = map {$outz->[$_][0]} (0 .. $n - 1);
                 
        }
        # copy solution to output array
        @{$self->[1][1]}[@so] = @sol;
        # return
        return(0);
}
# set object contents from parameter hash
# parameters: (ref_to_object, ref_to_parameter_hash)
sub _new_from_hash {
        # get parameters
        my ($self, $hash) = @_;
        # local variables
        my ($in, $out, $fit);
        # if 'input' is defined
        if (defined($in = $hash->{'input'})) {
                 
                # if 'input' is a vector (21 values or less, some may be 'undef')
                if (ICC::Shared::is_vector($in) && 21 >= @{$in}) {
                         
                        # set object 'input' array
                        $self->[1][0] = [@{$in}];
                         
                # if 'input' is an integer between 1 and 20
                } elsif (Scalar::Util::looks_like_number($in) && $in == int($in) && $in > 0 && $in < 21) {
                         
                        # make identity curve
                        $self->[1][0] = [map {$_/$in} (0 .. $in)];
                         
                } else {
                         
                        # error
                        croak("invalid 'input' values");
                         
                }
                 
        } else {
                 
                # disable input
                $self->[1][0] = [];
                 
        }
        # if 'output' is defined
        if (defined($out = $hash->{'output'})) {
                 
                # if 'output' is a vector (21 values or less, some may be 'undef')
                if (ICC::Shared::is_vector($out) && 21 >= @{$out}) {
                         
                        # set object 'output' array
                        $self->[1][1] = [@{$out}];
                         
                # if 'output' is an integer between 1 and 20
                } elsif (Scalar::Util::looks_like_number($out) && $out == int($out) && $out > 0 && $out < 21) {
                         
                        # make identity curve
                        $self->[1][1] = [map {$_/$out} (0 .. $out)];
                         
                } else {
                         
                        # error
                        croak("invalid 'output' values");
                         
                }
                 
        } else {
                 
                # disable output
                $self->[1][1] = [];
                 
        }
        # if 'fit' is defined and not empty
        if (defined($fit = $hash->{'fit'}) && @{$fit}) {
                 
                # verify 'fit' vectors
                (ICC::Shared::is_num_vector($fit->[0]) && 1 <= @{$fit->[0]}) or croak("invalid 'fit' input values");
                (ICC::Shared::is_num_vector($fit->[1]) && 1 <= @{$fit->[1]}) or croak("invalid 'fit' output values");
                (@{$fit->[0]} == @{$fit->[1]}) or croak("'fit' input and output vectors are different size");
                 
                # fit object to data
                (_fit($self, $fit, $hash)) && croak("'fit' operation failed");
                 
        } else {
                 
                # verify object parameters are all defined
                (@{$self->[1][0]} == grep {defined($_)} @{$self->[1][0]}) or croak("some 'input' values are 'undef'");
                (@{$self->[1][1]} == grep {defined($_)} @{$self->[1][1]}) or croak("some 'output' values are 'undef'");
                 
        }
        # add segment min/max y-values
        _minmax($self);
}
# write bern object to ICC profile
# note: writes an equivalent curv object
# parameters: (ref_to_object, ref_to_parent_object, file_handle, ref_to_tag_table_entry)
sub _writeICCbern {
        # get parameters
        my ($self, $parent, $fh, $tag) = @_;
        # local variables
        my ($n, $up, @table);
        # get count value (defaults to 4096, the maximum allowed in an 'mft2' tag)
        $n = $self->[0]{'curv_points'} // 4096;
        # validate number of table entries
        ($n == int($n) && $n >= 2) or carp('invalid number of table entries');
        # array upper index
        $up = $n - 1;
        # for each table entry
        for my $i (0 .. $up) {
                 
                # transform value
                $table[$i] = transform($self, $i/$up);
                 
        }
        # seek start of tag
        seek($fh, $tag->[1], 0);
        # write tag type signature and count
        print $fh pack('a4 x4 N', 'curv', $n);
        # write array
        print $fh pack('n*', map {$_ < 0 ? 0 : ($_ > 1 ? 65535 : $_ * 65535 + 0.5)} @table);
}
1;